<s> from __future__ import division , print_function , absolute_import <EOL> import math <EOL> import numpy as np <EOL> import scipy . linalg <EOL> from scipy . misc import doccer <EOL> from scipy . special import gammaln , psi , multigammaln <EOL> from scipy . _lib . _util import check_random_state <EOL> from scipy . linalg . blas import drot <EOL> __all__ = [ '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' , <EOL> '<STR_LIT>' ] <EOL> _LOG_2PI = np . log ( <NUM_LIT:2> * np . pi ) <EOL> _LOG_2 = np . log ( <NUM_LIT:2> ) <EOL> _LOG_PI = np . log ( np . pi ) <EOL> _doc_random_state = """<STR_LIT>""" <EOL> def _squeeze_output ( out ) : <EOL> """<STR_LIT>""" <EOL> out = out . squeeze ( ) <EOL> if out . ndim == <NUM_LIT:0> : <EOL> out = out [ ( ) ] <EOL> return out <EOL> def _eigvalsh_to_eps ( spectrum , cond = None , rcond = None ) : <EOL> """<STR_LIT>""" <EOL> if rcond is not None : <EOL> cond = rcond <EOL> if cond in [ None , - <NUM_LIT:1> ] : <EOL> t = spectrum . dtype . char . lower ( ) <EOL> factor = { '<STR_LIT:f>' : <NUM_LIT> , '<STR_LIT:d>' : <NUM_LIT> } <EOL> cond = factor [ t ] * np . finfo ( t ) . eps <EOL> eps = cond * np . max ( abs ( spectrum ) ) <EOL> return eps <EOL> def _pinv_1d ( v , eps = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> return np . array ( [ <NUM_LIT:0> if abs ( x ) <= eps else <NUM_LIT:1> / x for x in v ] , dtype = float ) <EOL> class _PSD ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , M , cond = None , rcond = None , lower = True , <EOL> check_finite = True , allow_singular = True ) : <EOL> s , u = scipy . linalg . eigh ( M , lower = lower , check_finite = check_finite ) <EOL> eps = _eigvalsh_to_eps ( s , cond , rcond ) <EOL> if np . min ( s ) < - eps : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> d = s [ s > eps ] <EOL> if len ( d ) < len ( s ) and not allow_singular : <EOL> raise np . linalg . LinAlgError ( '<STR_LIT>' ) <EOL> s_pinv = _pinv_1d ( s , eps ) <EOL> U = np . multiply ( u , np . sqrt ( s_pinv ) ) <EOL> self . rank = len ( d ) <EOL> self . U = U <EOL> self . log_pdet = np . sum ( np . log ( d ) ) <EOL> self . _pinv = None <EOL> @ property <EOL> def pinv ( self ) : <EOL> if self . _pinv is None : <EOL> self . _pinv = np . dot ( self . U , self . U . T ) <EOL> return self . _pinv <EOL> class multi_rv_generic ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( multi_rv_generic , self ) . __init__ ( ) <EOL> self . _random_state = check_random_state ( seed ) <EOL> @ property <EOL> def random_state ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . _random_state <EOL> @ random_state . setter <EOL> def random_state ( self , seed ) : <EOL> self . _random_state = check_random_state ( seed ) <EOL> def _get_random_state ( self , random_state ) : <EOL> if random_state is not None : <EOL> return check_random_state ( random_state ) <EOL> else : <EOL> return self . _random_state <EOL> class multi_rv_frozen ( object ) : <EOL> """<STR_LIT>""" <EOL> @ property <EOL> def random_state ( self ) : <EOL> return self . _dist . _random_state <EOL> @ random_state . setter <EOL> def random_state ( self , seed ) : <EOL> self . _dist . _random_state = check_random_state ( seed ) <EOL> _mvn_doc_default_callparams = """<STR_LIT>""" <EOL> _mvn_doc_callparams_note = """<STR_LIT>""" <EOL> _mvn_doc_frozen_callparams = "<STR_LIT>" <EOL> _mvn_doc_frozen_callparams_note = """<STR_LIT>""" <EOL> mvn_docdict_params = { <EOL> '<STR_LIT>' : _mvn_doc_default_callparams , <EOL> '<STR_LIT>' : _mvn_doc_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> mvn_docdict_noparams = { <EOL> '<STR_LIT>' : _mvn_doc_frozen_callparams , <EOL> '<STR_LIT>' : _mvn_doc_frozen_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> class multivariate_normal_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( multivariate_normal_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ , mvn_docdict_params ) <EOL> def __call__ ( self , mean = None , cov = <NUM_LIT:1> , allow_singular = False , seed = None ) : <EOL> """<STR_LIT>""" <EOL> return multivariate_normal_frozen ( mean , cov , <EOL> allow_singular = allow_singular , <EOL> seed = seed ) <EOL> def _process_parameters ( self , dim , mean , cov ) : <EOL> """<STR_LIT>""" <EOL> if dim is None : <EOL> if mean is None : <EOL> if cov is None : <EOL> dim = <NUM_LIT:1> <EOL> else : <EOL> cov = np . asarray ( cov , dtype = float ) <EOL> if cov . ndim < <NUM_LIT:2> : <EOL> dim = <NUM_LIT:1> <EOL> else : <EOL> dim = cov . shape [ <NUM_LIT:0> ] <EOL> else : <EOL> mean = np . asarray ( mean , dtype = float ) <EOL> dim = mean . size <EOL> else : <EOL> if not np . isscalar ( dim ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if mean is None : <EOL> mean = np . zeros ( dim ) <EOL> mean = np . asarray ( mean , dtype = float ) <EOL> if cov is None : <EOL> cov = <NUM_LIT:1.0> <EOL> cov = np . asarray ( cov , dtype = float ) <EOL> if dim == <NUM_LIT:1> : <EOL> mean . shape = ( <NUM_LIT:1> , ) <EOL> cov . shape = ( <NUM_LIT:1> , <NUM_LIT:1> ) <EOL> if mean . ndim != <NUM_LIT:1> or mean . shape [ <NUM_LIT:0> ] != dim : <EOL> raise ValueError ( "<STR_LIT>" % dim ) <EOL> if cov . ndim == <NUM_LIT:0> : <EOL> cov = cov * np . eye ( dim ) <EOL> elif cov . ndim == <NUM_LIT:1> : <EOL> cov = np . diag ( cov ) <EOL> elif cov . ndim == <NUM_LIT:2> and cov . shape != ( dim , dim ) : <EOL> rows , cols = cov . shape <EOL> if rows != cols : <EOL> msg = ( "<STR_LIT>" <EOL> "<STR_LIT>" % str ( cov . shape ) ) <EOL> else : <EOL> msg = ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> msg = msg % ( str ( cov . shape ) , len ( mean ) ) <EOL> raise ValueError ( msg ) <EOL> elif cov . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % cov . ndim ) <EOL> return dim , mean , cov <EOL> def _process_quantiles ( self , x , dim ) : <EOL> """<STR_LIT>""" <EOL> x = np . asarray ( x , dtype = float ) <EOL> if x . ndim == <NUM_LIT:0> : <EOL> x = x [ np . newaxis ] <EOL> elif x . ndim == <NUM_LIT:1> : <EOL> if dim == <NUM_LIT:1> : <EOL> x = x [ : , np . newaxis ] <EOL> else : <EOL> x = x [ np . newaxis , : ] <EOL> return x <EOL> def _logpdf ( self , x , mean , prec_U , log_det_cov , rank ) : <EOL> """<STR_LIT>""" <EOL> dev = x - mean <EOL> maha = np . sum ( np . square ( np . dot ( dev , prec_U ) ) , axis = - <NUM_LIT:1> ) <EOL> return - <NUM_LIT:0.5> * ( rank * _LOG_2PI + log_det_cov + maha ) <EOL> def logpdf ( self , x , mean , cov , allow_singular = False ) : <EOL> """<STR_LIT>""" <EOL> dim , mean , cov = self . _process_parameters ( None , mean , cov ) <EOL> x = self . _process_quantiles ( x , dim ) <EOL> psd = _PSD ( cov , allow_singular = allow_singular ) <EOL> out = self . _logpdf ( x , mean , psd . U , psd . log_pdet , psd . rank ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x , mean , cov , allow_singular = False ) : <EOL> """<STR_LIT>""" <EOL> dim , mean , cov = self . _process_parameters ( None , mean , cov ) <EOL> x = self . _process_quantiles ( x , dim ) <EOL> psd = _PSD ( cov , allow_singular = allow_singular ) <EOL> out = np . exp ( self . _logpdf ( x , mean , psd . U , psd . log_pdet , psd . rank ) ) <EOL> return _squeeze_output ( out ) <EOL> def rvs ( self , mean = None , cov = <NUM_LIT:1> , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> dim , mean , cov = self . _process_parameters ( None , mean , cov ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> out = random_state . multivariate_normal ( mean , cov , size ) <EOL> return _squeeze_output ( out ) <EOL> def entropy ( self , mean = None , cov = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> dim , mean , cov = self . _process_parameters ( None , mean , cov ) <EOL> _ , logdet = np . linalg . slogdet ( <NUM_LIT:2> * np . pi * np . e * cov ) <EOL> return <NUM_LIT:0.5> * logdet <EOL> multivariate_normal = multivariate_normal_gen ( ) <EOL> class multivariate_normal_frozen ( multi_rv_frozen ) : <EOL> def __init__ ( self , mean = None , cov = <NUM_LIT:1> , allow_singular = False , seed = None ) : <EOL> """<STR_LIT>""" <EOL> self . _dist = multivariate_normal_gen ( seed ) <EOL> self . dim , self . mean , self . cov = self . _dist . _process_parameters ( <EOL> None , mean , cov ) <EOL> self . cov_info = _PSD ( self . cov , allow_singular = allow_singular ) <EOL> def logpdf ( self , x ) : <EOL> x = self . _dist . _process_quantiles ( x , self . dim ) <EOL> out = self . _dist . _logpdf ( x , self . mean , self . cov_info . U , <EOL> self . cov_info . log_pdet , self . cov_info . rank ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x ) : <EOL> return np . exp ( self . logpdf ( x ) ) <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> return self . _dist . rvs ( self . mean , self . cov , size , random_state ) <EOL> def entropy ( self ) : <EOL> """<STR_LIT>""" <EOL> log_pdet = self . cov_info . log_pdet <EOL> rank = self . cov_info . rank <EOL> return <NUM_LIT:0.5> * ( rank * ( _LOG_2PI + <NUM_LIT:1> ) + log_pdet ) <EOL> for name in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> method = multivariate_normal_gen . __dict__ [ name ] <EOL> method_frozen = multivariate_normal_frozen . __dict__ [ name ] <EOL> method_frozen . __doc__ = doccer . docformat ( method . __doc__ , mvn_docdict_noparams ) <EOL> method . __doc__ = doccer . docformat ( method . __doc__ , mvn_docdict_params ) <EOL> _matnorm_doc_default_callparams = """<STR_LIT>""" <EOL> _matnorm_doc_callparams_note = """<STR_LIT>""" <EOL> _matnorm_doc_frozen_callparams = "<STR_LIT>" <EOL> _matnorm_doc_frozen_callparams_note = """<STR_LIT>""" <EOL> matnorm_docdict_params = { <EOL> '<STR_LIT>' : _matnorm_doc_default_callparams , <EOL> '<STR_LIT>' : _matnorm_doc_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> matnorm_docdict_noparams = { <EOL> '<STR_LIT>' : _matnorm_doc_frozen_callparams , <EOL> '<STR_LIT>' : _matnorm_doc_frozen_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> class matrix_normal_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( matrix_normal_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ , matnorm_docdict_params ) <EOL> def __call__ ( self , mean = None , rowcov = <NUM_LIT:1> , colcov = <NUM_LIT:1> , seed = None ) : <EOL> """<STR_LIT>""" <EOL> return matrix_normal_frozen ( mean , rowcov , colcov , seed = seed ) <EOL> def _process_parameters ( self , mean , rowcov , colcov ) : <EOL> """<STR_LIT>""" <EOL> if mean is not None : <EOL> mean = np . asarray ( mean , dtype = float ) <EOL> meanshape = mean . shape <EOL> if len ( meanshape ) != <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if np . any ( meanshape == <NUM_LIT:0> ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> rowcov = np . asarray ( rowcov , dtype = float ) <EOL> if rowcov . ndim == <NUM_LIT:0> : <EOL> if mean is not None : <EOL> rowcov = rowcov * np . identity ( meanshape [ <NUM_LIT:0> ] ) <EOL> else : <EOL> rowcov = rowcov * np . identity ( <NUM_LIT:1> ) <EOL> elif rowcov . ndim == <NUM_LIT:1> : <EOL> rowcov = np . diag ( rowcov ) <EOL> rowshape = rowcov . shape <EOL> if len ( rowshape ) != <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if rowshape [ <NUM_LIT:0> ] != rowshape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if rowshape [ <NUM_LIT:0> ] == <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> numrows = rowshape [ <NUM_LIT:0> ] <EOL> colcov = np . asarray ( colcov , dtype = float ) <EOL> if colcov . ndim == <NUM_LIT:0> : <EOL> if mean is not None : <EOL> colcov = colcov * np . identity ( meanshape [ <NUM_LIT:1> ] ) <EOL> else : <EOL> colcov = colcov * np . identity ( <NUM_LIT:1> ) <EOL> elif colcov . ndim == <NUM_LIT:1> : <EOL> colcov = np . diag ( colcov ) <EOL> colshape = colcov . shape <EOL> if len ( colshape ) != <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if colshape [ <NUM_LIT:0> ] != colshape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if colshape [ <NUM_LIT:0> ] == <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> numcols = colshape [ <NUM_LIT:0> ] <EOL> if mean is not None : <EOL> if meanshape [ <NUM_LIT:0> ] != numrows : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> if meanshape [ <NUM_LIT:1> ] != numcols : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> else : <EOL> mean = np . zeros ( ( numrows , numcols ) ) <EOL> dims = ( numrows , numcols ) <EOL> return dims , mean , rowcov , colcov <EOL> def _process_quantiles ( self , X , dims ) : <EOL> """<STR_LIT>""" <EOL> X = np . asarray ( X , dtype = float ) <EOL> if X . ndim == <NUM_LIT:2> : <EOL> X = X [ np . newaxis , : ] <EOL> if X . shape [ - <NUM_LIT:2> : ] != dims : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> return X <EOL> def _logpdf ( self , dims , X , mean , row_prec_rt , log_det_rowcov , <EOL> col_prec_rt , log_det_colcov ) : <EOL> """<STR_LIT>""" <EOL> numrows , numcols = dims <EOL> roll_dev = np . rollaxis ( X - mean , axis = - <NUM_LIT:1> , start = <NUM_LIT:0> ) <EOL> scale_dev = np . tensordot ( col_prec_rt . T , <EOL> np . dot ( roll_dev , row_prec_rt ) , <NUM_LIT:1> ) <EOL> maha = np . sum ( np . sum ( np . square ( scale_dev ) , axis = - <NUM_LIT:1> ) , axis = <NUM_LIT:0> ) <EOL> return - <NUM_LIT:0.5> * ( numrows * numcols * _LOG_2PI + numcols * log_det_rowcov <EOL> + numrows * log_det_colcov + maha ) <EOL> def logpdf ( self , X , mean = None , rowcov = <NUM_LIT:1> , colcov = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> dims , mean , rowcov , colcov = self . _process_parameters ( mean , rowcov , <EOL> colcov ) <EOL> X = self . _process_quantiles ( X , dims ) <EOL> rowpsd = _PSD ( rowcov , allow_singular = False ) <EOL> colpsd = _PSD ( colcov , allow_singular = False ) <EOL> out = self . _logpdf ( dims , X , mean , rowpsd . U , rowpsd . log_pdet , colpsd . U , <EOL> colpsd . log_pdet ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , X , mean = None , rowcov = <NUM_LIT:1> , colcov = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> return np . exp ( self . logpdf ( X , mean , rowcov , colcov ) ) <EOL> def rvs ( self , mean = None , rowcov = <NUM_LIT:1> , colcov = <NUM_LIT:1> , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> size = int ( size ) <EOL> dims , mean , rowcov , colcov = self . _process_parameters ( mean , rowcov , <EOL> colcov ) <EOL> rowchol = scipy . linalg . cholesky ( rowcov , lower = True ) <EOL> colchol = scipy . linalg . cholesky ( colcov , lower = True ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> std_norm = random_state . standard_normal ( size = ( dims [ <NUM_LIT:1> ] , size , dims [ <NUM_LIT:0> ] ) ) <EOL> roll_rvs = np . tensordot ( colchol , np . dot ( std_norm , rowchol . T ) , <NUM_LIT:1> ) <EOL> out = np . rollaxis ( roll_rvs . T , axis = <NUM_LIT:1> , start = <NUM_LIT:0> ) + mean [ np . newaxis , : , : ] <EOL> if size == <NUM_LIT:1> : <EOL> out = out . reshape ( mean . shape ) <EOL> return out <EOL> matrix_normal = matrix_normal_gen ( ) <EOL> class matrix_normal_frozen ( multi_rv_frozen ) : <EOL> def __init__ ( self , mean = None , rowcov = <NUM_LIT:1> , colcov = <NUM_LIT:1> , seed = None ) : <EOL> """<STR_LIT>""" <EOL> self . _dist = matrix_normal_gen ( seed ) <EOL> self . dims , self . mean , self . rowcov , self . colcov = self . _dist . _process_parameters ( mean , rowcov , colcov ) <EOL> self . rowpsd = _PSD ( self . rowcov , allow_singular = False ) <EOL> self . colpsd = _PSD ( self . colcov , allow_singular = False ) <EOL> def logpdf ( self , X ) : <EOL> X = self . _dist . _process_quantiles ( X , self . dims ) <EOL> out = self . _dist . _logpdf ( self . dims , X , self . mean , self . rowpsd . U , <EOL> self . rowpsd . log_pdet , self . colpsd . U , <EOL> self . colpsd . log_pdet ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , X ) : <EOL> return np . exp ( self . logpdf ( X ) ) <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> return self . _dist . rvs ( self . mean , self . rowcov , self . colcov , size , <EOL> random_state ) <EOL> for name in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> method = matrix_normal_gen . __dict__ [ name ] <EOL> method_frozen = matrix_normal_frozen . __dict__ [ name ] <EOL> method_frozen . __doc__ = doccer . docformat ( method . __doc__ , matnorm_docdict_noparams ) <EOL> method . __doc__ = doccer . docformat ( method . __doc__ , matnorm_docdict_params ) <EOL> _dirichlet_doc_default_callparams = """<STR_LIT>""" <EOL> _dirichlet_doc_frozen_callparams = "<STR_LIT>" <EOL> _dirichlet_doc_frozen_callparams_note = """<STR_LIT>""" <EOL> dirichlet_docdict_params = { <EOL> '<STR_LIT>' : _dirichlet_doc_default_callparams , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> dirichlet_docdict_noparams = { <EOL> '<STR_LIT>' : _dirichlet_doc_frozen_callparams , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> def _dirichlet_check_parameters ( alpha ) : <EOL> alpha = np . asarray ( alpha ) <EOL> if np . min ( alpha ) <= <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> elif alpha . ndim != <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( alpha . shape , ) ) <EOL> return alpha <EOL> def _dirichlet_check_input ( alpha , x ) : <EOL> x = np . asarray ( x ) <EOL> if x . shape [ <NUM_LIT:0> ] + <NUM_LIT:1> != alpha . shape [ <NUM_LIT:0> ] and x . shape [ <NUM_LIT:0> ] != alpha . shape [ <NUM_LIT:0> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" % ( alpha . shape , x . shape ) ) <EOL> if x . shape [ <NUM_LIT:0> ] != alpha . shape [ <NUM_LIT:0> ] : <EOL> xk = np . array ( [ <NUM_LIT:1> - np . sum ( x , <NUM_LIT:0> ) ] ) <EOL> if xk . ndim == <NUM_LIT:1> : <EOL> x = np . append ( x , xk ) <EOL> elif xk . ndim == <NUM_LIT:2> : <EOL> x = np . vstack ( ( x , xk ) ) <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> if np . min ( x ) <= <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if np . max ( x ) > <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if ( np . abs ( np . sum ( x , <NUM_LIT:0> ) - <NUM_LIT:1.0> ) > <NUM_LIT> ) . any ( ) : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % np . sum ( x , <NUM_LIT:0> ) ) <EOL> return x <EOL> def _lnB ( alpha ) : <EOL> r"""<STR_LIT>""" <EOL> return np . sum ( gammaln ( alpha ) ) - gammaln ( np . sum ( alpha ) ) <EOL> class dirichlet_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( dirichlet_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ , dirichlet_docdict_params ) <EOL> def __call__ ( self , alpha , seed = None ) : <EOL> return dirichlet_frozen ( alpha , seed = seed ) <EOL> def _logpdf ( self , x , alpha ) : <EOL> """<STR_LIT>""" <EOL> lnB = _lnB ( alpha ) <EOL> return - lnB + np . sum ( ( np . log ( x . T ) * ( alpha - <NUM_LIT:1> ) ) . T , <NUM_LIT:0> ) <EOL> def logpdf ( self , x , alpha ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> x = _dirichlet_check_input ( alpha , x ) <EOL> out = self . _logpdf ( x , alpha ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x , alpha ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> x = _dirichlet_check_input ( alpha , x ) <EOL> out = np . exp ( self . _logpdf ( x , alpha ) ) <EOL> return _squeeze_output ( out ) <EOL> def mean ( self , alpha ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> out = alpha / ( np . sum ( alpha ) ) <EOL> return _squeeze_output ( out ) <EOL> def var ( self , alpha ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> alpha0 = np . sum ( alpha ) <EOL> out = ( alpha * ( alpha0 - alpha ) ) / ( ( alpha0 * alpha0 ) * ( alpha0 + <NUM_LIT:1> ) ) <EOL> return out <EOL> def entropy ( self , alpha ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> alpha0 = np . sum ( alpha ) <EOL> lnB = _lnB ( alpha ) <EOL> K = alpha . shape [ <NUM_LIT:0> ] <EOL> out = lnB + ( alpha0 - K ) * scipy . special . psi ( alpha0 ) - np . sum ( <EOL> ( alpha - <NUM_LIT:1> ) * scipy . special . psi ( alpha ) ) <EOL> return _squeeze_output ( out ) <EOL> def rvs ( self , alpha , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> alpha = _dirichlet_check_parameters ( alpha ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> return random_state . dirichlet ( alpha , size = size ) <EOL> dirichlet = dirichlet_gen ( ) <EOL> class dirichlet_frozen ( multi_rv_frozen ) : <EOL> def __init__ ( self , alpha , seed = None ) : <EOL> self . alpha = _dirichlet_check_parameters ( alpha ) <EOL> self . _dist = dirichlet_gen ( seed ) <EOL> def logpdf ( self , x ) : <EOL> return self . _dist . logpdf ( x , self . alpha ) <EOL> def pdf ( self , x ) : <EOL> return self . _dist . pdf ( x , self . alpha ) <EOL> def mean ( self ) : <EOL> return self . _dist . mean ( self . alpha ) <EOL> def var ( self ) : <EOL> return self . _dist . var ( self . alpha ) <EOL> def entropy ( self ) : <EOL> return self . _dist . entropy ( self . alpha ) <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> return self . _dist . rvs ( self . alpha , size , random_state ) <EOL> for name in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> method = dirichlet_gen . __dict__ [ name ] <EOL> method_frozen = dirichlet_frozen . __dict__ [ name ] <EOL> method_frozen . __doc__ = doccer . docformat ( <EOL> method . __doc__ , dirichlet_docdict_noparams ) <EOL> method . __doc__ = doccer . docformat ( method . __doc__ , dirichlet_docdict_params ) <EOL> _wishart_doc_default_callparams = """<STR_LIT>""" <EOL> _wishart_doc_callparams_note = "<STR_LIT>" <EOL> _wishart_doc_frozen_callparams = "<STR_LIT>" <EOL> _wishart_doc_frozen_callparams_note = """<STR_LIT>""" <EOL> wishart_docdict_params = { <EOL> '<STR_LIT>' : _wishart_doc_default_callparams , <EOL> '<STR_LIT>' : _wishart_doc_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> wishart_docdict_noparams = { <EOL> '<STR_LIT>' : _wishart_doc_frozen_callparams , <EOL> '<STR_LIT>' : _wishart_doc_frozen_callparams_note , <EOL> '<STR_LIT>' : _doc_random_state <EOL> } <EOL> class wishart_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( wishart_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ , wishart_docdict_params ) <EOL> def __call__ ( self , df = None , scale = None , seed = None ) : <EOL> """<STR_LIT>""" <EOL> return wishart_frozen ( df , scale , seed ) <EOL> def _process_parameters ( self , df , scale ) : <EOL> if scale is None : <EOL> scale = <NUM_LIT:1.0> <EOL> scale = np . asarray ( scale , dtype = float ) <EOL> if scale . ndim == <NUM_LIT:0> : <EOL> scale = scale [ np . newaxis , np . newaxis ] <EOL> elif scale . ndim == <NUM_LIT:1> : <EOL> scale = np . diag ( scale ) <EOL> elif scale . ndim == <NUM_LIT:2> and not scale . shape [ <NUM_LIT:0> ] == scale . shape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> % str ( scale . shape ) ) <EOL> elif scale . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % scale . ndim ) <EOL> dim = scale . shape [ <NUM_LIT:0> ] <EOL> if df is None : <EOL> df = dim <EOL> elif not np . isscalar ( df ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> elif df < dim : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % df ) <EOL> return dim , df , scale <EOL> def _process_quantiles ( self , x , dim ) : <EOL> """<STR_LIT>""" <EOL> x = np . asarray ( x , dtype = float ) <EOL> if x . ndim == <NUM_LIT:0> : <EOL> x = x * np . eye ( dim ) [ : , : , np . newaxis ] <EOL> if x . ndim == <NUM_LIT:1> : <EOL> if dim == <NUM_LIT:1> : <EOL> x = x [ np . newaxis , np . newaxis , : ] <EOL> else : <EOL> x = np . diag ( x ) [ : , : , np . newaxis ] <EOL> elif x . ndim == <NUM_LIT:2> : <EOL> if not x . shape [ <NUM_LIT:0> ] == x . shape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> % str ( x . shape ) ) <EOL> x = x [ : , : , np . newaxis ] <EOL> elif x . ndim == <NUM_LIT:3> : <EOL> if not x . shape [ <NUM_LIT:0> ] == x . shape [ <NUM_LIT:1> ] : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" % str ( x . shape ) ) <EOL> elif x . ndim > <NUM_LIT:3> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" % x . ndim ) <EOL> if not x . shape [ <NUM_LIT:0> : <NUM_LIT:2> ] == ( dim , dim ) : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' % ( ( dim , dim ) , x . shape [ <NUM_LIT:0> : <NUM_LIT:2> ] ) ) <EOL> return x <EOL> def _process_size ( self , size ) : <EOL> size = np . asarray ( size ) <EOL> if size . ndim == <NUM_LIT:0> : <EOL> size = size [ np . newaxis ] <EOL> elif size . ndim > <NUM_LIT:1> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' <EOL> '<STR_LIT>' % str ( tuple ( size ) ) ) <EOL> n = size . prod ( ) <EOL> shape = tuple ( size ) <EOL> return n , shape <EOL> def _logpdf ( self , x , dim , df , scale , log_det_scale , C ) : <EOL> """<STR_LIT>""" <EOL> log_det_x = np . zeros ( x . shape [ - <NUM_LIT:1> ] ) <EOL> scale_inv_x = np . zeros ( x . shape ) <EOL> tr_scale_inv_x = np . zeros ( x . shape [ - <NUM_LIT:1> ] ) <EOL> for i in range ( x . shape [ - <NUM_LIT:1> ] ) : <EOL> _ , log_det_x [ i ] = self . _cholesky_logdet ( x [ : , : , i ] ) <EOL> scale_inv_x [ : , : , i ] = scipy . linalg . cho_solve ( ( C , True ) , x [ : , : , i ] ) <EOL> tr_scale_inv_x [ i ] = scale_inv_x [ : , : , i ] . trace ( ) <EOL> out = ( ( <NUM_LIT:0.5> * ( df - dim - <NUM_LIT:1> ) * log_det_x - <NUM_LIT:0.5> * tr_scale_inv_x ) - <EOL> ( <NUM_LIT:0.5> * df * dim * _LOG_2 + <NUM_LIT:0.5> * df * log_det_scale + <EOL> multigammaln ( <NUM_LIT:0.5> * df , dim ) ) ) <EOL> return out <EOL> def logpdf ( self , x , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> x = self . _process_quantiles ( x , dim ) <EOL> C , log_det_scale = self . _cholesky_logdet ( scale ) <EOL> out = self . _logpdf ( x , dim , df , scale , log_det_scale , C ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x , df , scale ) : <EOL> """<STR_LIT>""" <EOL> return np . exp ( self . logpdf ( x , df , scale ) ) <EOL> def _mean ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> return df * scale <EOL> def mean ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _mean ( dim , df , scale ) <EOL> return _squeeze_output ( out ) <EOL> def _mode ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> if df >= dim + <NUM_LIT:1> : <EOL> out = ( df - dim - <NUM_LIT:1> ) * scale <EOL> else : <EOL> out = None <EOL> return out <EOL> def mode ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _mode ( dim , df , scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def _var ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> var = scale ** <NUM_LIT:2> <EOL> diag = scale . diagonal ( ) <EOL> var += np . outer ( diag , diag ) <EOL> var *= df <EOL> return var <EOL> def var ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _var ( dim , df , scale ) <EOL> return _squeeze_output ( out ) <EOL> def _standard_rvs ( self , n , shape , dim , df , random_state ) : <EOL> """<STR_LIT>""" <EOL> n_tril = dim * ( dim - <NUM_LIT:1> ) // <NUM_LIT:2> <EOL> covariances = random_state . normal ( <EOL> size = n * n_tril ) . reshape ( shape + ( n_tril , ) ) <EOL> variances = np . r_ [ [ random_state . chisquare ( df - ( i + <NUM_LIT:1> ) + <NUM_LIT:1> , size = n ) ** <NUM_LIT:0.5> <EOL> for i in range ( dim ) ] ] . reshape ( ( dim , ) + shape [ : : - <NUM_LIT:1> ] ) . T <EOL> A = np . zeros ( shape + ( dim , dim ) ) <EOL> size_idx = tuple ( [ slice ( None , None , None ) ] * len ( shape ) ) <EOL> tril_idx = np . tril_indices ( dim , k = - <NUM_LIT:1> ) <EOL> A [ size_idx + tril_idx ] = covariances <EOL> diag_idx = np . diag_indices ( dim ) <EOL> A [ size_idx + diag_idx ] = variances <EOL> return A <EOL> def _rvs ( self , n , shape , dim , df , C , random_state ) : <EOL> """<STR_LIT>""" <EOL> random_state = self . _get_random_state ( random_state ) <EOL> A = self . _standard_rvs ( n , shape , dim , df , random_state ) <EOL> for index in np . ndindex ( shape ) : <EOL> CA = np . dot ( C , A [ index ] ) <EOL> A [ index ] = np . dot ( CA , CA . T ) <EOL> return A <EOL> def rvs ( self , df , scale , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> n , shape = self . _process_size ( size ) <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> C = scipy . linalg . cholesky ( scale , lower = True ) <EOL> out = self . _rvs ( n , shape , dim , df , C , random_state ) <EOL> return _squeeze_output ( out ) <EOL> def _entropy ( self , dim , df , log_det_scale ) : <EOL> """<STR_LIT>""" <EOL> return ( <EOL> <NUM_LIT:0.5> * ( dim + <NUM_LIT:1> ) * log_det_scale + <EOL> <NUM_LIT:0.5> * dim * ( dim + <NUM_LIT:1> ) * _LOG_2 + <EOL> multigammaln ( <NUM_LIT:0.5> * df , dim ) - <EOL> <NUM_LIT:0.5> * ( df - dim - <NUM_LIT:1> ) * np . sum ( <EOL> [ psi ( <NUM_LIT:0.5> * ( df + <NUM_LIT:1> - ( i + <NUM_LIT:1> ) ) ) for i in range ( dim ) ] <EOL> ) + <EOL> <NUM_LIT:0.5> * df * dim <EOL> ) <EOL> def entropy ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> _ , log_det_scale = self . _cholesky_logdet ( scale ) <EOL> return self . _entropy ( dim , df , log_det_scale ) <EOL> def _cholesky_logdet ( self , scale ) : <EOL> """<STR_LIT>""" <EOL> c_decomp = scipy . linalg . cholesky ( scale , lower = True ) <EOL> logdet = <NUM_LIT:2> * np . sum ( np . log ( c_decomp . diagonal ( ) ) ) <EOL> return c_decomp , logdet <EOL> wishart = wishart_gen ( ) <EOL> class wishart_frozen ( multi_rv_frozen ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , df , scale , seed = None ) : <EOL> self . _dist = wishart_gen ( seed ) <EOL> self . dim , self . df , self . scale = self . _dist . _process_parameters ( <EOL> df , scale ) <EOL> self . C , self . log_det_scale = self . _dist . _cholesky_logdet ( self . scale ) <EOL> def logpdf ( self , x ) : <EOL> x = self . _dist . _process_quantiles ( x , self . dim ) <EOL> out = self . _dist . _logpdf ( x , self . dim , self . df , self . scale , <EOL> self . log_det_scale , self . C ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x ) : <EOL> return np . exp ( self . logpdf ( x ) ) <EOL> def mean ( self ) : <EOL> out = self . _dist . _mean ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) <EOL> def mode ( self ) : <EOL> out = self . _dist . _mode ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def var ( self ) : <EOL> out = self . _dist . _var ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> n , shape = self . _dist . _process_size ( size ) <EOL> out = self . _dist . _rvs ( n , shape , self . dim , self . df , <EOL> self . C , random_state ) <EOL> return _squeeze_output ( out ) <EOL> def entropy ( self ) : <EOL> return self . _dist . _entropy ( self . dim , self . df , self . log_det_scale ) <EOL> for name in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> method = wishart_gen . __dict__ [ name ] <EOL> method_frozen = wishart_frozen . __dict__ [ name ] <EOL> method_frozen . __doc__ = doccer . docformat ( <EOL> method . __doc__ , wishart_docdict_noparams ) <EOL> method . __doc__ = doccer . docformat ( method . __doc__ , wishart_docdict_params ) <EOL> from numpy import asarray_chkfinite , asarray <EOL> from scipy . linalg . misc import LinAlgError <EOL> from scipy . linalg . lapack import get_lapack_funcs <EOL> def _cho_inv_batch ( a , check_finite = True ) : <EOL> """<STR_LIT>""" <EOL> if check_finite : <EOL> a1 = asarray_chkfinite ( a ) <EOL> else : <EOL> a1 = asarray ( a ) <EOL> if len ( a1 . shape ) < <NUM_LIT:2> or a1 . shape [ - <NUM_LIT:2> ] != a1 . shape [ - <NUM_LIT:1> ] : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> potrf , potri = get_lapack_funcs ( ( '<STR_LIT>' , '<STR_LIT>' ) , ( a1 , ) ) <EOL> tril_idx = np . tril_indices ( a . shape [ - <NUM_LIT:2> ] , k = - <NUM_LIT:1> ) <EOL> triu_idx = np . triu_indices ( a . shape [ - <NUM_LIT:2> ] , k = <NUM_LIT:1> ) <EOL> for index in np . ndindex ( a1 . shape [ : - <NUM_LIT:2> ] ) : <EOL> a1 [ index ] , info = potrf ( a1 [ index ] , lower = True , overwrite_a = False , <EOL> clean = False ) <EOL> if info > <NUM_LIT:0> : <EOL> raise LinAlgError ( "<STR_LIT>" <EOL> % info ) <EOL> if info < <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' % - info ) <EOL> a1 [ index ] , info = potri ( a1 [ index ] , lower = True , overwrite_c = False ) <EOL> if info > <NUM_LIT:0> : <EOL> raise LinAlgError ( "<STR_LIT>" ) <EOL> if info < <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' % - info ) <EOL> a1 [ index ] [ triu_idx ] = a1 [ index ] [ tril_idx ] <EOL> return a1 <EOL> class invwishart_gen ( wishart_gen ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( invwishart_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ , wishart_docdict_params ) <EOL> def __call__ ( self , df = None , scale = None , seed = None ) : <EOL> """<STR_LIT>""" <EOL> return invwishart_frozen ( df , scale , seed ) <EOL> def _logpdf ( self , x , dim , df , scale , log_det_scale ) : <EOL> """<STR_LIT>""" <EOL> log_det_x = np . zeros ( x . shape [ - <NUM_LIT:1> ] ) <EOL> x_inv = np . copy ( x ) . T <EOL> if dim > <NUM_LIT:1> : <EOL> _cho_inv_batch ( x_inv ) <EOL> else : <EOL> x_inv = <NUM_LIT:1.> / x_inv <EOL> tr_scale_x_inv = np . zeros ( x . shape [ - <NUM_LIT:1> ] ) <EOL> for i in range ( x . shape [ - <NUM_LIT:1> ] ) : <EOL> C , lower = scipy . linalg . cho_factor ( x [ : , : , i ] , lower = True ) <EOL> log_det_x [ i ] = <NUM_LIT:2> * np . sum ( np . log ( C . diagonal ( ) ) ) <EOL> tr_scale_x_inv [ i ] = np . dot ( scale , x_inv [ i ] ) . trace ( ) <EOL> out = ( ( <NUM_LIT:0.5> * df * log_det_scale - <NUM_LIT:0.5> * tr_scale_x_inv ) - <EOL> ( <NUM_LIT:0.5> * df * dim * _LOG_2 + <NUM_LIT:0.5> * ( df + dim + <NUM_LIT:1> ) * log_det_x ) - <EOL> multigammaln ( <NUM_LIT:0.5> * df , dim ) ) <EOL> return out <EOL> def logpdf ( self , x , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> x = self . _process_quantiles ( x , dim ) <EOL> _ , log_det_scale = self . _cholesky_logdet ( scale ) <EOL> out = self . _logpdf ( x , dim , df , scale , log_det_scale ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x , df , scale ) : <EOL> """<STR_LIT>""" <EOL> return np . exp ( self . logpdf ( x , df , scale ) ) <EOL> def _mean ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> if df > dim + <NUM_LIT:1> : <EOL> out = scale / ( df - dim - <NUM_LIT:1> ) <EOL> else : <EOL> out = None <EOL> return out <EOL> def mean ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _mean ( dim , df , scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def _mode ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> return scale / ( df + dim + <NUM_LIT:1> ) <EOL> def mode ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _mode ( dim , df , scale ) <EOL> return _squeeze_output ( out ) <EOL> def _var ( self , dim , df , scale ) : <EOL> """<STR_LIT>""" <EOL> if df > dim + <NUM_LIT:3> : <EOL> var = ( df - dim + <NUM_LIT:1> ) * scale ** <NUM_LIT:2> <EOL> diag = scale . diagonal ( ) <EOL> var += ( df - dim - <NUM_LIT:1> ) * np . outer ( diag , diag ) <EOL> var /= ( df - dim ) * ( df - dim - <NUM_LIT:1> ) ** <NUM_LIT:2> * ( df - dim - <NUM_LIT:3> ) <EOL> else : <EOL> var = None <EOL> return var <EOL> def var ( self , df , scale ) : <EOL> """<STR_LIT>""" <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> out = self . _var ( dim , df , scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def _rvs ( self , n , shape , dim , df , C , random_state ) : <EOL> """<STR_LIT>""" <EOL> random_state = self . _get_random_state ( random_state ) <EOL> A = super ( invwishart_gen , self ) . _standard_rvs ( n , shape , dim , <EOL> df , random_state ) <EOL> eye = np . eye ( dim ) <EOL> trtrs = get_lapack_funcs ( ( '<STR_LIT>' ) , ( A , ) ) <EOL> for index in np . ndindex ( A . shape [ : - <NUM_LIT:2> ] ) : <EOL> CA = np . dot ( C , A [ index ] ) <EOL> if dim > <NUM_LIT:1> : <EOL> CA , info = trtrs ( CA , eye , lower = True ) <EOL> if info > <NUM_LIT:0> : <EOL> raise LinAlgError ( "<STR_LIT>" ) <EOL> if info < <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' % - info ) <EOL> else : <EOL> CA = <NUM_LIT:1.> / CA <EOL> A [ index ] = np . dot ( CA . T , CA ) <EOL> return A <EOL> def rvs ( self , df , scale , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> n , shape = self . _process_size ( size ) <EOL> dim , df , scale = self . _process_parameters ( df , scale ) <EOL> eye = np . eye ( dim ) <EOL> L , lower = scipy . linalg . cho_factor ( scale , lower = True ) <EOL> inv_scale = scipy . linalg . cho_solve ( ( L , lower ) , eye ) <EOL> C = scipy . linalg . cholesky ( inv_scale , lower = True ) <EOL> out = self . _rvs ( n , shape , dim , df , C , random_state ) <EOL> return _squeeze_output ( out ) <EOL> def entropy ( self ) : <EOL> raise AttributeError <EOL> invwishart = invwishart_gen ( ) <EOL> class invwishart_frozen ( multi_rv_frozen ) : <EOL> def __init__ ( self , df , scale , seed = None ) : <EOL> """<STR_LIT>""" <EOL> self . _dist = invwishart_gen ( seed ) <EOL> self . dim , self . df , self . scale = self . _dist . _process_parameters ( <EOL> df , scale <EOL> ) <EOL> C , lower = scipy . linalg . cho_factor ( self . scale , lower = True ) <EOL> self . log_det_scale = <NUM_LIT:2> * np . sum ( np . log ( C . diagonal ( ) ) ) <EOL> eye = np . eye ( self . dim ) <EOL> self . inv_scale = scipy . linalg . cho_solve ( ( C , lower ) , eye ) <EOL> self . C = scipy . linalg . cholesky ( self . inv_scale , lower = True ) <EOL> def logpdf ( self , x ) : <EOL> x = self . _dist . _process_quantiles ( x , self . dim ) <EOL> out = self . _dist . _logpdf ( x , self . dim , self . df , self . scale , <EOL> self . log_det_scale ) <EOL> return _squeeze_output ( out ) <EOL> def pdf ( self , x ) : <EOL> return np . exp ( self . logpdf ( x ) ) <EOL> def mean ( self ) : <EOL> out = self . _dist . _mean ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def mode ( self ) : <EOL> out = self . _dist . _mode ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) <EOL> def var ( self ) : <EOL> out = self . _dist . _var ( self . dim , self . df , self . scale ) <EOL> return _squeeze_output ( out ) if out is not None else out <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> n , shape = self . _dist . _process_size ( size ) <EOL> out = self . _dist . _rvs ( n , shape , self . dim , self . df , <EOL> self . C , random_state ) <EOL> return _squeeze_output ( out ) <EOL> def entropy ( self ) : <EOL> raise AttributeError <EOL> for name in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] : <EOL> method = invwishart_gen . __dict__ [ name ] <EOL> method_frozen = wishart_frozen . __dict__ [ name ] <EOL> method_frozen . __doc__ = doccer . docformat ( <EOL> method . __doc__ , wishart_docdict_noparams ) <EOL> method . __doc__ = doccer . docformat ( method . __doc__ , wishart_docdict_params ) <EOL> class special_ortho_group_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( special_ortho_group_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ ) <EOL> def __call__ ( self , dim = None , seed = None ) : <EOL> """<STR_LIT>""" <EOL> return special_ortho_group_frozen ( dim , seed = seed ) <EOL> def _process_parameters ( self , dim ) : <EOL> """<STR_LIT>""" <EOL> if dim is None or not np . isscalar ( dim ) or dim <= <NUM_LIT:1> or dim != int ( dim ) : <EOL> raise ValueError ( """<STR_LIT>""" ) <EOL> return dim <EOL> def rvs ( self , dim , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> size = int ( size ) <EOL> if size > <NUM_LIT:1> : <EOL> return np . array ( [ self . rvs ( dim , size = <NUM_LIT:1> , random_state = random_state ) <EOL> for i in range ( size ) ] ) <EOL> dim = self . _process_parameters ( dim ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> H = np . eye ( dim ) <EOL> D = np . ones ( ( dim , ) ) <EOL> for n in range ( <NUM_LIT:1> , dim ) : <EOL> x = random_state . normal ( size = ( dim - n + <NUM_LIT:1> , ) ) <EOL> D [ n - <NUM_LIT:1> ] = np . sign ( x [ <NUM_LIT:0> ] ) <EOL> x [ <NUM_LIT:0> ] -= D [ n - <NUM_LIT:1> ] * np . sqrt ( ( x * x ) . sum ( ) ) <EOL> Hx = ( np . eye ( dim - n + <NUM_LIT:1> ) <EOL> - <NUM_LIT> * np . outer ( x , x ) / ( x * x ) . sum ( ) ) <EOL> mat = np . eye ( dim ) <EOL> mat [ n - <NUM_LIT:1> : , n - <NUM_LIT:1> : ] = Hx <EOL> H = np . dot ( H , mat ) <EOL> D [ - <NUM_LIT:1> ] = ( - <NUM_LIT:1> ) ** ( <NUM_LIT:1> - ( dim % <NUM_LIT:2> ) ) * D . prod ( ) <EOL> H = ( D * H . T ) . T <EOL> return H <EOL> special_ortho_group = special_ortho_group_gen ( ) <EOL> class special_ortho_group_frozen ( multi_rv_frozen ) : <EOL> def __init__ ( self , dim = None , seed = None ) : <EOL> """<STR_LIT>""" <EOL> self . _dist = special_ortho_group_gen ( seed ) <EOL> self . dim = self . _dist . _process_parameters ( dim ) <EOL> def rvs ( self , size = <NUM_LIT:1> , random_state = None ) : <EOL> return self . _dist . rvs ( self . dim , size , random_state ) <EOL> class ortho_group_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( ortho_group_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ ) <EOL> def _process_parameters ( self , dim ) : <EOL> """<STR_LIT>""" <EOL> if dim is None or not np . isscalar ( dim ) or dim <= <NUM_LIT:1> or dim != int ( dim ) : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> return dim <EOL> def rvs ( self , dim , size = <NUM_LIT:1> , random_state = None ) : <EOL> """<STR_LIT>""" <EOL> size = int ( size ) <EOL> if size > <NUM_LIT:1> : <EOL> return np . array ( [ self . rvs ( dim , size = <NUM_LIT:1> , random_state = random_state ) <EOL> for i in range ( size ) ] ) <EOL> dim = self . _process_parameters ( dim ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> H = np . eye ( dim ) <EOL> for n in range ( <NUM_LIT:1> , dim ) : <EOL> x = random_state . normal ( size = ( dim - n + <NUM_LIT:1> , ) ) <EOL> D = np . sign ( x [ <NUM_LIT:0> ] ) <EOL> x [ <NUM_LIT:0> ] += D * np . sqrt ( ( x * x ) . sum ( ) ) <EOL> Hx = - D * ( np . eye ( dim - n + <NUM_LIT:1> ) <EOL> - <NUM_LIT> * np . outer ( x , x ) / ( x * x ) . sum ( ) ) <EOL> mat = np . eye ( dim ) <EOL> mat [ n - <NUM_LIT:1> : , n - <NUM_LIT:1> : ] = Hx <EOL> H = np . dot ( H , mat ) <EOL> return H <EOL> ortho_group = ortho_group_gen ( ) <EOL> class random_correlation_gen ( multi_rv_generic ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , seed = None ) : <EOL> super ( random_correlation_gen , self ) . __init__ ( seed ) <EOL> self . __doc__ = doccer . docformat ( self . __doc__ ) <EOL> def _process_parameters ( self , eigs , tol ) : <EOL> eigs = np . asarray ( eigs , dtype = float ) <EOL> dim = eigs . size <EOL> if eigs . ndim != <NUM_LIT:1> or eigs . shape [ <NUM_LIT:0> ] != dim or dim <= <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if np . fabs ( np . sum ( eigs ) - dim ) > tol : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> for x in eigs : <EOL> if x < - tol : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> return dim , eigs <EOL> def _givens_to_1 ( self , aii , ajj , aij ) : <EOL> """<STR_LIT>""" <EOL> aiid = aii - <NUM_LIT:1.> <EOL> ajjd = ajj - <NUM_LIT:1.> <EOL> if ajjd == <NUM_LIT:0> : <EOL> return <NUM_LIT:0.> , <NUM_LIT:1.> <EOL> dd = math . sqrt ( max ( aij ** <NUM_LIT:2> - aiid * ajjd , <NUM_LIT:0> ) ) <EOL> t = ( aij + math . copysign ( dd , aij ) ) / ajjd <EOL> c = <NUM_LIT:1.> / math . sqrt ( <NUM_LIT:1.> + t * t ) <EOL> if c == <NUM_LIT:0> : <EOL> s = <NUM_LIT:1.0> <EOL> else : <EOL> s = c * t <EOL> return c , s <EOL> def _to_corr ( self , m ) : <EOL> """<STR_LIT>""" <EOL> if not ( m . flags . c_contiguous and m . dtype == np . float64 and m . shape [ <NUM_LIT:0> ] == m . shape [ <NUM_LIT:1> ] ) : <EOL> raise ValueError ( ) <EOL> d = m . shape [ <NUM_LIT:0> ] <EOL> for i in range ( d - <NUM_LIT:1> ) : <EOL> if m [ i , i ] == <NUM_LIT:1> : <EOL> continue <EOL> elif m [ i , i ] > <NUM_LIT:1> : <EOL> for j in range ( i + <NUM_LIT:1> , d ) : <EOL> if m [ j , j ] < <NUM_LIT:1> : <EOL> break <EOL> else : <EOL> for j in range ( i + <NUM_LIT:1> , d ) : <EOL> if m [ j , j ] > <NUM_LIT:1> : <EOL> break <EOL> c , s = self . _givens_to_1 ( m [ i , i ] , m [ j , j ] , m [ i , j ] ) <EOL> mv = m . ravel ( ) <EOL> drot ( mv , mv , c , - s , n = d , <EOL> offx = i * d , incx = <NUM_LIT:1> , offy = j * d , incy = <NUM_LIT:1> , <EOL> overwrite_x = True , overwrite_y = True ) <EOL> drot ( mv , mv , c , - s , n = d , <EOL> offx = i , incx = d , offy = j , incy = d , <EOL> overwrite_x = True , overwrite_y = True ) <EOL> return m <EOL> def rvs ( self , eigs , random_state = None , tol = <NUM_LIT> , diag_tol = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> dim , eigs = self . _process_parameters ( eigs , tol = tol ) <EOL> random_state = self . _get_random_state ( random_state ) <EOL> m = ortho_group . rvs ( dim , random_state = random_state ) <EOL> m = np . dot ( np . dot ( m , np . diag ( eigs ) ) , m . T ) <EOL> m = self . _to_corr ( m ) <EOL> if abs ( m . diagonal ( ) - <NUM_LIT:1> ) . max ( ) > diag_tol : <EOL> raise RuntimeError ( "<STR_LIT>" ) <EOL> return m <EOL> random_correlation = random_correlation_gen ( ) </s>
