Home > wafo > trgauss > mvnormprb.m

# mvnormprb

## PURPOSE

Multivariate Normal probability by Genz' algorithm.

## SYNOPSIS

[value,err,inform,exTime] = mvnormprb(correl,A,B,abseps,releps,maxpts,method);

## DESCRIPTION

``` MVNORMPRB Multivariate Normal probability by Genz' algorithm.

CALL [value,error,inform]=mvnormprb(correl,A,B,abseps,releps,maxpts,method);

VALUE  REAL estimated value for the integral
ERROR  REAL estimated absolute error, with 99% confidence level.
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 > NMAX or N < 1. where NMAX depends on the
integration method

CORREL = Positive semidefinite correlation matrix
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
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)
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

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]

N = 5; rho=0.3; NIT=3; Nt=N; indI=[0 N];
B_lo=-10; B_up=0; m=1.2*ones(N,1);
Sc=(ones(N)-eye(N))*rho+eye(N);
E = rind(Sc,m,[],Nt,NIT,[],indI,B_lo,B_up) % exact prob. 0.00195
A = [-inf -inf -inf -inf -inf],
B = [0 0 0 0 0]-m'
[val,err,inform] = mvnormprb(Sc,A,B);

## CROSS-REFERENCE INFORMATION

This function calls:
 mexGenzMvnPrb Computes multivariate normal probability by Genz' algorithm mexmvnprb Computes multivariate normal probability by Genz' algorithm mexmvnprb2 Computes multivariate normal probability by Genz' algorithm clock Current date and time as date vector. error Display message and abort function. etime Elapsed time. tril Extract lower triangular part.
This function is called by:

## SOURCE CODE

```001 function [value,err,inform,exTime] = mvnormprb(correl,A,B,abseps,releps,maxpts,method);
002 %MVNORMPRB Multivariate Normal probability by Genz' algorithm.
003 %
004 %  CALL [value,error,inform]=mvnormprb(correl,A,B,abseps,releps,maxpts,method);
005 %
006 %     VALUE  REAL estimated value for the integral
007 %     ERROR  REAL estimated absolute error, with 99% confidence level.
008 %     INFORM INTEGER, termination status parameter:
009 %            if INFORM = 0, normal completion with ERROR < EPS;
010 %            if INFORM = 1, completion with ERROR > EPS and MAXPTS
011 %                           function vaules used; increase MAXPTS to
012 %                           decrease ERROR;
013 %            if INFORM = 2, N > NMAX or N < 1. where NMAX depends on the
014 %                           integration method
015 %
016 %     CORREL = Positive semidefinite correlation matrix
017 %     A         = vector of lower integration limits.
018 %     B         = vector of upper integration limits.
019 %     ABSEPS = absolute error tolerance.
020 %     RELEPS = relative error tolerance.
021 %     MAXPTS = maximum number of function values allowed. This
022 %              parameter can be used to limit the time. A sensible
023 %              strategy is to start with MAXPTS = 1000*N, and then
024 %              increase MAXPTS if ERROR is too large.
025 %     METHOD = integer defining the integration method
026 %             -1 KRBVRC randomized Korobov rules for the first 20
027 %                variables, randomized Richtmeyer rules for the rest,
028 %                NMAX = 500
029 %              0 KRBVRC, NMAX = 100 (default)
031 %              2 KROBOV Randomized KOROBOV rules,              NMAX = 100
032 %              3 RCRUDE Crude Monte-Carlo Algorithm with simple
033 %                antithetic variates and weighted results on restart
034 %              4 SPHMVN Monte-Carlo algorithm by Deak (1980),  NMAX = 100
035 %
036 % Example:% Compute the probability that X1<0,X2<0,X3<0,X4<0,X5<0,
037 %           % Xi are zero-mean Gaussian variables with variances one
038 %           % and correlations Cov(X(i),X(j))=0.3:
039 %           % indI=[0 5], and barriers B_lo=[-inf 0], B_lo=[0  inf]
040 %           % gives H_lo = [-inf -inf -inf -inf -inf]  H_lo = [0 0 0 0 0]
041 %
042 %    N = 5; rho=0.3; NIT=3; Nt=N; indI=[0 N];
043 %    B_lo=-10; B_up=0; m=1.2*ones(N,1);
044 %    Sc=(ones(N)-eye(N))*rho+eye(N);
045 %    E = rind(Sc,m,[],Nt,NIT,[],indI,B_lo,B_up) % exact prob. 0.00195
046 %    A = [-inf -inf -inf -inf -inf],
047 %    B = [0 0 0 0 0]-m'
048 %    [val,err,inform] = mvnormprb(Sc,A,B);
049 %
051
052 %History
053 % By pab 2002
054 error(nargchk(3,7,nargin))
055 [m,n] = size(correl);
056 Na = length(A);
057 Nb = length(B);
058 if (m~=n | m~=Na | m~=Nb)
059    error('Size of input is inconsistent!')
060 end
061
062 if nargin<4 | isempty(abseps), abseps = 1e-4; end
063 if nargin<5 | isempty(releps), releps = 1e-3; end
064 if nargin<6 | isempty(maxpts), maxpts = 1000*n; end
065 if nargin<7 | isempty(method), method = 0; end
066
067 maxpts = max(round(maxpts),10*n);
068
069 %            array of correlation coefficients; the correlation
070 %            coefficient in row I column J of the correlation matrix
071 %            should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
072 %            The correlation matrix must be positive semidefinite.
073
074 D = diag(correl);
075 if (any(D~=1))
076    error('This is not a correlation matrix')
077 end
078
079 % Make sure integration limits are finite
080 A = min(max(A,-100),100);
081 B = max(min(B,100),-100);
082 L = correl(find(tril(correl,-1)));    % return only off diagonal elements
083 %CALL the mexroutine
084 t0 = clock;
085 if ((method==0) & (n<=100)),
086   %NMAX = 100
087   [value, err,inform] = mexmvnprb(L,A,B,abseps,releps,maxpts);
088 elseif ( (method<0) | ((method<=0) & (n>100)) ),
089   % NMAX = 500
090   [value, err,inform] = mexmvnprb2(L,A,B,abseps,releps,maxpts);
091 else
092   [value, err,inform] = mexGenzMvnPrb(L,A,B,abseps,releps,maxpts,method);
093 end
094 exTime = etime(clock,t0);
095
096
097 return
098 if m>100 | m<1
099    value = 0;
100    err = 1;
101    inform = 2;
102    return
103 end```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003