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
Computes multivariate normal probability by Genz' algorithm | |
Computes multivariate normal probability by Genz' algorithm | |
Computes multivariate normal probability by Genz' algorithm | |
Current date and time as date vector. | |
Display message and abort function. | |
Elapsed time. | |
Extract lower triangular part. |
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
Comments or corrections to the WAFO group