|
|
|
|
@ -1,6 +1,7 @@
|
|
|
|
|
from numpy import pi, r_, minimum, maximum, atleast_1d, atleast_2d, mod, ones, floor, \
|
|
|
|
|
random, eye, nonzero, where, repeat, sqrt, exp, inf, diag, zeros, sin, arcsin, nan #@UnresolvedImport
|
|
|
|
|
from numpy import triu #@UnresolvedImport
|
|
|
|
|
from numpy import (pi, r_, minimum, maximum, atleast_1d, atleast_2d, mod, ones,
|
|
|
|
|
floor, random, eye, nonzero, where, repeat, sqrt, exp, inf,
|
|
|
|
|
diag, zeros, sin, arcsin, nan) # @UnresolvedImport
|
|
|
|
|
from numpy import triu # @UnresolvedImport
|
|
|
|
|
from scipy.special import ndtr as cdfnorm, ndtri as invnorm
|
|
|
|
|
from scipy.special import erfc
|
|
|
|
|
from wafo import mvn
|
|
|
|
|
@ -10,10 +11,13 @@ import wafo.rindmod as rindmod
|
|
|
|
|
import warnings
|
|
|
|
|
from wafo.misc import common_shape
|
|
|
|
|
|
|
|
|
|
__all__ = ['Rind', 'rindmod', 'mvnprdmod', 'mvn', 'cdflomax' , 'prbnormtndpc',
|
|
|
|
|
'prbnormndpc', 'prbnormnd', 'cdfnorm2d', 'prbnorm2d','cdfnorm','invnorm',
|
|
|
|
|
'test_docstring']
|
|
|
|
|
__all__ = ['Rind', 'rindmod', 'mvnprdmod', 'mvn', 'cdflomax', 'prbnormtndpc',
|
|
|
|
|
'prbnormndpc', 'prbnormnd', 'cdfnorm2d', 'prbnorm2d', 'cdfnorm',
|
|
|
|
|
'invnorm', 'test_docstring']
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
class Rind(object):
|
|
|
|
|
|
|
|
|
|
'''
|
|
|
|
|
RIND Computes multivariate normal expectations
|
|
|
|
|
|
|
|
|
|
@ -24,7 +28,8 @@ class Rind(object):
|
|
|
|
|
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.
|
|
|
|
|
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)
|
|
|
|
|
@ -50,7 +55,7 @@ class Rind(object):
|
|
|
|
|
"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.
|
|
|
|
|
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.
|
|
|
|
|
@ -126,10 +131,10 @@ class Rind(object):
|
|
|
|
|
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)
|
|
|
|
|
Methodology And Computing In Applied Probability, Volume 8, Number 1,
|
|
|
|
|
pp. 65-91(27)
|
|
|
|
|
'''
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def __init__(self, **kwds):
|
|
|
|
|
'''
|
|
|
|
|
Parameters
|
|
|
|
|
@ -149,25 +154,28 @@ class Rind(object):
|
|
|
|
|
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)
|
|
|
|
|
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.
|
|
|
|
|
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)
|
|
|
|
|
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.)
|
|
|
|
|
(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
|
|
|
|
|
@ -187,8 +195,8 @@ class Rind(object):
|
|
|
|
|
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)
|
|
|
|
|
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
|
|
|
|
|
@ -263,7 +271,7 @@ class Rind(object):
|
|
|
|
|
|
|
|
|
|
self.releps = min(self.abseps, 1.0e-2)
|
|
|
|
|
|
|
|
|
|
if self.method == 0 :
|
|
|
|
|
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)
|
|
|
|
|
@ -278,17 +286,17 @@ class Rind(object):
|
|
|
|
|
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);
|
|
|
|
|
# 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']
|
|
|
|
|
'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
|
|
|
|
|
rindmod.set_constants(*constants) # @UndefinedVariable
|
|
|
|
|
|
|
|
|
|
def __call__(self, cov, m, ab, bb, indI=None, xc=None, nt=None, **kwds):
|
|
|
|
|
if any(kwds):
|
|
|
|
|
@ -317,7 +325,7 @@ class Rind(object):
|
|
|
|
|
|
|
|
|
|
Ex, indI = atleast_1d(m, indI)
|
|
|
|
|
if self.seed is None:
|
|
|
|
|
seed = int(floor(random.rand(1) * 1e10)) #@UndefinedVariable
|
|
|
|
|
seed = int(floor(random.rand(1) * 1e10)) # @UndefinedVariable
|
|
|
|
|
else:
|
|
|
|
|
seed = int(self.seed)
|
|
|
|
|
|
|
|
|
|
@ -336,7 +344,9 @@ class Rind(object):
|
|
|
|
|
Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]])
|
|
|
|
|
Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]])
|
|
|
|
|
ind2 = indI + 1
|
|
|
|
|
return rindmod.rind(BIG, Ex, xc, nt, ind2, Blo, Bup, infin, seed) #@UndefinedVariable
|
|
|
|
|
|
|
|
|
|
return rindmod.rind(BIG, Ex, xc, nt, ind2, Blo, Bup, infin, seed) # @UndefinedVariable @IgnorePep8
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def test_rind():
|
|
|
|
|
''' Small test function
|
|
|
|
|
@ -346,17 +356,17 @@ def test_rind():
|
|
|
|
|
Bup = -1.2
|
|
|
|
|
indI = [-1, n - 1] # Barriers
|
|
|
|
|
# A = np.repeat(Blo, n)
|
|
|
|
|
# B = np.repeat(Bup, n) # Integration limits
|
|
|
|
|
# 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)
|
|
|
|
|
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
|
|
|
|
|
_E1 = rind(triu(Sc), m, A, B) # same as E0
|
|
|
|
|
|
|
|
|
|
xc = zeros((0, 1))
|
|
|
|
|
infinity = 37
|
|
|
|
|
@ -365,7 +375,8 @@ def test_rind():
|
|
|
|
|
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)
|
|
|
|
|
_E3 = rind(Sc, m, Blo, Bup, indI, xc, nt=1)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def cdflomax(x, alpha, m0):
|
|
|
|
|
'''
|
|
|
|
|
@ -423,9 +434,10 @@ def cdflomax(x, alpha, 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 structure.
|
|
|
|
|
Return Multivariate normal or T probability with product correlation.
|
|
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
|
----------
|
|
|
|
|
@ -435,13 +447,14 @@ def prbnormtndpc(rho, a, b, D=None, df=0, abseps=1e-4, IERC=0, HNC=0.24):
|
|
|
|
|
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.
|
|
|
|
|
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 intuitive error control based on halving the intervals (default)
|
|
|
|
|
0 if error control based on halving the intervals (default)
|
|
|
|
|
HNC = start interval width of simpson rule (default 0.24)
|
|
|
|
|
|
|
|
|
|
Returns
|
|
|
|
|
@ -453,7 +466,7 @@ def prbnormtndpc(rho, a, b, D=None, df=0, abseps=1e-4, IERC=0, HNC=0.24):
|
|
|
|
|
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 computed exceeds maximum number of evaluation points
|
|
|
|
|
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
|
|
|
|
|
@ -465,12 +478,12 @@ def prbnormtndpc(rho, a, b, D=None, df=0, abseps=1e-4, IERC=0, HNC=0.24):
|
|
|
|
|
Example:
|
|
|
|
|
--------
|
|
|
|
|
>>> import wafo.gaussian as wg
|
|
|
|
|
>>> rho2 = np.random.rand(2);
|
|
|
|
|
>>> a2 = np.zeros(2);
|
|
|
|
|
>>> b2 = np.repeat(np.inf,2);
|
|
|
|
|
>>> rho2 = np.random.rand(2)
|
|
|
|
|
>>> a2 = np.zeros(2)
|
|
|
|
|
>>> b2 = np.repeat(np.inf,2)
|
|
|
|
|
>>> [val2,err2, ift2] = wg.prbnormtndpc(rho2,a2,b2)
|
|
|
|
|
>>> g2 = lambda x : 0.25+np.arcsin(x[0]*x[1])/(2*pi)
|
|
|
|
|
>>> E2 = g2(rho2) #% exact value
|
|
|
|
|
>>> E2 = g2(rho2) # exact value
|
|
|
|
|
>>> np.abs(E2-val2)<err2
|
|
|
|
|
True
|
|
|
|
|
|
|
|
|
|
@ -500,9 +513,12 @@ def prbnormtndpc(rho, a, b, D=None, df=0, abseps=1e-4, IERC=0, HNC=0.24):
|
|
|
|
|
# Make sure integration limits are finite
|
|
|
|
|
A = np.clip(a - D, -100, 100)
|
|
|
|
|
B = np.clip(b - D, -100, 100)
|
|
|
|
|
return mvnprdmod.prbnormtndpc(rho, A, B, df, abseps, IERC, HNC) #@UndefinedVariable
|
|
|
|
|
|
|
|
|
|
def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True, usebreakpoints=False):
|
|
|
|
|
return mvnprdmod.prbnormtndpc(rho, A, B, df, abseps, IERC, HNC) # @UndefinedVariable @IgnorePep8
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True,
|
|
|
|
|
usebreakpoints=False):
|
|
|
|
|
'''
|
|
|
|
|
Return Multivariate Normal probabilities with product correlation
|
|
|
|
|
|
|
|
|
|
@ -527,9 +543,9 @@ def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True, usebreakpo
|
|
|
|
|
Example:
|
|
|
|
|
-------
|
|
|
|
|
>>> import wafo.gaussian as wg
|
|
|
|
|
>>> rho2 = np.random.rand(2);
|
|
|
|
|
>>> a2 = np.zeros(2);
|
|
|
|
|
>>> b2 = np.repeat(np.inf,2);
|
|
|
|
|
>>> 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*pi)
|
|
|
|
|
>>> E2 = g2(rho2) #% exact value
|
|
|
|
|
@ -561,38 +577,43 @@ def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True, usebreakpo
|
|
|
|
|
|
|
|
|
|
'''
|
|
|
|
|
# Call fortran implementation
|
|
|
|
|
val, err, ier = mvnprdmod.prbnormndpc(rho, a, b, abserr, relerr, usebreakpoints, usesimpson); #@UndefinedVariable
|
|
|
|
|
val, err, ier = mvnprdmod.prbnormndpc(rho, a, b, abserr, relerr, usebreakpoints, usesimpson) # @UndefinedVariable @IgnorePep8
|
|
|
|
|
|
|
|
|
|
if ier > 0:
|
|
|
|
|
warnings.warn('Abnormal termination ier = %d\n\n%s' % (ier, _ERRORMESSAGE[ier]))
|
|
|
|
|
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.
|
|
|
|
|
'''
|
|
|
|
|
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.'''
|
|
|
|
|
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.'''
|
|
|
|
|
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.'''
|
|
|
|
|
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.'''
|
|
|
|
|
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
|
|
|
|
|
@ -600,31 +621,30 @@ _ERRORMESSAGE[6] = '''the input is invalid because:
|
|
|
|
|
4) limit < npts2.'''
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
'''
|
|
|
|
|
Multivariate Normal probability by Genz' algorithm.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Parameters
|
|
|
|
|
CORREL = Positive semidefinite correlation matrix
|
|
|
|
|
A = vector of lower integration limits.
|
|
|
|
|
B = vector of upper integration limits.
|
|
|
|
|
A = vector of lower integration limits.
|
|
|
|
|
B = vector of upper integration limits.
|
|
|
|
|
ABSEPS = absolute error tolerance.
|
|
|
|
|
RELEPS = relative error tolerance.
|
|
|
|
|
MAXPTS = 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.
|
|
|
|
|
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.
|
|
|
|
|
METHOD = integer defining the integration method
|
|
|
|
|
-1 KRBVRC randomized Korobov rules for the first 20
|
|
|
|
|
variables, randomized Richtmeyer rules for the rest,
|
|
|
|
|
NMAX = 500
|
|
|
|
|
0 KRBVRC, NMAX = 100 (default)
|
|
|
|
|
1 SADAPT Subregion Adaptive integration method, NMAX = 20
|
|
|
|
|
2 KROBOV Randomized KOROBOV rules, NMAX = 100
|
|
|
|
|
3 RCRUDE Crude Monte-Carlo Algorithm with simple
|
|
|
|
|
antithetic variates and weighted results on restart
|
|
|
|
|
4 SPHMVN Monte-Carlo algorithm by Deak (1980), NMAX = 100
|
|
|
|
|
-1 KRBVRC randomized Korobov rules for the first 20 variables,
|
|
|
|
|
randomized Richtmeyer rules for the rest, NMAX = 500
|
|
|
|
|
0 KRBVRC, NMAX = 100 (default)
|
|
|
|
|
1 SADAPT Subregion Adaptive integration method, NMAX = 20
|
|
|
|
|
2 KROBOV Randomized KOROBOV rules, NMAX = 100
|
|
|
|
|
3 RCRUDE Crude Monte-Carlo Algorithm with simple
|
|
|
|
|
antithetic variates and weighted results on restart
|
|
|
|
|
4 SPHMVN Monte-Carlo algorithm by Deak (1980), NMAX = 100
|
|
|
|
|
Returns
|
|
|
|
|
-------
|
|
|
|
|
VALUE REAL estimated value for the integral
|
|
|
|
|
@ -667,7 +687,6 @@ def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
prbnormndpc, Rind
|
|
|
|
|
'''
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
m, n = correl.shape
|
|
|
|
|
Na = len(a)
|
|
|
|
|
Nb = len(b)
|
|
|
|
|
@ -677,7 +696,7 @@ def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
if maxpts is None:
|
|
|
|
|
maxpts = 1000 * n
|
|
|
|
|
|
|
|
|
|
maxpts = max(round(maxpts), 10 * n);
|
|
|
|
|
maxpts = max(round(maxpts), 10 * n)
|
|
|
|
|
|
|
|
|
|
# % array of correlation coefficients; the correlation
|
|
|
|
|
# % coefficient in row I column J of the correlation matrix
|
|
|
|
|
@ -692,16 +711,16 @@ def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
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
|
|
|
|
|
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) #@UndefinedVariable
|
|
|
|
|
err, val, inform = mvn.mvndst(A, B, infin, L, maxpts, abseps, releps) # @UndefinedVariable @IgnorePep8
|
|
|
|
|
|
|
|
|
|
return val, err, inform
|
|
|
|
|
|
|
|
|
|
#CALL the mexroutine
|
|
|
|
|
# CALL the mexroutine
|
|
|
|
|
# t0 = clock;
|
|
|
|
|
# if ((method==0) && (n<=100)),
|
|
|
|
|
# %NMAX = 100
|
|
|
|
|
@ -715,25 +734,29 @@ def prbnormnd(correl, a, b, abseps=1e-4, releps=1e-3, maxpts=None, method=0):
|
|
|
|
|
# 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]
|
|
|
|
|
# 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):
|
|
|
|
|
'''
|
|
|
|
|
@ -798,8 +821,8 @@ def cdfnorm2d(b1, b2, r):
|
|
|
|
|
|
|
|
|
|
bvn = where(abs(r) > 1, nan, 0.0)
|
|
|
|
|
|
|
|
|
|
two = 2.e0;
|
|
|
|
|
twopi = 6.283185307179586e0;
|
|
|
|
|
two = 2.e0
|
|
|
|
|
twopi = 6.283185307179586e0
|
|
|
|
|
|
|
|
|
|
hk = h * k
|
|
|
|
|
|
|
|
|
|
@ -813,90 +836,97 @@ def cdfnorm2d(b1, b2, r):
|
|
|
|
|
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));
|
|
|
|
|
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];
|
|
|
|
|
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));
|
|
|
|
|
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);
|
|
|
|
|
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[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));
|
|
|
|
|
k1, = nonzero((0.925 <= abs(r)) & (abs(r) <= 1))
|
|
|
|
|
if len(k1) > 0:
|
|
|
|
|
k2, = nonzero(r[k1] < 0);
|
|
|
|
|
k2, = nonzero(r[k1] < 0)
|
|
|
|
|
if len(k2) > 0:
|
|
|
|
|
k12 = k1[k2];
|
|
|
|
|
k[k12] = -k[k12];
|
|
|
|
|
hk[k12] = -hk[k12];
|
|
|
|
|
k12 = k1[k2]
|
|
|
|
|
k[k12] = -k[k12]
|
|
|
|
|
hk[k12] = -hk[k12]
|
|
|
|
|
|
|
|
|
|
k3, = nonzero(abs(r[k1]) < 1);
|
|
|
|
|
k3, = nonzero(abs(r[k1]) < 1)
|
|
|
|
|
if len(k3) > 0:
|
|
|
|
|
k13 = k1[k3];
|
|
|
|
|
a2 = (1 - r[k13]) * (1 + r[k13]);
|
|
|
|
|
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);
|
|
|
|
|
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);
|
|
|
|
|
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);
|
|
|
|
|
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)
|
|
|
|
|
# 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) ;
|
|
|
|
|
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]))))
|
|
|
|
|
(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;
|
|
|
|
|
bvn[k3] = -bvn[k3] / twopi
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
k7, = nonzero(r[k1] > 0);
|
|
|
|
|
k7, = nonzero(r[k1] > 0)
|
|
|
|
|
if len(k7):
|
|
|
|
|
k17 = k1[k7]
|
|
|
|
|
bvn[k17] += fi(-np.maximum(h[k17], k[k17]));
|
|
|
|
|
bvn[k17] += fi(-np.maximum(h[k17], k[k17]))
|
|
|
|
|
|
|
|
|
|
k8, = nonzero(r[k1] < 0);
|
|
|
|
|
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]));
|
|
|
|
|
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
|
|
|
|
|
@ -932,7 +962,7 @@ def prbnorm2d(a, b, r):
|
|
|
|
|
infinity = 37
|
|
|
|
|
lower = np.asarray(a)
|
|
|
|
|
upper = np.asarray(b)
|
|
|
|
|
if np.all((lower <= -infinity) & (infinity<=upper)):
|
|
|
|
|
if np.all((lower <= -infinity) & (infinity <= upper)):
|
|
|
|
|
return 1.0
|
|
|
|
|
if (lower >= upper).any():
|
|
|
|
|
return 0.0
|
|
|
|
|
@ -941,37 +971,43 @@ def prbnorm2d(a, b, r):
|
|
|
|
|
|
|
|
|
|
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) :
|
|
|
|
|
- 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) :
|
|
|
|
|
elif (infin[0] == 0 and infin[1] == 1):
|
|
|
|
|
return bvd(-upper[0], lower[1], -correl)
|
|
|
|
|
elif (infin[0] == 1 and infin[1] == 1):
|
|
|
|
|
elif (infin[0] == 1 and infin[1] == 1):
|
|
|
|
|
return bvd(lower[0], lower[1], correl)
|
|
|
|
|
elif (infin[0] == 0 and infin[1] == 0) :
|
|
|
|
|
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)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def test_docstrings():
|
|
|
|
|
import doctest
|
|
|
|
|
doctest.testmod()
|
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
|
test_docstrings()
|
|
|
|
|
#if __name__ == '__main__':
|
|
|
|
|
# if __name__ == '__main__':
|
|
|
|
|
# if False: #True: #
|
|
|
|
|
# test_rind()
|
|
|
|
|
# else:
|
|
|
|
|
|
Per.Andreas.Brodtkorb
10 years ago