[docs]classSkewNorm_gen(distributions.rv_continuous):'''univariate Skew-Normal distribution of Azzalini class follows scipy.stats.distributions pattern but with __init__ '''def__init__(self):#super(SkewNorm_gen,self).__init__(distributions.rv_continuous.__init__(self,name='Skew Normal distribution',shapes='alpha',extradoc=''' ''')def_argcheck(self,alpha):return1#(alpha >= 0)def_rvs(self,alpha):# see http://azzalini.stat.unipd.it/SN/faq.htmldelta=alpha/np.sqrt(1+alpha**2)u0=stats.norm.rvs(size=self._size)u1=delta*u0+np.sqrt(1-delta**2)*stats.norm.rvs(size=self._size)returnnp.where(u0>0,u1,-u1)def_munp(self,n,alpha):# use pdf integration with _mom0_sc if only _pdf is defined.# default stats calculation uses ppf, which is much slowerreturnself._mom0_sc(n,alpha)def_pdf(self,x,alpha):# 2*normpdf(x)*normcdf(alpha*x)return2.0/np.sqrt(2*np.pi)*np.exp(-x**2/2.0)*special.ndtr(alpha*x)def_stats_skip(self,x,alpha,moments='mvsk'):#skip for now to force moment integration as checkpass

skewnorm=SkewNorm_gen()# generated the same way as distributions in stats.distributions

[docs]defpdf_moments_st(cnt):"""Return the Gaussian expanded pdf function given the list of central moments (first one is mean). version of scipy.stats, any changes ? the scipy.stats version has a bug and returns normal distribution """N=len(cnt)ifN<2:raiseValueError("At least two moments must be given to ""approximate the pdf.")totp=poly1d(1)sig=sqrt(cnt[1])mu=cnt[0]ifN>2:Dvals=_hermnorm(N+1)forkinrange(3,N+1):# Find CkCk=0.0forninrange((k-3)/2):m=k-2*nifm%2:# m is oddmomdiff=cnt[m-1]else:momdiff=cnt[m-1]-sig*sig*scipy.factorial2(m-1)Ck+=Dvals[k][m]/sig**m*momdiff# Add to totpraiseSystemErrorprint(Dvals)print(Ck)totp=totp+Ck*Dvals[k]defthisfunc(x):xn=(x-mu)/sigreturntotp(xn)*exp(-xn*xn/2.0)/sqrt(2*np.pi)/sigreturnthisfunc,totp

[docs]defpdf_mvsk(mvsk):"""Return the Gaussian expanded pdf function given the list of 1st, 2nd moment and skew and Fisher (excess) kurtosis. Parameters ---------- mvsk : list of mu, mc2, skew, kurt distribution is matched to these four moments Returns ------- pdffunc : function function that evaluates the pdf(x), where x is the non-standardized random variable. Notes ----- Changed so it works only if four arguments are given. Uses explicit formula, not loop. This implements a Gram-Charlier expansion of the normal distribution where the first 2 moments coincide with those of the normal distribution but skew and kurtosis can deviate from it. In the Gram-Charlier distribution it is possible that the density becomes negative. This is the case when the deviation from the normal distribution is too large. References ---------- https://en.wikipedia.org/wiki/Edgeworth_series Johnson N.L., S. Kotz, N. Balakrishnan: Continuous Univariate Distributions, Volume 1, 2nd ed., p.30 """N=len(mvsk)ifN<4:raiseValueError("Four moments must be given to ""approximate the pdf.")mu,mc2,skew,kurt=mvsktotp=poly1d(1)sig=sqrt(mc2)ifN>2:Dvals=_hermnorm(N+1)C3=skew/6.0C4=kurt/24.0# Note: Hermite polynomial for order 3 in _hermnorm is negative# instead of positivetotp=totp-C3*Dvals[3]+C4*Dvals[4]defpdffunc(x):xn=(x-mu)/sigreturntotp(xn)*np.exp(-xn*xn/2.0)/np.sqrt(2*np.pi)/sigreturnpdffunc

