MVNORMPCPRB Multivariate Normal probabilities with product correlation CALL: [value, error] = mvnormpcprb(rho,a,b,tol) value = value of integral error = estimated absolute error 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 MVNORMPCPRB 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: rho2 = rand(1,2); a2 = zeros(1,2); b2 = repmat(inf,1,2); tol = 1e-4; [val2,err2] = mvnormpcprb(rho2,a2,b2,tol); g2 = inline('0.25+asin(x(1)*x(2))/(2*pi)'); E2 = g2(rho2) % exact value rho3 = rand(1,3); a3 = zeros(1,3); b3 = repmat(inf,1,3); tol = 1e-4; [val3,err] = mvnormpcprb(rho3,a3,b3,tol); g3 = inline('0.5-sum(sort(acos([x(1)*x(2),x(1)*x(3),x(2)*x(3)])))/(4*pi)'); E3 = g3(rho3) % Exact value See also mvnortpcprb, mvnormprb, rind
Computes multivariate normal probability | |
Normal cumulative distribution function | |
Inverse of the Normal distribution function | |
Current date and time as date vector. | |
QAGPE Automatic integrator, general-purpose, | |
Elapsed time. | |
Write formatted data to string. | |
Convert string matrix to numeric array. | |
Set unique. | |
MATLAB version number. | |
Display warning message; disable or enable warning messages. |
001 function [val, err,ier,extime]= mvnormpcprb(rho,a,b,tol,useSimpson,useBreakPoints) 002 %MVNORMPCPRB Multivariate Normal probabilities with product correlation 003 % 004 % CALL: [value, error] = mvnormpcprb(rho,a,b,tol) 005 % 006 % value = value of integral 007 % error = estimated absolute error 008 % rho = vector defining the correlation structure, i.e., 009 % corr(Xi,Xj) = rho(i)*rho(j) for i~=j 010 % = 1 for i==j 011 % -1 <= rho <= 1 012 % a,b = lower and upper integration limits respectively. 013 % tol = requested absolute tolerance 014 % 015 % MVNORMPCPRB calculates multivariate normal probability 016 % with product correlation structure for rectangular regions. 017 % The accuracy is up to almost double precision, i.e., about 1e-14. 018 % 019 % Example: 020 % rho2 = rand(1,2); 021 % a2 = zeros(1,2); 022 % b2 = repmat(inf,1,2); 023 % tol = 1e-4; 024 % [val2,err2] = mvnormpcprb(rho2,a2,b2,tol); 025 % g2 = inline('0.25+asin(x(1)*x(2))/(2*pi)'); 026 % E2 = g2(rho2) % exact value 027 % 028 % rho3 = rand(1,3); 029 % a3 = zeros(1,3); 030 % b3 = repmat(inf,1,3); 031 % tol = 1e-4; 032 % [val3,err] = mvnormpcprb(rho3,a3,b3,tol); 033 % g3 = inline('0.5-sum(sort(acos([x(1)*x(2),x(1)*x(3),x(2)*x(3)])))/(4*pi)'); 034 % E3 = g3(rho3) % Exact value 035 % 036 % See also mvnortpcprb, mvnormprb, rind 037 038 %Reference 039 % P. A. Brodtkorb (2004), 040 % Evaluating multinormal probabilites with product correlation structure. 041 % In Lund university report series 042 % and in the Dr.Ing thesis: 043 % The probability of Occurrence of dangerous Wave Situations at Sea. 044 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 045 % Trondheim, Norway. 046 047 % History 048 % by pab 10.06.2003 049 if nargin<4|isempty(tol) 050 tol = 1e-4; 051 end 052 if nargin<5 | isempty(useSimpson) 053 useSimpson = 1; 054 % ind = find(abs(rho)<1); 055 % if any(ind) 056 % if all(abs(rho(ind))<1-1e-3) 057 % useSimpson = 0; 058 % end 059 % else 060 % useSimpson = 0; 061 % end 062 end 063 if nargin<5 | isempty(useBreakPoints) 064 useBreakPoints = 0; 065 end 066 067 068 tol = tol; 069 reltol = tol; %*1e-1; 070 071 trace = 0; 072 x_trace = []; 073 y_trace = []; 074 075 a = a(:); 076 b = b(:); 077 rho = rho(:); 078 %if (useSimpson) 079 % useBreakPoints = 0; 080 %else 081 % useBreakPoints = 0; 082 %end 083 084 verstr = version; 085 ind = find(verstr=='.'); 086 versionNumber = str2num(verstr(1:ind(2)-1)); 087 t1 = clock; 088 % Call fortran implementation 089 if versionNumber<=5.3 090 [val,err,ier] = mvnprodcorrprbmex53(rho,a,b,tol,reltol,... 091 useBreakPoints,useSimpson); 092 else 093 [val,err,ier] = mvnprodcorrprbmex(rho,a,b,tol,reltol,... 094 useBreakPoints,useSimpson); 095 end 096 extime = etime(clock,t1); 097 if (ier>0) 098 warning(sprintf('Abnormal termination ier = %d',ier)) 099 disp(getErrorMessage(ier)) 100 end 101 return 102 % Old matlab calls kept just in case 103 104 if any(b<=a) 105 val = 0; 106 err = 0; 107 return 108 end 109 %remove insignificant variables 110 ind = find(((a==-inf) & (b == inf))); 111 if any(ind) 112 rho(ind) = []; 113 a(ind) = []; 114 b(ind) = []; 115 % base(ind) = []; 116 if isempty(a) 117 val = 1; 118 err = 0; 119 return 120 end 121 end 122 123 val = min(wnormcdf(b)-wnormcdf(a)); 124 if length(a)==1 125 err = eps; 126 return 127 end 128 129 den = sqrt(1-rho).*sqrt(1+rho); 130 %ind = find(base); 131 %if any(ind) 132 % r0 = rho(ind); 133 % r = abs(r0); 134 % den(ind) = sqrt(r).*sqrt(2-r); 135 % rho(ind) = (2*(r0>=0)-1).*(1-r); 136 %end 137 138 139 140 val0 = 1; 141 ind = find(rho==0); 142 if any(ind) 143 val0 = prod(wnormcdf(b(ind))-wnormcdf(a(ind))); 144 a(ind) = []; 145 b(ind) = []; 146 rho(ind)= []; 147 den(ind) = []; 148 if (length(a)==0), 149 val = val0; 150 err = eps; 151 return 152 end 153 end 154 xCutOff = abs(max(wnorminv(0.05*reltol),-37)); 155 xlo = -xCutOff; %abs(max(wnorminv(0.0001*reltol),-37)) 156 xup = -xlo; 157 158 [xlo, xup] = c1c2(xlo,xup,rho,a,b,xCutOff,den); 159 160 ind = find(abs(abs(rho)-1)<=eps); 161 if any(ind) 162 a(ind) = []; 163 b(ind) = []; 164 rho(ind)= []; 165 den(ind) = []; 166 if length(a)==0 167 val = (wnormcdf(xup)-wnormcdf(xlo))*val0; 168 err = sqrt(eps); 169 return 170 end 171 end 172 173 174 175 limits = [xlo,xup]; 176 isLimitsNarrowed = ((-7 < xlo) | (xup<7)); 177 if (isLimitsNarrowed & useBreakPoints) 178 179 %xCut = max(abs(xlo),abs(xup)) 180 xCut = 2*min(den); 181 %[xlo1,xup1] = c1c2(xlo,xup,rho,a,b,xCut,den); 182 %limits = unique([limits xlo1, xup1]); 183 [xlo1,xup1] = c1c2(xlo,xup,rho,a,b,0,den); 184 [xlo2,xup2] = c1c2(xlo,xup,rho,a,b,-xCut,den); 185 [xlo3,xup3] = c1c2(xlo,xup,rho,a,b,xCut,den); 186 limits = unique([limits xlo1, xup1 xlo2, xup2 xlo3,xup3]) 187 188 %x = linspace(xlo,xup,256); 189 %f = intfun(x,rho,a,b,den); 190 %[fMax,ixMax] = max(f); 191 % $$$ if 0 %(fMax>0), 192 % $$$ limits = unique([limits,x(ixMax)]); 193 % $$$ indLo = find(f(1:ixMax)<=0); 194 % $$$ if any(indLo) 195 % $$$ ind = find(x(indLo(end))<limits); 196 % $$$ limits = unique([x(indLo(end)) limits(ind)]); 197 % $$$ end 198 % $$$ indUp = find(f(ixMax+1:end) <=0); 199 % $$$ if any(indUp) 200 % $$$ ind = find(limits<x(ixMax+indUp(1))); 201 % $$$ limits = unique([limits(ind) x(ixMax+indUp(1))]); 202 % $$$ end 203 % $$$ end 204 end 205 206 %limits 207 %limits = unique((limits+2)-2); 208 %ix = find(sqrt(eps)*10 < diff(limits)); 209 %limits = [limits(1) limits(ix+1)] 210 %limits 211 if (useSimpson) 212 [val,err,ier] = adaptivesimpson2('intfun', xlo,xup,reltol,rho,a,b,den); 213 else 214 maxSubdivisions = 100; 215 [val, err,neval,ier,alist,blist,rlist,elist,pts,iord, ... 216 level,ndin, last] = ... 217 dqagpe('intfun',xlo,xup,limits(2:end-1),reltol,0,... 218 maxSubdivisions,rho,a,b,den); 219 end 220 221 val = val*val0; 222 err = err*val0; 223 224 return 225 226 227 228 function [xMin,xMax ] = c1c2(xMin,xMax,rho,a,b,xCutOff,den) 229 %C1C2 uses the regression equation to limit the integration limits 230 % 231 232 if nargin<7 233 den = sqrt(1-rho.^2); 234 end 235 %xCutOff = xMax; 236 k = find(rho>0); 237 if any(k) 238 xMax = max(xMin, min(xMax,min((b(k)+den(k)*xCutOff)./rho(k)))); 239 xMin = min(xMax, max(xMin,max((a(k)-den(k)*xCutOff)./rho(k)))) ; 240 end 241 k1 = find(rho<0); 242 if any(k1) 243 xMax = max(xMin,min(xMax,min((a(k1)-den(k1)*xCutOff)./rho(k1)))); 244 xMin = min(xMax,max(xMin,max((b(k1)+den(k1)*xCutOff)./rho(k1)))); 245 end 246 247 return 248 249 250 function msg = getErrorMessage(ier) 251 msg = ''; 252 switch ier 253 case 1, 254 msg = sprintf('%s\n',... 255 'maximum number of subdivisions allowed',... 256 'has been achieved. one can allow more',... 257 'subdivisions by increasing the value of',... 258 'limit (and taking the according dimension',... 259 'adjustments into account). however, if',... 260 'this yields no improvement it is advised',... 261 'to analyze the integrand in order to',... 262 'determine the integration difficulties. if',... 263 'the position of a local difficulty can be',... 264 'determined (i.e. singularity,',... 265 'discontinuity within the interval), it',... 266 'should be supplied to the routine as an',... 267 'element of the vector points. If necessary',... 268 'an appropriate special-purpose integrator',... 269 'must be used, which is designed for',... 270 'handling the type of difficulty involved.'); 271 case 2, 272 msg = sprintf('%s\n',... 273 'the occurrence of roundoff error is',... 274 'detected, which prevents the requested',... 275 'tolerance from being achieved.',... 276 'the error may be under-estimated.'); 277 case 3, 278 msg = sprintf('%s\n',... 279 'extremely bad integrand behaviour occurs',... 280 'at some points of the integration interval.'); 281 282 case 4, 283 msg = sprintf('%s\n',... 284 'the algorithm does not converge.',... 285 'roundoff error is detected in the',... 286 'extrapolation table. it is presumed that',... 287 'the requested tolerance cannot be',... 288 'achieved, and that the returned result is',... 289 'the best which can be obtained.'); 290 case 5, 291 msg = sprintf('%s\n',... 292 'the integral is probably divergent, or',... 293 'slowly convergent. It must be noted that',... 294 'divergence can occur with any other value of ier>0.'); 295 case 6, 296 msg = sprintf('%s\n',... 297 'the input is invalid because:',... 298 '1) npts2 < 2',... 299 '2) break points are specified outside the integration range',... 300 '3) (epsabs<=0 and epsrel<max(50*rel.mach.acc.,0.5d-28))',... 301 '4) limit < npts2.') 302 end 303 return
Comments or corrections to the WAFO group