from __future__ import absolute_import from numpy import (r_, minimum, maximum, atleast_1d, atleast_2d, mod, ones, floor, random, eye, nonzero, where, repeat, sqrt, exp, inf, diag, zeros, sin, arcsin, nan) from numpy import triu from scipy.special import ndtr as cdfnorm, ndtri as invnorm from scipy.special import erfc import warnings import numpy as np try: from wafo import mvn # @UnresolvedImport except ImportError: warnings.warn('mvn not found. Check its compilation.') mvn = None try: from wafo import mvnprdmod # @UnresolvedImport except ImportError: warnings.warn('mvnprdmod not found. Check its compilation.') mvnprdmod = None try: from wafo import rindmod # @UnresolvedImport except ImportError: warnings.warn('rindmod not found. Check its compilation.') rindmod = None __all__ = ['Rind', 'rindmod', 'mvnprdmod', 'mvn', 'cdflomax', 'prbnormtndpc', 'prbnormndpc', 'prbnormnd', 'cdfnorm2d', 'prbnorm2d', 'cdfnorm', 'invnorm'] class Rind(object): """ RIND Computes multivariate normal expectations Parameters ---------- S : array-like, shape Ntdc x Ntdc Covariance matrix of X=[Xt,Xd,Xc] (Ntdc = Nt+Nd+Nc) m : array-like, size Ntdc expectation of X=[Xt,Xd,Xc] Blo, Bup : array-like, shape Mb x Nb Lower and upper barriers used to compute the integration limits, Hlo and Hup, respectively. indI : array-like, length Ni vector of indices to the different barriers in the indicator function. (NB! restriction indI(1)=-1, indI(NI)=Nt+Nd, Ni = Nb+1) (default indI = 0:Nt+Nd) xc : values to condition on (default xc = zeros(0,1)), size Nc x Nx Nt : size of Xt (default Nt = Ntdc - Nc) Returns ------- val: ndarray, size Nx expectation/density as explained below err, terr : ndarray, size Nx estimated sampling error and estimated truncation error, respectively. (err is with 99 confidence level) Notes ----- RIND computes multivariate normal expectations, i.e., E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix)) where "Indicator" = I{ Hlo(i) < X(i) < Hup(i), i = 1:N_t+N_d } "Jacobian" = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is "Jacobian" = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)| "condition" = Xc=xc(:,ix), ix=1,...,Nx. X = [Xt, Xd, Xc], a stochastic vector of Multivariate Gaussian variables where Xt,Xd and Xc have the length Nt,Nd and Nc, respectively (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10) Multivariate probability is computed if Nd = 0. If Mb>> import wafo.gaussian as wg >>> import numpy as np >>> n = 5 >>> Blo =-np.inf; Bup=-1.2; indI=np.array([-1, n-1], dtype=int) # Barriers >>> m = np.zeros(n); rho = 0.3; >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n) >>> rind = wg.Rind() >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI) # exact prob. 0.001946 >>> A = np.repeat(Blo,n); B = np.repeat(Bup,n) # Integration limits >>> E1 = rind(np.triu(Sc),m,A,B) #same as E0 Compute expectation E( abs(X1*X2*...*X5) ) >>> xc = np.zeros((0,1)) >>> infinity = 37 >>> dev = np.sqrt(np.diag(Sc)) # std >>> ind = np.nonzero(indI[1:])[0] >>> Bup, Blo = np.atleast_2d(Bup,Blo) >>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]]) >>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]]) >>> val, err, terr = rind(Sc,m,Blo,Bup,indI, xc, nt=0) >>> np.allclose(val, 0.05494076, rtol=4e-2) True >>> err[0] < 5e-3, terr[0] < 1e-7 (True, True) Compute expectation E( X1^{+}*X2^{+} ) with random correlation coefficient,Cov(X1,X2) = rho2. >>> m2 = [0, 0]; rho2 = np.random.rand(1) >>> Sc2 = [[1, rho2], [rho2 ,1]] >>> Blo2 = 0; Bup2 = np.inf; indI2 = [-1, 1] >>> rind2 = wg.Rind(method=1) >>> def g2(x): ... return (x*(np.pi/2+np.arcsin(x))+np.sqrt(1-x**2))/(2*np.pi) >>> E2 = g2(rho2) # exact value >>> E3 = rind(Sc2,m2,Blo2,Bup2,indI2,nt=0) >>> E4 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0) >>> E5 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0,abseps=1e-4) See also -------- prbnormnd, prbnormndpc References ---------- Podgorski et al. (2000) "Exact distributions for apparent waves in irregular seas" Ocean Engineering, Vol 27, no 1, pp979-1016. P. A. Brodtkorb (2004), Numerical evaluation of multinormal expectations In Lund university report series and in the Dr.Ing thesis: The probability of Occurrence of dangerous Wave Situations at Sea. Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, Trondheim, Norway. Per A. Brodtkorb (2006) "Evaluating Nearly Singular Multinormal Expectations with Application to Wave Distributions", Methodology And Computing In Applied Probability, Volume 8, Number 1, pp. 65-91(27) """ def __init__(self, **kwds): """ Parameters ---------- method : integer, optional defining the integration method 0 Integrate by Gauss-Legendre quadrature (Podgorski et al. 1999) 1 Integrate by SADAPT for Ndim<9 and by KRBVRC otherwise 2 Integrate by SADAPT for Ndim<20 and by KRBVRC otherwise 3 Integrate by KRBVRC by Genz (1993) (Fast Ndim<101) (default) 4 Integrate by KROBOV by Genz (1992) (Fast Ndim<101) 5 Integrate by RCRUDE by Genz (1992) (Slow Ndim<1001) 6 Integrate by SOBNIED (Fast Ndim<1041) 7 Integrate by DKBVRC by Genz (2003) (Fast Ndim<1001) xcscale : real scalar, optional scales the conditinal probability density, i.e., f_{Xc} = exp(-0.5*Xc*inv(Sxc)*Xc + XcScale) (default XcScale=0) abseps, releps : real scalars, optional absolute and relative error tolerance. (default abseps=0, releps=1e-3) coveps : real scalar, optional error tolerance in Cholesky factorization (default 1e-13) maxpts, minpts : scalar integers, optional maximum and minimum number of function values allowed. The parameter, maxpts, 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. (Only for METHOD~=0) (default maxpts=40000, minpts=0) seed : scalar integer, optional seed to the random generator used in the integrations (Only for METHOD~=0)(default floor(rand*1e9)) nit : scalar integer, optional maximum number of Xt variables to integrate. This parameter can be used to limit the time. If NIT is less than the rank of the covariance matrix, the returned result is a upper bound for the true value of the integral. (default 1000) xcutoff : real scalar, optional cut off value where the marginal normal distribution is truncated. (Depends on requested accuracy. A value between 4 and 5 is reasonable.) xsplit : real scalar parameter controlling performance of quadrature integration: if Hup>=xCutOff AND Hlo<-XSPLIT OR Hup>=XSPLIT AND Hlo<=-xCutOff then do a different integration to increase speed in rind2 and rindnit. This give slightly different results if XSPILT>=xCutOff => do the same integration always (Only for METHOD==0)(default XSPLIT = 1.5) quadno : scalar integer Quadrature formulae number used in integration of Xd variables. This number implicitly determines number of nodes used. (Only for METHOD==0) speed : scalar integer defines accuracy of calculations by choosing different parameters, possible values: 1,2...,9 (9 fastest, default []). If not speed is None the parameters, ABSEPS, RELEPS, COVEPS, XCUTOFF, MAXPTS and QUADNO will be set according to INITOPTIONS. nc1c2 : scalar integer, optional number of times to use the regression equation to restrict integration area. Nc1c2 = 1,2 is recommended. (default 2) (note: works only for method >0) """ self.method = 3 self.xcscale = 0 self.abseps = 0 self.releps = 1e-3, self.coveps = 1e-10 self.maxpts = 40000 self.minpts = 0 self.seed = None self.nit = 1000, self.xcutoff = None self.xsplit = 1.5 self.quadno = 2 self.speed = None self.nc1c2 = 2 self.__dict__.update(**kwds) self.initialize(self.speed) self.set_constants() def initialize(self, speed=None): """ Initializes member variables according to speed. Parameter --------- speed : scalar integer defining accuracy of calculations. Valid numbers: 1,2,...,10 (1=slowest and most accurate,10=fastest, but less accuracy) Member variables initialized according to speed: ----------------------------------------------- speed : Integer defining accuracy of calculations. abseps : Absolute error tolerance. releps : Relative error tolerance. covep : Error tolerance in Cholesky factorization. xcutoff : Truncation limit of the normal CDF maxpts : Maximum number of function values allowed. quadno : Quadrature formulae used in integration of Xd(i) implicitly determining # nodes """ if speed is None: return self.speed = min(max(speed, 1), 13) self.maxpts = 10000 self.quadno = r_[1:4] + (10 - min(speed, 9)) + (speed == 1) if speed in (11, 12, 13): self.abseps = 1e-1 elif speed == 10: self.abseps = 1e-2 elif speed in (7, 8, 9): self.abseps = 1e-2 elif speed in (4, 5, 6): self.maxpts = 20000 self.abseps = 1e-3 elif speed in (1, 2, 3): self.maxpts = 30000 self.abseps = 1e-4 if speed < 12: tmp = max(abs(11 - abs(speed)), 1) expon = mod(tmp + 1, 3) + 1 self.coveps = self.abseps * ((1.0e-1) ** expon) elif speed < 13: self.coveps = 0.1 else: self.coveps = 0.5 self.releps = min(self.abseps, 1.0e-2) if self.method == 0: # This gives approximately the same accuracy as when using # RINDDND and RINDNIT # xCutOff= MIN(MAX(xCutOff+0.5d0,4.d0),5.d0) self.abseps = self.abseps * 1.0e-1 trunc_error = 0.05 * max(0, self.abseps) self.xcutoff = max(min(abs(invnorm(trunc_error)), 7), 1.2) self.abseps = max(self.abseps - trunc_error, 0) def set_constants(self): if self.xcutoff is None: trunc_error = 0.1 * self.abseps self.nc1c2 = max(1, self.nc1c2) xcut = abs(invnorm(trunc_error / (self.nc1c2 * 2))) self.xcutoff = max(min(xcut, 8.5), 1.2) # self.abseps = max(self.abseps- truncError,0); # self.releps = max(self.releps- truncError,0); if self.method > 0: names = ['method', 'xcscale', 'abseps', 'releps', 'coveps', 'maxpts', 'minpts', 'nit', 'xcutoff', 'nc1c2', 'quadno', 'xsplit'] constants = [getattr(self, name) for name in names] constants[0] = mod(constants[0], 10) rindmod.set_constants(*constants) # @UndefinedVariable def __call__(self, cov, m, ab, bb, indI=None, xc=None, nt=None, **kwds): if any(kwds): self.__dict__.update(**kwds) self.set_constants() if xc is None: xc = zeros((0, 1)) big, b_lo, b_up, xc = atleast_2d(cov, ab, bb, xc) b_lo = b_lo.copy() b_up = b_up.copy() Ntdc = big.shape[0] Nc = xc.shape[0] if nt is None: nt = Ntdc - Nc unused_Mb, Nb = b_lo.shape Nd = Ntdc - nt - Nc Ntd = nt + Nd if indI is None: if Nb != Ntd: raise ValueError('Inconsistent size of b_lo and b_up') indI = r_[-1:Ntd] Ex, indI = atleast_1d(m, indI) if self.seed is None: seed = int(floor(random.rand(1) * 1e10)) # @UndefinedVariable else: seed = int(self.seed) # INFIN = INTEGER, array of integration limits flags: size 1 x Nb # if INFIN(I) < 0, Ith limits are (-infinity, infinity); # if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)]; # if INFIN(I) = 1, Ith limits are [Hlo(I), infinity); # if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)]. infinity = 37 dev = sqrt(diag(big)) # std ind = nonzero(indI[1:] > -1)[0] infin = repeat(2, len(indI) - 1) infin[ind] = (2 - (b_up[0, ind] > infinity * dev[indI[ind + 1]]) - 2 * (b_lo[0, ind] < -infinity * dev[indI[ind + 1]])) b_up[0, ind] = minimum(b_up[0, ind], infinity * dev[indI[ind + 1]]) b_lo[0, ind] = maximum(b_lo[0, ind], -infinity * dev[indI[ind + 1]]) ind2 = indI + 1 return rindmod.rind(big, Ex, xc, nt, ind2, b_lo, b_up, infin, seed) def test_rind(): """ Small test function """ n = 5 Blo = -inf Bup = -1.2 indI = [-1, n - 1] # Barriers # A = np.repeat(Blo, n) # B = np.repeat(Bup, n) # Integration limits m = zeros(n) rho = 0.3 Sc = (ones((n, n)) - eye(n)) * rho + eye(n) rind = Rind() E0 = rind(Sc, m, Blo, Bup, indI) # exact prob. 0.001946 A) print(E0) A = repeat(Blo, n) B = repeat(Bup, n) # Integration limits _E1 = rind(triu(Sc), m, A, B) # same as E0 xc = zeros((0, 1)) infinity = 37 dev = sqrt(diag(Sc)) # std ind = nonzero(indI[1:])[0] Bup, Blo = atleast_2d(Bup, Blo) Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]]) Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]]) _E3 = rind(Sc, m, Blo, Bup, indI, xc, nt=1) def cdflomax(x, alpha, m0): """ Return CDF for local maxima for a zero-mean Gaussian process Parameters ---------- x : array-like evaluation points alpha, m0 : real scalars irregularity factor and zero-order spectral moment (variance of the process), respectively. Returns ------- prb : ndarray distribution function evaluated at x Notes ----- The cdf is calculated from an explicit expression involving the standard-normal cdf. This relation is sometimes written as a convolution M = sqrt(m0)*( sqrt(1-a^2)*Z + a*R ) where M denotes local maximum, Z is a standard normal r.v., R is a standard Rayleigh r.v., and "=" means equality in distribution. Note that all local maxima of the process are considered, not only crests of waves. Example ------- >>> import pylab >>> import wafo.gaussian as wg >>> import wafo.spectrum.models as wsm >>> import wafo.objects as wo >>> import wafo.stats as ws >>> S = wsm.Jonswap(Hm0=10).tospecdata(); >>> xs = S.sim(10000) >>> ts = wo.mat2timeseries(xs) >>> tp = ts.turning_points() >>> mM = tp.cycle_pairs() >>> m0 = S.moment(1)[0] >>> alpha = S.characteristic('alpha')[0] >>> x = np.linspace(-10,10,200); >>> mcdf = ws.edf(mM.data) h = mcdf.plot(), pylab.plot(x,wg.cdflomax(x,alpha,m0)) See also -------- spec2mom, spec2bw """ c1 = 1.0 / (sqrt(1 - alpha ** 2)) * x / sqrt(m0) c2 = alpha * c1 return cdfnorm(c1) - alpha * exp(-x ** 2 / 2 / m0) * cdfnorm(c2) def prbnormtndpc(rho, a, b, d=None, df=0, abseps=1e-4, ierc=0, hnc=0.24): """ Return Multivariate normal or T probability with product correlation. Parameters ---------- rho : array-like vector of coefficients defining the correlation coefficient by: correlation(I,J) = rho[i]*rho[j]) for J!=I where -1 < rho[i] < 1 a,b : array-like vector of lower and upper integration limits, respectively. Note: any values greater the 37 in magnitude, are considered as infinite values. d : array-like vector of means (default zeros(size(rho))) df = Degrees of freedom, NDF<=0 gives normal probabilities (default) abseps = absolute error tolerance. (default 1e-4) ierc = 1 if strict error control based on fourth derivative 0 if error control based on halving the intervals (default) hnc = start interval width of simpson rule (default 0.24) Returns ------- value = estimated value for the integral bound = bound on the error of the approximation inform = INTEGER, termination status parameter: 0, if normal completion with ERROR < EPS; 1, if N > 1000 or N < 1. 2, IF any abs(rho)>=1 4, if ANY(b(I)<=A(i)) 5, if number of terms exceeds maximum number of evaluation points 6, if fault accurs in normal subroutines 7, if subintervals are too narrow or too many 8, if bounds exceeds abseps PRBNORMTNDPC calculates multivariate normal or student T probability with product correlation structure for rectangular regions. The accuracy is as best around single precision, i.e., about 1e-7. Example: -------- >>> import wafo.gaussian as wg >>> rho2 = np.random.rand(2) >>> a2 = np.zeros(2) >>> b2 = np.repeat(np.inf,2) >>> [val2,err2, ift2] = wg.prbnormtndpc(rho2,a2,b2) >>> def g2(x): ... return 0.25+np.arcsin(x[0]*x[1])/(2*np.pi) >>> E2 = g2(rho2) # exact value >>> np.abs(E2-val2)>> rho3 = np.random.rand(3) >>> a3 = np.zeros(3) >>> b3 = np.repeat(inf,3) >>> [val3, err3, ift3] = wg.prbnormtndpc(rho3,a3,b3) >>> def g3(x): ... return 0.5-sum(np.sort(np.arccos([x[0]*x[1], ... x[0]*x[2],x[1]*x[2]])))/(4*np.pi) >>> E3 = g3(rho3) # Exact value >>> np.abs(E3-val3) < 5 * err2 True See also -------- prbnormndpc, prbnormnd, Rind Reference --------- Charles Dunnett (1989) "Multivariate normal probability integrals with product correlation structure", Applied statistics, Vol 38,No3, (Algorithm AS 251) """ if d is None: d = zeros(len(rho)) # Make sure integration limits are finite aa = np.clip(a - d, -100, 100) bb = np.clip(b - d, -100, 100) return mvnprdmod.prbnormtndpc(rho, aa, bb, df, abseps, ierc, hnc) def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True, usebreakpoints=False): """ Return Multivariate Normal probabilities with product correlation Parameters ---------- rho = vector defining the correlation structure, i.e., corr(Xi,Xj) = rho(i)*rho(j) for i~=j = 1 for i==j -1 <= rho <= 1 a,b = lower and upper integration limits respectively. tol = requested absolute tolerance Returns ------- value = value of integral error = estimated absolute error PRBNORMNDPC calculates multivariate normal probability with product correlation structure for rectangular regions. The accuracy is up to almost double precision, i.e., about 1e-14. Example: ------- >>> import wafo.gaussian as wg >>> rho2 = np.random.rand(2) >>> a2 = np.zeros(2) >>> b2 = np.repeat(np.inf,2) >>> [val2,err2, ift2] = wg.prbnormndpc(rho2,a2,b2) >>> g2 = lambda x : 0.25+np.arcsin(x[0]*x[1])/(2*np.pi) >>> E2 = g2(rho2) #% exact value >>> np.abs(E2-val2)>> rho3 = np.random.rand(3) >>> a3 = np.zeros(3) >>> b3 = np.repeat(inf,3) >>> [val3,err3, ift3] = wg.prbnormndpc(rho3,a3,b3) >>> def g3(x): ... return 0.5-sum(np.sort(np.arccos([x[0]*x[1], ... x[0]*x[2],x[1]*x[2]])))/(4*np.pi) >>> E3 = g3(rho3) # Exact value >>> np.abs(E3-val3) 0: warnings.warn('Abnormal termination ier = %d\n\n%s' % (ier, _ERRORMESSAGE[ier])) return val, err, ier _ERRORMESSAGE = {} _ERRORMESSAGE[0] = '' _ERRORMESSAGE[1] = """ Maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yields no improvement it is advised to analyze the integrand in order to determine the integration difficulties. if the position of a local difficulty can be determined (i.e. singularity discontinuity within the interval), it should be supplied to the routine as an element of the vector points. If necessary an appropriate special-purpose integrator must be used, which is designed for handling the type of difficulty involved. """ _ERRORMESSAGE[2] = """ the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be under-estimated.""" _ERRORMESSAGE[3] = """ Extremely bad integrand behaviour occurs at some points of the integration interval.""" _ERRORMESSAGE[4] = """ The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained. """ _ERRORMESSAGE[5] = """ The integral is probably divergent, or slowly convergent. It must be noted that divergence can occur with any other value of ier>0. """ _ERRORMESSAGE[6] = """the input is invalid because: 1) npts2 < 2 2) break points are specified outside the integration range 3) (epsabs<=0 and epsrel EPS and MAXPTS function vaules used; increase MAXPTS to decrease ERROR; if INFORM = 2, N > NMAX or N < 1. where NMAX depends on the integration method Example ------- Compute the probability that X1<0,X2<0,X3<0,X4<0,X5<0, Xi are zero-mean Gaussian variables with variances one and correlations Cov(X(i),X(j))=0.3: indI=[0 5], and barriers B_lo=[-inf 0], B_lo=[0 inf] gives H_lo = [-inf -inf -inf -inf -inf] H_lo = [0 0 0 0 0] >>> Et = 0.001946 # # exact prob. >>> n = 5; nt = n >>> Blo =-np.inf; Bup=0; indI=[-1, n-1] # Barriers >>> m = 1.2*np.ones(n); rho = 0.3; >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n) >>> rind = Rind() >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI, nt=nt) >>> A = np.repeat(Blo,n) >>> B = np.repeat(Bup,n)-m >>> val, err, inform = prbnormnd(Sc,A,B) >>> np.allclose([val, err, inform], ... [0.0019456719705212067, 1.0059406844578488e-05, 0]) True >>> np.allclose(np.abs(val-Et) < err0+terr0, True) True >>> 'val = %2.5f' % val 'val = 0.00195' See also -------- prbnormndpc, Rind """ m, n = correl.shape Na = len(a) Nb = len(b) if (m != n or m != Na or m != Nb): raise ValueError('Size of input is inconsistent!') if maxpts is None: maxpts = 1000 * n maxpts = max(round(maxpts), 10 * n) # % 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. D = np.diag(correl) if (any(D != 1)): raise ValueError('This is not a correlation matrix') # Make sure integration limits are finite A = np.clip(a, -100, 100) B = np.clip(b, -100, 100) ix = np.where(np.triu(np.ones((m, m)), 1) != 0) L = correl[ix].ravel() # % return only off diagonal elements infinity = 37 infin = np.repeat(2, n) - (B > infinity) - 2 * (A < -infinity) err, val, inform = mvn.mvndst(A, B, infin, L, maxpts, abseps, releps) return val, err, inform # CALL the mexroutine # t0 = clock; # if ((method==0) && (n<=100)), # %NMAX = 100 # [value, err,inform] = mexmvnprb(L,A,B,abseps,releps,maxpts); # elseif ( (method<0) || ((method<=0) && (n>100)) ), # % NMAX = 500 # [value, err,inform] = mexmvnprb2(L,A,B,abseps,releps,maxpts); # else # [value, err,inform] = mexGenzMvnPrb(L,A,B,abseps,releps,maxpts,method); # end # exTime = etime(clock,t0); # ' # gauss legendre points and weights, n = 6 _W6 = [0.1713244923791705e+00, 0.3607615730481384e+00, 0.4679139345726904e+00] _X6 = [-0.9324695142031522e+00, - 0.6612093864662647e+00, -0.2386191860831970e+00] # gauss legendre points and weights, n = 12 _W12 = [0.4717533638651177e-01, 0.1069393259953183e+00, 0.1600783285433464e+00, 0.2031674267230659e+00, 0.2334925365383547e+00, 0.2491470458134029e+00] _X12 = [-0.9815606342467191e+00, -0.9041172563704750e+00, -0.7699026741943050e+00, - 0.5873179542866171e+00, -0.3678314989981802e+00, -0.1252334085114692e+00] # gauss legendre points and weights, n = 20 _W20 = [0.1761400713915212e-01, 0.4060142980038694e-01, 0.6267204833410906e-01, 0.8327674157670475e-01, 0.1019301198172404e+00, 0.1181945319615184e+00, 0.1316886384491766e+00, 0.1420961093183821e+00, 0.1491729864726037e+00, 0.1527533871307259e+00] _X20 = [-0.9931285991850949e+00, -0.9639719272779138e+00, - 0.9122344282513259e+00, -0.8391169718222188e+00, - 0.7463319064601508e+00, -0.6360536807265150e+00, - 0.5108670019508271e+00, -0.3737060887154196e+00, - 0.2277858511416451e+00, -0.7652652113349733e-01] def cdfnorm2d(b1, b2, r): """ Returnc Bivariate Normal cumulative distribution function Parameters ---------- b1, b2 : array-like upper integration limits r : real scalar correlation coefficient (-1 <= r <= 1). Returns ------- bvn : ndarray distribution function evaluated at b1, b2. Notes ----- CDFNORM2D computes bivariate normal probabilities, i.e., the probability Prob(X1 <= B1 and X2 <= B2) with an absolute error less than 1e-15. This function is based on the method described by Drezner, z and G.O. Wesolowsky, (1989), with major modifications for double precision, and for |r| close to 1. Example ------- >>> import wafo.gaussian as wg >>> x = np.linspace(-5,5,20) >>> [B1,B2] = np.meshgrid(x, x) >>> r = 0.3; >>> F = wg.cdfnorm2d(B1,B2,r) surf(x,x,F) See also -------- cdfnorm Reference --------- Drezner, z and g.o. Wesolowsky, (1989), "On the computation of the bivariate normal integral", Journal of statist. comput. simul. 35, pp. 101-107, """ # Translated into Python # Per A. Brodtkorb # # Original code # by alan genz # department of mathematics # washington state university # pullman, wa 99164-3113 # email : alangenz@wsu.edu b1, b2, r = np.broadcast_arrays(b1, b2, r) cshape = b1.shape h, k, r = -b1.ravel(), -b2.ravel(), r.ravel() bvn = where(abs(r) > 1, nan, 0.0) two = 2.e0 twopi = 6.283185307179586e0 hk = h * k k0, = nonzero(abs(r) < 0.925e0) if len(k0) > 0: hs = (h[k0] ** 2 + k[k0] ** 2) / two asr = arcsin(r[k0]) k1, = nonzero(r[k0] >= 0.75) if len(k1) > 0: k01 = k0[k1] for i in range(10): for sign in - 1, 1: sn = sin(asr[k1] * (sign * _X20[i] + 1) / 2) bvn[k01] = bvn[k01] + _W20[i] * \ exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn)) k1, = nonzero((0.3 <= r[k0]) & (r[k0] < 0.75)) if len(k1) > 0: k01 = k0[k1] for i in range(6): for sign in - 1, 1: sn = sin(asr[k1] * (sign * _X12[i] + 1) / 2) bvn[k01] = bvn[k01] + _W12[i] * \ exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn)) k1, = nonzero(r[k0] < 0.3) if len(k1) > 0: k01 = k0[k1] for i in range(3): for sign in - 1, 1: sn = sin(asr[k1] * (sign * _X6[i] + 1) / 2) bvn[k01] = bvn[k01] + _W6[i] * \ exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn)) bvn[k0] *= asr / (two * twopi) bvn[k0] += fi(-h[k0]) * fi(-k[k0]) k1, = nonzero((0.925 <= abs(r)) & (abs(r) <= 1)) if len(k1) > 0: k2, = nonzero(r[k1] < 0) if len(k2) > 0: k12 = k1[k2] k[k12] = -k[k12] hk[k12] = -hk[k12] k3, = nonzero(abs(r[k1]) < 1) if len(k3) > 0: k13 = k1[k3] a2 = (1 - r[k13]) * (1 + r[k13]) a = sqrt(a2) b = abs(h[k13] - k[k13]) bs = b * b c = (4.e0 - hk[k13]) / 8.e0 d = (12.e0 - hk[k13]) / 16.e0 asr = -(bs / a2 + hk[k13]) / 2.e0 k4, = nonzero(asr > -100.e0) if len(k4) > 0: bvn[k13[k4]] = (a[k4] * exp(asr[k4]) * (1 - c[k4] * (bs[k4] - a2[k4]) * (1 - d[k4] * bs[k4] / 5) / 3 + c[k4] * d[k4] * a2[k4] ** 2 / 5)) k5, = nonzero(hk[k13] < 100.e0) if len(k5) > 0: # b = sqrt(bs); k135 = k13[k5] bvn[k135] = bvn[k135] - exp(-hk[k135] / 2) * sqrt(twopi) * \ fi(-b[k5] / a[k5]) * b[k5] * \ (1 - c[k5] * bs[k5] * (1 - d[k5] * bs[k5] / 5) / 3) a /= two for i in range(10): for sign in - 1, 1: xs = (a * (sign * _X20[i] + 1)) ** 2 rs = sqrt(1 - xs) asr = -(bs / xs + hk[k13]) / 2 k6, = nonzero(asr > -100.e0) if len(k6) > 0: k136 = k13[k6] bvn[k136] += (a[k6] * _W20[i] * exp(asr[k6]) * (exp(-hk[k136] * (1 - rs[k6]) / (2 * (1 + rs[k6]))) / rs[k6] - (1 + c[k6] * xs[k6] * (1 + d[k6] * xs[k6])))) bvn[k3] = -bvn[k3] / twopi k7, = nonzero(r[k1] > 0) if len(k7): k17 = k1[k7] bvn[k17] += fi(-np.maximum(h[k17], k[k17])) k8, = nonzero(r[k1] < 0) if len(k8) > 0: k18 = k1[k8] bvn[k18] = -bvn[k18] + np.maximum(0, fi(-h[k18]) - fi(-k[k18])) bvn.shape = cshape return bvn def fi(x): return 0.5 * (erfc((-x) / sqrt(2))) def prbnorm2d(a, b, r): """ Returns Bivariate Normal probability Parameters --------- a, b : array-like, size==2 vector with lower and upper integration limits, respectively. r : real scalar correlation coefficient Returns ------- prb : real scalar computed probability Prob(A[0] <= X1 <= B[0] and A[1] <= X2 <= B[1]) with an absolute error less than 1e-15. Example ------- >>> import wafo.gaussian as wg >>> a = [-1, -2] >>> b = [1, 1] >>> r = 0.3 >>> np.allclose(wg.prbnorm2d(a,b,r), 0.56659121350428077) True See also -------- cdfnorm2d, cdfnorm, prbnormndpc """ infinity = 37 lower = np.asarray(a) upper = np.asarray(b) if np.all((lower <= -infinity) & (infinity <= upper)): return 1.0 if (lower >= upper).any(): return 0.0 correl = r infin = np.repeat(2, 2) - (upper > infinity) - 2 * (lower < -infinity) if np.all(infin == 2): return (bvd(lower[0], lower[1], correl) - bvd(upper[0], lower[1], correl) - bvd(lower[0], upper[1], correl) + bvd(upper[0], upper[1], correl)) elif (infin[0] == 2 and infin[1] == 1): return (bvd(lower[0], lower[1], correl) - bvd(upper[0], lower[1], correl)) elif (infin[0] == 1 and infin[1] == 2): return (bvd(lower[0], lower[1], correl) - bvd(lower[0], upper[1], correl)) elif (infin[0] == 2 and infin[1] == 0): return (bvd(-upper[0], -upper[1], correl) - bvd(-lower[0], -upper[1], correl)) elif (infin[0] == 0 and infin[1] == 2): return (bvd(-upper[0], -upper[1], correl) - bvd(-upper[0], -lower[1], correl)) elif (infin[0] == 1 and infin[1] == 0): return bvd(lower[0], -upper[1], -correl) elif (infin[0] == 0 and infin[1] == 1): return bvd(-upper[0], lower[1], -correl) elif (infin[0] == 1 and infin[1] == 1): return bvd(lower[0], lower[1], correl) elif (infin[0] == 0 and infin[1] == 0): return bvd(-upper[0], -upper[1], correl) return 1 def bvd(lo, up, r): return cdfnorm2d(-lo, -up, r) if __name__ == '__main__': from wafo.testing import test_docstrings test_docstrings(__file__) # if __name__ == '__main__': # if False: #True: # # test_rind() # else: # import doctest # doctest.testmod()