[docs]defpdf_moments(cnt):"""Return the Gaussian expanded pdf function given the list of central moments (first one is mean). Changed so it works only if four arguments are given. Uses explicit formula, not loop. Notes ----- This implements a Gram-Charlier expansion of the normal distribution where the first 2 moments coincide with those of the normal distribution but skew and kurtosis can deviate from it. In the Gram-Charlier distribution it is possible that the density becomes negative. This is the case when the deviation from the normal distribution is too large. References ---------- https://en.wikipedia.org/wiki/Edgeworth_series Johnson N.L., S. Kotz, N. Balakrishnan: Continuous Univariate Distributions, Volume 1, 2nd ed., p.30 """N=len(cnt)ifN<2:raiseValueError("At least two moments must be given to ""approximate the pdf.")mc,mc2,mc3,mc4=cntskew=mc3/mc2**1.5kurt=mc4/mc2**2.0-3.0# Fisher kurtosis, excess kurtosistotp=poly1d(1)sig=sqrt(cnt[1])mu=cnt[0]ifN>2:Dvals=_hermnorm(N+1)## for k in range(3,N+1):## # Find Ck## Ck = 0.0## for n in range((k-3)/2):## m = k-2*n## if m % 2: # m is odd## momdiff = cnt[m-1]## else:## momdiff = cnt[m-1] - sig*sig*scipy.factorial2(m-1)## Ck += Dvals[k][m] / sig**m * momdiff## # Add to totp## raise## print Dvals## print Ck## totp = totp + Ck*Dvals[k]C3=skew/6.0C4=kurt/24.0totp=totp-C3*Dvals[3]+C4*Dvals[4]defthisfunc(x):xn=(x-mu)/sigreturntotp(xn)*np.exp(-xn*xn/2.0)/np.sqrt(2*np.pi)/sigreturnthisfunc

[docs]classNormExpan_gen(distributions.rv_continuous):'''Gram-Charlier Expansion of Normal distribution class follows scipy.stats.distributions pattern but with __init__ '''def__init__(self,args,**kwds):#todo: replace with super calldistributions.rv_continuous.__init__(self,name='Normal Expansion distribution',shapes=' ',extradoc=''' The distribution is defined as the Gram-Charlier expansion of the normal distribution using the first four moments. The pdf is given by pdf(x) = (1+ skew/6.0 * H(xc,3) + kurt/24.0 * H(xc,4))*normpdf(xc) where xc = (x-mu)/sig is the standardized value of the random variable and H(xc,3) and H(xc,4) are Hermite polynomials Note: This distribution has to be parametrized during initialization and instantiation, and does not have a shape parameter after instantiation (similar to frozen distribution except for location and scale.) Location and scale can be used as with other distributions, however note, that they are relative to the initialized distribution. ''')#print args, kwdsmode=kwds.get('mode','sample')ifmode=='sample':mu,sig,sk,kur=stats.describe(args)[2:]self.mvsk=(mu,sig,sk,kur)cnt=mvsk2mc((mu,sig,sk,kur))elifmode=='mvsk':cnt=mvsk2mc(args)self.mvsk=argselifmode=='centmom':cnt=argsself.mvsk=mc2mvsk(cnt)else:raiseValueError("mode must be 'mvsk' or centmom")self.cnt=cnt#self.mvsk = (mu,sig,sk,kur)#self._pdf = pdf_moments(cnt)self._pdf=pdf_mvsk(self.mvsk)def_munp(self,n):# use pdf integration with _mom0_sc if only _pdf is defined.# default stats calculation uses ppfreturnself._mom0_sc(n)def_stats_skip(self):# skip for now to force numerical integration of pdf for testingreturnself.mvsk

