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  
               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 
  
  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);   
  
  See also  mvnormpcprb, rind

CROSS-REFERENCE INFORMATION ^

This function calls: 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) 
030 %              1 SADAPT Subregion Adaptive integration method, NMAX = 20  
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 % 
050 % See also  mvnormpcprb, rind 
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