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 
  
  See also  mvnortpcprb, mvnormprb, rind

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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

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