## copied from nonlinear_transform_gen.py''' A class for the distribution of a non-linear monotonic transformation of a continuous random variablesimplest usage:example: create log-gamma distribution, i.e. y = log(x), where x is gamma distributed (also available in scipy.stats) loggammaexpg = Transf_gen(stats.gamma, np.log, np.exp)example: what is the distribution of the discount factor y=1/(1+x) where interest rate x is normally distributed with N(mux,stdx**2)')? (just to come up with a story that implies a nice transformation) invnormalg = Transf_gen(stats.norm, inversew, inversew_inv, decr=True, a=-np.inf)This class does not work well for distributions with difficult shapes, e.g. 1/x where x is standard normal, because of the singularity and jump at zero.Note: I'm working from my version of scipy.stats.distribution. But this script runs under scipy 0.6.0 (checked with numpy: 1.2.0rc2 and python 2.4)This is not yet thoroughly tested, polished or optimizedTODO: * numargs handling is not yet working properly, numargs needs to be specified (default = 0 or 1) * feeding args and kwargs to underlying distribution is untested and incomplete * distinguish args and kwargs for the transformed and the underlying distribution - currently all args and no kwargs are transmitted to underlying distribution - loc and scale only work for transformed, but not for underlying distribution - possible to separate args for transformation and underlying distribution parameters * add _rvs as method, will be faster in many casesCreated on Tuesday, October 28, 2008, 12:40:37 PMAuthor: josef-pktdLicense: BSD'''defget_u_argskwargs(**kwargs):#Todo: What's this? wrong spacing, used in Transf_gen TransfTwo_genu_kwargs=dict((k.replace('u_','',1),v)fork,viniteritems(kwargs)ifk.startswith('u_'))u_args=u_kwargs.pop('u_args',None)returnu_args,u_kwargsclassTransf_gen(distributions.rv_continuous):'''a class for non-linear monotonic transformation of a continuous random variable '''def__init__(self,kls,func,funcinv,*args,**kwargs):#print args#print kwargsself.func=funcself.funcinv=funcinv#explicit for self.__dict__.update(kwargs)#need to set numargs because inspection does not workself.numargs=kwargs.pop('numargs',0)#print self.numargsname=kwargs.pop('name','transfdist')longname=kwargs.pop('longname','Non-linear transformed distribution')extradoc=kwargs.pop('extradoc',None)a=kwargs.pop('a',-np.inf)b=kwargs.pop('b',np.inf)self.decr=kwargs.pop('decr',False)#defines whether it is a decreasing (True)# or increasing (False) monotonic transformationself.u_args,self.u_kwargs=get_u_argskwargs(**kwargs)self.kls=kls#(self.u_args, self.u_kwargs)# possible to freeze the underlying distributionsuper(Transf_gen,self).__init__(a=a,b=b,name=name,longname=longname,extradoc=extradoc)def_rvs(self,*args,**kwargs):self.kls._size=self._sizereturnself.funcinv(self.kls._rvs(*args))def_cdf(self,x,*args,**kwargs):#print argsifnotself.decr:returnself.kls._cdf(self.funcinv(x),*args,**kwargs)#note scipy _cdf only take *args not *kwargselse:return1.0-self.kls._cdf(self.funcinv(x),*args,**kwargs)def_ppf(self,q,*args,**kwargs):ifnotself.decr:returnself.func(self.kls._ppf(q,*args,**kwargs))else:returnself.func(self.kls._ppf(1-q,*args,**kwargs))definverse(x):returnnp.divide(1.0,x)mux,stdx=0.05,0.1mux,stdx=9.0,1.0definversew(x):return1.0/(1+mux+x*stdx)definversew_inv(x):return(1.0/x-1.0-mux)/stdx#.np.divide(1.0,x)-10defidentit(x):returnxinvdnormalg=Transf_gen(stats.norm,inversew,inversew_inv,decr=True,#a=-np.inf,numargs=0,name='discf',longname='normal-based discount factor',extradoc='\ndistribution of discount factor y=1/(1+x)) with x N(0.05,0.1**2)')lognormalg=Transf_gen(stats.norm,np.exp,np.log,numargs=2,a=0,name='lnnorm',longname='Exp transformed normal',extradoc='\ndistribution of y = exp(x), with x standard normal''precision for moment andstats is not very high, 2-3 decimals')loggammaexpg=Transf_gen(stats.gamma,np.log,np.exp,numargs=1)## copied form nonlinear_transform_short.py'''univariate distribution of a non-linear monotonic transformation of arandom variable'''classExpTransf_gen(distributions.rv_continuous):'''Distribution based on log/exp transformation the constructor can be called with a distribution class and generates the distribution of the transformed random variable '''def__init__(self,kls,*args,**kwargs):#print args#print kwargs#explicit for self.__dict__.update(kwargs)if'numargs'inkwargs:self.numargs=kwargs['numargs']else:self.numargs=1if'name'inkwargs:name=kwargs['name']else:name='Log transformed distribution'if'a'inkwargs:a=kwargs['a']else:a=0super(ExpTransf_gen,self).__init__(a=0,name=name)self.kls=klsdef_cdf(self,x,*args):pass#print argsreturnself.kls.cdf(np.log(x),*args)def_ppf(self,q,*args):returnnp.exp(self.kls.ppf(q,*args))classLogTransf_gen(distributions.rv_continuous):'''Distribution based on log/exp transformation the constructor can be called with a distribution class and generates the distribution of the transformed random variable '''def__init__(self,kls,*args,**kwargs):#explicit for self.__dict__.update(kwargs)if'numargs'inkwargs:self.numargs=kwargs['numargs']else:self.numargs=1if'name'inkwargs:name=kwargs['name']else:name='Log transformed distribution'if'a'inkwargs:a=kwargs['a']else:a=0super(LogTransf_gen,self).__init__(a=a,name=name)self.kls=klsdef_cdf(self,x,*args):#print argsreturnself.kls._cdf(np.exp(x),*args)def_ppf(self,q,*args):returnnp.log(self.kls._ppf(q,*args))## copied from transformtwo.py'''Created on Apr 28, 2009@author: Josef Perktold'''''' A class for the distribution of a non-linear u-shaped or hump shaped transformation of acontinuous random variableThis is a companion to the distributions of non-linear monotonic transformation to the casewhen the inverse mapping is a 2-valued correspondence, for example for absolute value or squaresimplest usage:example: create squared distribution, i.e. y = x**2, where x is normal or t distributedThis class does not work well for distributions with difficult shapes, e.g. 1/x where x is standard normal, because of the singularity and jump at zero.This verifies for normal - chi2, normal - halfnorm, foldnorm, and t - FTODO: * numargs handling is not yet working properly, numargs needs to be specified (default = 0 or 1) * feeding args and kwargs to underlying distribution works in t distribution example * distinguish args and kwargs for the transformed and the underlying distribution - currently all args and no kwargs are transmitted to underlying distribution - loc and scale only work for transformed, but not for underlying distribution - possible to separate args for transformation and underlying distribution parameters * add _rvs as method, will be faster in many cases'''classTransfTwo_gen(distributions.rv_continuous):'''Distribution based on a non-monotonic (u- or hump-shaped transformation) the constructor can be called with a distribution class, and functions that define the non-linear transformation. and generates the distribution of the transformed random variable Note: the transformation, it's inverse and derivatives need to be fully specified: func, funcinvplus, funcinvminus, derivplus, derivminus. Currently no numerical derivatives or inverse are calculated This can be used to generate distribution instances similar to the distributions in scipy.stats. '''#a class for non-linear non-monotonic transformation of a continuous random variabledef__init__(self,kls,func,funcinvplus,funcinvminus,derivplus,derivminus,*args,**kwargs):#print args#print kwargsself.func=funcself.funcinvplus=funcinvplusself.funcinvminus=funcinvminusself.derivplus=derivplusself.derivminus=derivminus#explicit for self.__dict__.update(kwargs)#need to set numargs because inspection does not workself.numargs=kwargs.pop('numargs',0)#print self.numargsname=kwargs.pop('name','transfdist')longname=kwargs.pop('longname','Non-linear transformed distribution')extradoc=kwargs.pop('extradoc',None)a=kwargs.pop('a',-np.inf)# attached to self in superb=kwargs.pop('b',np.inf)# self.a, self.b would be overwrittenself.shape=kwargs.pop('shape',False)#defines whether it is a `u` shaped or `hump' shaped# transformationself.u_args,self.u_kwargs=get_u_argskwargs(**kwargs)self.kls=kls#(self.u_args, self.u_kwargs)# possible to freeze the underlying distributionsuper(TransfTwo_gen,self).__init__(a=a,b=b,name=name,shapes=kls.shapes,longname=longname,extradoc=extradoc)# add enough info for self.freeze() to be able to reconstruct the instancetry:self._ctor_param.update(dict(kls=kls,func=func,funcinvplus=funcinvplus,funcinvminus=funcinvminus,derivplus=derivplus,derivminus=derivminus,shape=self.shape))exceptAttributeError:# scipy < 0.14 does not have this, ignore and do nothingpassdef_rvs(self,*args):self.kls._size=self._size#size attached to self, not function argumentreturnself.func(self.kls._rvs(*args))def_pdf(self,x,*args,**kwargs):#print argsifself.shape=='u':signpdf=1elifself.shape=='hump':signpdf=-1else:raiseValueError('shape can only be `u` or `hump`')returnsignpdf*(self.derivplus(x)*self.kls._pdf(self.funcinvplus(x),*args,**kwargs)-self.derivminus(x)*self.kls._pdf(self.funcinvminus(x),*args,**kwargs))#note scipy _cdf only take *args not *kwargsdef_cdf(self,x,*args,**kwargs):#print argsifself.shape=='u':returnself.kls._cdf(self.funcinvplus(x),*args,**kwargs)- \
self.kls._cdf(self.funcinvminus(x),*args,**kwargs)#note scipy _cdf only take *args not *kwargselse:return1.0-self._sf(x,*args,**kwargs)def_sf(self,x,*args,**kwargs):#print argsifself.shape=='hump':returnself.kls._cdf(self.funcinvplus(x),*args,**kwargs)- \
self.kls._cdf(self.funcinvminus(x),*args,**kwargs)#note scipy _cdf only take *args not *kwargselse:return1.0-self._cdf(x,*args,**kwargs)def_munp(self,n,*args,**kwargs):returnself._mom0_sc(n,*args)# ppf might not be possible in general case?# should be possible in symmetric case# def _ppf(self, q, *args, **kwargs):# if self.shape == 'u':# return self.func(self.kls._ppf(q,*args, **kwargs))# elif self.shape == 'hump':# return self.func(self.kls._ppf(1-q,*args, **kwargs))#TODO: rename these functions to have unique namesclassSquareFunc(object):'''class to hold quadratic function with inverse function and derivative using instance methods instead of class methods, if we want extension to parametrized function '''definverseplus(self,x):returnnp.sqrt(x)definverseminus(self,x):return0.0-np.sqrt(x)defderivplus(self,x):return0.5/np.sqrt(x)defderivminus(self,x):return0.0-0.5/np.sqrt(x)defsquarefunc(self,x):returnnp.power(x,2)sqfunc=SquareFunc()squarenormalg=TransfTwo_gen(stats.norm,sqfunc.squarefunc,sqfunc.inverseplus,sqfunc.inverseminus,sqfunc.derivplus,sqfunc.derivminus,shape='u',a=0.0,b=np.inf,numargs=0,name='squarenorm',longname='squared normal distribution',extradoc='\ndistribution of the square of a normal random variable'+\
' y=x**2 with x N(0.0,1)')#u_loc=l, u_scale=s)squaretg=TransfTwo_gen(stats.t,sqfunc.squarefunc,sqfunc.inverseplus,sqfunc.inverseminus,sqfunc.derivplus,sqfunc.derivminus,shape='u',a=0.0,b=np.inf,numargs=1,name='squarenorm',longname='squared t distribution',extradoc='\ndistribution of the square of a t random variable'+\
' y=x**2 with x t(dof,0.0,1)')definverseplus(x):returnnp.sqrt(-x)definverseminus(x):return0.0-np.sqrt(-x)defderivplus(x):return0.0-0.5/np.sqrt(-x)defderivminus(x):return0.5/np.sqrt(-x)defnegsquarefunc(x):return-np.power(x,2)negsquarenormalg=TransfTwo_gen(stats.norm,negsquarefunc,inverseplus,inverseminus,derivplus,derivminus,shape='hump',a=-np.inf,b=0.0,numargs=0,name='negsquarenorm',longname='negative squared normal distribution',extradoc='\ndistribution of the negative square of a normal random variable'+\
' y=-x**2 with x N(0.0,1)')#u_loc=l, u_scale=s)definverseplus(x):returnxdefinverseminus(x):return0.0-xdefderivplus(x):return1.0defderivminus(x):return0.0-1.0defabsfunc(x):returnnp.abs(x)absnormalg=TransfTwo_gen(stats.norm,np.abs,inverseplus,inverseminus,derivplus,derivminus,shape='u',a=0.0,b=np.inf,numargs=0,name='absnorm',longname='absolute of normal distribution',extradoc='\ndistribution of the absolute value of a normal random variable'+\
' y=abs(x) with x N(0,1)')#copied from mvncdf.py'''multivariate normal probabilities and cumulative distribution functiona wrapper for scipy.stats.kde.mvndst SUBROUTINE MVNDST( N, LOWER, UPPER, INFIN, CORREL, MAXPTS, & ABSEPS, RELEPS, ERROR, VALUE, INFORM )** A subroutine for computing multivariate normal probabilities.* This subroutine uses an algorithm given in the paper* "Numerical Computation of Multivariate Normal Probabilities", in* J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by* Alan Genz* Department of Mathematics* Washington State University* Pullman, WA 99164-3113* Email : AlanGenz@wsu.edu** Parameters** N INTEGER, the number of variables.* LOWER REAL, array of lower integration limits.* UPPER REAL, array of upper integration limits.* INFIN INTEGER, array of integration limits flags:* if INFIN(I) < 0, Ith limits are (-infinity, infinity);* if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];* if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].* CORREL REAL, array of correlation coefficients; the correlation* coefficient in row I column J of the correlation matrix* should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.* The correlation matrix must be positive semidefinite.* MAXPTS INTEGER, maximum number of function values allowed. This* parameter can be used to limit the time. A sensible* strategy is to start with MAXPTS = 1000*N, and then* increase MAXPTS if ERROR is too large.* ABSEPS REAL absolute error tolerance.* RELEPS REAL relative error tolerance.* ERROR REAL estimated absolute error, with 99% confidence level.* VALUE REAL estimated value for the integral* INFORM INTEGER, termination status parameter:* if INFORM = 0, normal completion with ERROR < EPS;* if INFORM = 1, completion with ERROR > EPS and MAXPTS* function vaules used; increase MAXPTS to* decrease ERROR;* if INFORM = 2, N > 500 or N < 1.*>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[10.0,10.0],[0,0],[0.5])(2e-016, 1.0, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[100.0,100.0],[0,0],[0.0])(2e-016, 1.0, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[1.0,1.0],[0,0],[0.0])(2e-016, 0.70786098173714096, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.001,1.0],[0,0],[0.0])(2e-016, 0.42100802096993045, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.001,10.0],[0,0],[0.0])(2e-016, 0.50039894221391101, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.001,100.0],[0,0],[0.0])(2e-016, 0.50039894221391101, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.01,100.0],[0,0],[0.0])(2e-016, 0.5039893563146316, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.1,100.0],[0,0],[0.0])(2e-016, 0.53982783727702899, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.1,100.0],[2,2],[0.0])(2e-016, 0.019913918638514494, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.0,0.0],[0,0],[0.0])(2e-016, 0.25, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.0,0.0],[-1,0],[0.0])(2e-016, 0.5, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.0,0.0],[-1,0],[0.5])(2e-016, 0.5, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.0,0.0],[0,0],[0.5])(2e-016, 0.33333333333333337, 0)>>> scipy.stats.kde.mvn.mvndst([0.0,0.0],[0.0,0.0],[0,0],[0.99])(2e-016, 0.47747329317779391, 0)'''#from scipy.stats import kdeinformcode={0:'normal completion with ERROR < EPS',1:'''completion with ERROR > EPS and MAXPTS function values used; increase MAXPTS to decrease ERROR;''',2:'N > 500 or N < 1'}

