<s> """<STR_LIT>""" <EOL> from __future__ import ( absolute_import , division , print_function , <EOL> unicode_literals ) <EOL> import math <EOL> from math import cos , sin <EOL> from matplotlib . patches import Ellipse <EOL> import matplotlib . pyplot as plt <EOL> import numpy as np <EOL> import random <EOL> import scipy . linalg as linalg <EOL> import scipy . sparse as sp <EOL> import scipy . sparse . linalg as spln <EOL> import scipy . stats <EOL> from scipy . stats import norm <EOL> import warnings <EOL> def gaussian ( x , mean , var ) : <EOL> """<STR_LIT>""" <EOL> return ( np . exp ( ( - <NUM_LIT:0.5> * ( np . asarray ( x ) - mean ) ** <NUM_LIT:2> ) / var ) / <EOL> math . sqrt ( <NUM_LIT:2> * math . pi * var ) ) <EOL> def mul ( mean1 , var1 , mean2 , var2 ) : <EOL> """<STR_LIT>""" <EOL> mean = ( var1 * mean2 + var2 * mean1 ) / ( var1 + var2 ) <EOL> var = <NUM_LIT:1> / ( <NUM_LIT:1> / var1 + <NUM_LIT:1> / var2 ) <EOL> return ( mean , var ) <EOL> def add ( mean1 , var1 , mean2 , var2 ) : <EOL> """<STR_LIT>""" <EOL> return ( mean1 + mean2 , var1 + var2 ) <EOL> def multivariate_gaussian ( x , mu , cov ) : <EOL> """<STR_LIT>""" <EOL> x = np . array ( x , copy = False , ndmin = <NUM_LIT:1> ) . flatten ( ) <EOL> mu = np . array ( mu , copy = False , ndmin = <NUM_LIT:1> ) . flatten ( ) <EOL> nx = len ( mu ) <EOL> cov = _to_cov ( cov , nx ) <EOL> norm_coeff = nx * math . log ( <NUM_LIT:2> * math . pi ) + np . linalg . slogdet ( cov ) [ <NUM_LIT:1> ] <EOL> err = x - mu <EOL> if ( sp . issparse ( cov ) ) : <EOL> numerator = spln . spsolve ( cov , err ) . T . dot ( err ) <EOL> else : <EOL> numerator = np . linalg . solve ( cov , err ) . T . dot ( err ) <EOL> return math . exp ( - <NUM_LIT:0.5> * ( norm_coeff + numerator ) ) <EOL> def multivariate_multiply ( m1 , c1 , m2 , c2 ) : <EOL> """<STR_LIT>""" <EOL> C1 = np . asarray ( c1 ) <EOL> C2 = np . asarray ( c2 ) <EOL> M1 = np . asarray ( m1 ) <EOL> M2 = np . asarray ( m2 ) <EOL> sum_inv = np . linalg . inv ( C1 + C2 ) <EOL> C3 = np . dot ( C1 , sum_inv ) . dot ( C2 ) <EOL> M3 = ( np . dot ( C2 , sum_inv ) . dot ( M1 ) + <EOL> np . dot ( C1 , sum_inv ) . dot ( M2 ) ) <EOL> return M3 , C3 <EOL> def plot_discrete_cdf ( xs , ys , ax = None , xlabel = None , ylabel = None , <EOL> label = None ) : <EOL> """<STR_LIT>""" <EOL> if ax is None : <EOL> ax = plt . gca ( ) <EOL> if xs is None : <EOL> xs = range ( len ( ys ) ) <EOL> ys = np . cumsum ( ys ) <EOL> ax . plot ( xs , ys , label = label ) <EOL> ax . set_xlabel ( xlabel ) <EOL> ax . set_ylabel ( ylabel ) <EOL> return ax <EOL> def plot_gaussian_cdf ( mean = <NUM_LIT:0.> , variance = <NUM_LIT:1.> , <EOL> ax = None , <EOL> xlim = None , ylim = ( <NUM_LIT:0.> , <NUM_LIT:1.> ) , <EOL> xlabel = None , ylabel = None , <EOL> label = None ) : <EOL> """<STR_LIT>""" <EOL> if ax is None : <EOL> ax = plt . gca ( ) <EOL> sigma = math . sqrt ( variance ) <EOL> n = scipy . stats . norm ( mean , sigma ) <EOL> if xlim is None : <EOL> xlim = [ n . ppf ( <NUM_LIT> ) , n . ppf ( <NUM_LIT> ) ] <EOL> xs = np . arange ( xlim [ <NUM_LIT:0> ] , xlim [ <NUM_LIT:1> ] , ( xlim [ <NUM_LIT:1> ] - xlim [ <NUM_LIT:0> ] ) / <NUM_LIT> ) <EOL> cdf = n . cdf ( xs ) <EOL> ax . plot ( xs , cdf , label = label ) <EOL> ax . set_xlim ( xlim ) <EOL> ax . set_ylim ( ylim ) <EOL> ax . set_xlabel ( xlabel ) <EOL> ax . set_ylabel ( ylabel ) <EOL> return ax <EOL> def plot_gaussian_pdf ( mean = <NUM_LIT:0.> , variance = <NUM_LIT:1.> , <EOL> ax = None , <EOL> mean_line = False , <EOL> xlim = None , ylim = None , <EOL> xlabel = None , ylabel = None , <EOL> label = None ) : <EOL> """<STR_LIT>""" <EOL> if ax is None : <EOL> ax = plt . gca ( ) <EOL> sigma = math . sqrt ( variance ) <EOL> n = scipy . stats . norm ( mean , sigma ) <EOL> if xlim is None : <EOL> xlim = [ n . ppf ( <NUM_LIT> ) , n . ppf ( <NUM_LIT> ) ] <EOL> xs = np . arange ( xlim [ <NUM_LIT:0> ] , xlim [ <NUM_LIT:1> ] , ( xlim [ <NUM_LIT:1> ] - xlim [ <NUM_LIT:0> ] ) / <NUM_LIT> ) <EOL> ax . plot ( xs , n . pdf ( xs ) , label = label ) <EOL> ax . set_xlim ( xlim ) <EOL> if ylim is not None : <EOL> ax . set_ylim ( ylim ) <EOL> if mean_line : <EOL> plt . axvline ( mean ) <EOL> if xlabel is not None : <EOL> ax . set_xlabel ( xlabel ) <EOL> if ylabel is not None : <EOL> ax . set_ylabel ( ylabel ) <EOL> return ax <EOL> def plot_gaussian ( mean = <NUM_LIT:0.> , variance = <NUM_LIT:1.> , <EOL> ax = None , <EOL> mean_line = False , <EOL> xlim = None , <EOL> ylim = None , <EOL> xlabel = None , <EOL> ylabel = None , <EOL> label = None ) : <EOL> """<STR_LIT>""" <EOL> warnings . warn ( '<STR_LIT>' <EOL> '<STR_LIT>' <EOL> '<STR_LIT>' , <EOL> DeprecationWarning ) <EOL> return plot_gaussian_pdf ( mean , variance , ax , mean_line , xlim , ylim , xlabel , <EOL> ylabel , label ) <EOL> def covariance_ellipse ( P , deviations = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> U , s , v = linalg . svd ( P ) <EOL> orientation = math . atan2 ( U [ <NUM_LIT:1> , <NUM_LIT:0> ] , U [ <NUM_LIT:0> , <NUM_LIT:0> ] ) <EOL> width = deviations * math . sqrt ( s [ <NUM_LIT:0> ] ) <EOL> height = deviations * math . sqrt ( s [ <NUM_LIT:1> ] ) <EOL> assert width >= height <EOL> return ( orientation , width , height ) <EOL> def plot_covariance_ellipse ( mean , cov = None , variance = <NUM_LIT:1.0> , std = None , <EOL> ellipse = None , title = None , axis_equal = True , show_semiaxis = False , <EOL> facecolor = None , edgecolor = None , <EOL> fc = '<STR_LIT:none>' , ec = '<STR_LIT>' , <EOL> alpha = <NUM_LIT:1.0> , xlim = None , ylim = None , <EOL> ls = '<STR_LIT>' ) : <EOL> """<STR_LIT>""" <EOL> assert cov is None or ellipse is None <EOL> assert not ( cov is None and ellipse is None ) <EOL> if facecolor is None : <EOL> facecolor = fc <EOL> if edgecolor is None : <EOL> edgecolor = ec <EOL> if cov is not None : <EOL> ellipse = covariance_ellipse ( cov ) <EOL> if axis_equal : <EOL> plt . axis ( '<STR_LIT>' ) <EOL> if title is not None : <EOL> plt . title ( title ) <EOL> compute_std = False <EOL> if std is None : <EOL> std = variance <EOL> compute_std = True <EOL> if np . isscalar ( std ) : <EOL> std = [ std ] <EOL> if compute_std : <EOL> std = np . sqrt ( np . asarray ( std ) ) <EOL> ax = plt . gca ( ) <EOL> angle = np . degrees ( ellipse [ <NUM_LIT:0> ] ) <EOL> width = ellipse [ <NUM_LIT:1> ] * <NUM_LIT> <EOL> height = ellipse [ <NUM_LIT:2> ] * <NUM_LIT> <EOL> for sd in std : <EOL> e = Ellipse ( xy = mean , width = sd * width , height = sd * height , angle = angle , <EOL> facecolor = facecolor , <EOL> edgecolor = edgecolor , <EOL> alpha = alpha , <EOL> lw = <NUM_LIT:2> , ls = ls ) <EOL> ax . add_patch ( e ) <EOL> x , y = mean <EOL> plt . scatter ( x , y , marker = '<STR_LIT:+>' , color = edgecolor ) <EOL> if xlim is not None : <EOL> ax . set_xlim ( xlim ) <EOL> if ylim is not None : <EOL> ax . set_ylim ( ylim ) <EOL> if show_semiaxis : <EOL> a = ellipse [ <NUM_LIT:0> ] <EOL> h , w = height / <NUM_LIT:4> , width / <NUM_LIT:4> <EOL> plt . plot ( [ x , x + h * cos ( a + np . pi / <NUM_LIT:2> ) ] , [ y , y + h * sin ( a + np . pi / <NUM_LIT:2> ) ] ) <EOL> plt . plot ( [ x , x + w * cos ( a ) ] , [ y , y + w * sin ( a ) ] ) <EOL> def norm_cdf ( x_range , mu , var = <NUM_LIT:1> , std = None ) : <EOL> """<STR_LIT>""" <EOL> if std is None : <EOL> std = math . sqrt ( var ) <EOL> return abs ( norm . cdf ( x_range [ <NUM_LIT:0> ] , loc = mu , scale = std ) - <EOL> norm . cdf ( x_range [ <NUM_LIT:1> ] , loc = mu , scale = std ) ) <EOL> def _is_inside_ellipse ( x , y , ex , ey , orientation , width , height ) : <EOL> co = np . cos ( orientation ) <EOL> so = np . sin ( orientation ) <EOL> xx = x * co + y * so <EOL> yy = y * co - x * so <EOL> return ( xx / width ) ** <NUM_LIT:2> + ( yy / height ) ** <NUM_LIT:2> <= <NUM_LIT:1.> <EOL> return ( ( x - ex ) * co - ( y - ey ) * so ) ** <NUM_LIT:2> / width ** <NUM_LIT:2> + ( ( x - ex ) * so + ( y - ey ) * co ) ** <NUM_LIT:2> / height ** <NUM_LIT:2> <= <NUM_LIT:1> <EOL> def _to_cov ( x , n ) : <EOL> """<STR_LIT>""" <EOL> try : <EOL> x . shape <EOL> if type ( x ) != np . ndarray : <EOL> x = np . asarray ( x ) [ <NUM_LIT:0> ] <EOL> return x <EOL> except : <EOL> return np . eye ( n ) * x <EOL> def _do_plot_test ( ) : <EOL> from numpy . random import multivariate_normal <EOL> p = np . array ( [ [ <NUM_LIT:32> , <NUM_LIT:15> ] , [ <NUM_LIT> , <NUM_LIT> ] ] ) <EOL> x , y = multivariate_normal ( mean = ( <NUM_LIT:0> , <NUM_LIT:0> ) , cov = p , size = <NUM_LIT> ) . T <EOL> sd = <NUM_LIT:2> <EOL> a , w , h = covariance_ellipse ( p , sd ) <EOL> print ( np . degrees ( a ) , w , h ) <EOL> count = <NUM_LIT:0> <EOL> color = [ ] <EOL> for i in range ( len ( x ) ) : <EOL> if _is_inside_ellipse ( x [ i ] , y [ i ] , <NUM_LIT:0> , <NUM_LIT:0> , a , w , h ) : <EOL> color . append ( '<STR_LIT:b>' ) <EOL> count += <NUM_LIT:1> <EOL> else : <EOL> color . append ( '<STR_LIT:r>' ) <EOL> plt . scatter ( x , y , alpha = <NUM_LIT> , c = color ) <EOL> plt . axis ( '<STR_LIT>' ) <EOL> plot_covariance_ellipse ( mean = ( <NUM_LIT:0.> , <NUM_LIT:0.> ) , <EOL> cov = p , <EOL> std = sd , <EOL> facecolor = '<STR_LIT:none>' ) <EOL> print ( count / len ( x ) ) <EOL> def plot_std_vs_var ( ) : <EOL> plt . figure ( ) <EOL> x = ( <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> P = np . array ( [ [ <NUM_LIT:3> , <NUM_LIT:1> ] , [ <NUM_LIT:1> , <NUM_LIT:3> ] ] ) <EOL> plot_covariance_ellipse ( x , P , std = [ <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:3> ] , facecolor = '<STR_LIT:g>' , alpha = <NUM_LIT> ) <EOL> plot_covariance_ellipse ( x , P , variance = [ <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:3> ] , facecolor = '<STR_LIT:r>' , alpha = <NUM_LIT> ) <EOL> def rand_student_t ( df , mu = <NUM_LIT:0> , std = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> x = random . gauss ( <NUM_LIT:0> , std ) <EOL> y = <NUM_LIT> * random . gammavariate ( <NUM_LIT:0.5> * df , <NUM_LIT> ) <EOL> return x / ( math . sqrt ( y / df ) ) + mu <EOL> def NESS ( xs , est_xs , ps ) : <EOL> """<STR_LIT>""" <EOL> est_err = xs - est_xs <EOL> ness = [ ] <EOL> for x , p in zip ( est_err , ps ) : <EOL> ness . append ( np . dot ( x . T , linalg . inv ( p ) ) . dot ( x ) ) <EOL> return ness <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> ax = plot_gaussian_pdf ( <NUM_LIT:2> , <NUM_LIT:3> ) <EOL> plot_gaussian_cdf ( <NUM_LIT:2> , <NUM_LIT:3> , ax = ax ) <EOL> plt . show ( ) <EOL> ys = np . abs ( np . random . randn ( <NUM_LIT:100> ) ) <EOL> ys /= np . sum ( ys ) <EOL> plot_discrete_cdf ( xs = None , ys = ys ) <EOL> """<STR_LIT>""" </s>
