Home > wafo > trgauss > mvnormpcprb.m

# mvnormpcprb

## PURPOSE

Multivariate Normal probabilities with product correlation

## SYNOPSIS

[val, err,ier,extime]= mvnormpcprb(rho,a,b,tol,useSimpson,useBreakPoints)

## DESCRIPTION

``` 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

## CROSS-REFERENCE INFORMATION

This function calls:
 mvnprodcorrprbmex Computes multivariate normal probability wnormcdf Normal cumulative distribution function wnorminv Inverse of the Normal distribution function adaptivesimpson2 clock Current date and time as date vector. dqagpe QAGPE Automatic integrator, general-purpose, etime Elapsed time. mvnprodcorrprbmex53 sprintf Write formatted data to string. str2num Convert string matrix to numeric array. unique Set unique. version MATLAB version number. warning Display warning message; disable or enable warning messages.
This function is called by:

## SOURCE CODE

```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 %
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)
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```

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