[docs]defmvstdnormcdf(lower,upper,corrcoef,**kwds):'''standardized multivariate normal cumulative distribution function This is a wrapper for scipy.stats.kde.mvn.mvndst which calculates a rectangular integral over a standardized multivariate normal distribution. This function assumes standardized scale, that is the variance in each dimension is one, but correlation can be arbitrary, covariance = correlation matrix Parameters ---------- lower, upper : array_like, 1d lower and upper integration limits with length equal to the number of dimensions of the multivariate normal distribution. It can contain -np.inf or np.inf for open integration intervals corrcoef : float or array_like specifies correlation matrix in one of three ways, see notes optional keyword parameters to influence integration * maxpts : int, maximum number of function values allowed. This parameter can be used to limit the time. A sensible strategy is to start with `maxpts` = 1000*N, and then increase `maxpts` if ERROR is too large. * abseps : float absolute error tolerance. * releps : float relative error tolerance. Returns ------- cdfvalue : float value of the integral Notes ----- The correlation matrix corrcoef can be given in 3 different ways If the multivariate normal is two-dimensional than only the correlation coefficient needs to be provided. For general dimension the correlation matrix can be provided either as a one-dimensional array of the upper triangular correlation coefficients stacked by rows, or as full square correlation matrix See Also -------- mvnormcdf : cdf of multivariate normal distribution without standardization Examples -------- >>> print(mvstdnormcdf([-np.inf,-np.inf], [0.0,np.inf], 0.5)) 0.5 >>> corr = [[1.0, 0, 0.5],[0,1,0],[0.5,0,1]] >>> print(mvstdnormcdf([-np.inf,-np.inf,-100.0], [0.0,0.0,0.0], corr, abseps=1e-6)) 0.166666399198 >>> print(mvstdnormcdf([-np.inf,-np.inf,-100.0],[0.0,0.0,0.0],corr, abseps=1e-8)) something wrong completion with ERROR > EPS and MAXPTS function values used; increase MAXPTS to decrease ERROR; 1.048330348e-006 0.166666546218 >>> print(mvstdnormcdf([-np.inf,-np.inf,-100.0],[0.0,0.0,0.0], corr, \ maxpts=100000, abseps=1e-8)) 0.166666588293 '''n=len(lower)#do not know if converting to array is necessary,#but it makes ndim check possiblelower=np.array(lower)upper=np.array(upper)corrcoef=np.array(corrcoef)correl=np.zeros(int(n*(n-1)/2.0))#dtype necessary?if(lower.ndim!=1)or(upper.ndim!=1):raiseValueError('can handle only 1D bounds')iflen(upper)!=n:raiseValueError('bounds have different lengths')ifn==2andcorrcoef.size==1:correl=corrcoef#print 'case scalar rho', nelifcorrcoef.ndim==1andlen(corrcoef)==n*(n-1)/2.0:#print 'case flat corr', corrcoeff.shapecorrel=corrcoefelifcorrcoef.shape==(n,n):#print 'case square corr', correl.shapecorrel=corrcoef[np.tril_indices(n,-1)]# for ii in range(n):# for jj in range(ii):# correl[ jj + ((ii-2)*(ii-1))/2] = corrcoef[ii,jj]else:raiseValueError('corrcoef has incorrect dimension')if'maxpts'notinkwds:ifn>2:kwds['maxpts']=10000*nlowinf=np.isneginf(lower)uppinf=np.isposinf(upper)infin=2.0*np.ones(n)np.putmask(infin,lowinf,0)# infin.putmask(0,lowinf)np.putmask(infin,uppinf,1)#infin.putmask(1,uppinf)#this has to be lastnp.putmask(infin,lowinf*uppinf,-1)## #remove infs## np.putmask(lower,lowinf,-100)# infin.putmask(0,lowinf)## np.putmask(upper,uppinf,100) #infin.putmask(1,uppinf)#print lower,',',upper,',',infin,',',correl#print correl.shape#print kwds.items()error,cdfvalue,inform=scipy.stats.kde.mvn.mvndst(lower,upper,infin,correl,**kwds)ifinform:print('something wrong',informcode[inform],error)returncdfvalue

[docs]defmvnormcdf(upper,mu,cov,lower=None,**kwds):'''multivariate normal cumulative distribution function This is a wrapper for scipy.stats.kde.mvn.mvndst which calculates a rectangular integral over a multivariate normal distribution. Parameters ---------- lower, upper : array_like, 1d lower and upper integration limits with length equal to the number of dimensions of the multivariate normal distribution. It can contain -np.inf or np.inf for open integration intervals mu : array_lik, 1d list or array of means cov : array_like, 2d specifies covariance matrix optional keyword parameters to influence integration * maxpts : int, maximum number of function values allowed. This parameter can be used to limit the time. A sensible strategy is to start with `maxpts` = 1000*N, and then increase `maxpts` if ERROR is too large. * abseps : float absolute error tolerance. * releps : float relative error tolerance. Returns ------- cdfvalue : float value of the integral Notes ----- This function normalizes the location and scale of the multivariate normal distribution and then uses `mvstdnormcdf` to call the integration. See Also -------- mvstdnormcdf : location and scale standardized multivariate normal cdf '''upper=np.array(upper)iflowerisNone:lower=-np.ones(upper.shape)*np.infelse:lower=np.array(lower)cov=np.array(cov)stdev=np.sqrt(np.diag(cov))# standard deviation vector#do I need to make sure stdev is float and not int?#is this correct to normalize to corr?lower=(lower-mu)/stdevupper=(upper-mu)/stdevdivrow=np.atleast_2d(stdev)corr=cov/divrow/divrow.T#v/np.sqrt(np.atleast_2d(np.diag(covv)))/np.sqrt(np.atleast_2d(np.diag(covv))).Treturnmvstdnormcdf(lower,upper,corr,**kwds)