Home > wafo > trgauss > bvnormcdf.m

bvnormcdf

PURPOSE ^

Bivariate Normal cumulative distribution function

SYNOPSIS ^

bvn = bvnormcdf( b1, b2, r )

DESCRIPTION ^

 BVNORMCDF Bivariate Normal cumulative distribution function 
  
   CALL:  val = bvnormcdf( b1, b2, r ) 
  
     bvn   = distribution function evaluated at b1, b2.   
     b1,b2 = upper integration limits 
     r     = correlation coefficient  (-1 <= r <= 1). 
  
   BVNORMCDF computes bivariate normal probabilities, i.e., 
   the probability Prob(X1 <= B1 and X2 <= B2) with an absolute error 
   less than 1e-15. 
       
   This function is based on the method described by 
   drezner, z and g.o. wesolowsky, (1989), 
   on the computation of the bivariate normal integral, 
   journal of statist. comput. simul. 35, pp. 101-107, 
   with major modifications for double precision, and for |r| close to 1. 
  
  Example  
   x = linspace(-5,5,20);   
   [B1,B2] = meshgrid(x); 
   r  = 0.3; 
   F = bvnormcdf(B1,B2,r);   
   surf(x,x,F) 
  
  See also  wnormcdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function bvn = bvnormcdf( b1, b2, r )  
002 %BVNORMCDF Bivariate Normal cumulative distribution function 
003 % 
004 %  CALL:  val = bvnormcdf( b1, b2, r ) 
005 % 
006 %    bvn   = distribution function evaluated at b1, b2.   
007 %    b1,b2 = upper integration limits 
008 %    r     = correlation coefficient  (-1 <= r <= 1). 
009 % 
010 %  BVNORMCDF computes bivariate normal probabilities, i.e., 
011 %  the probability Prob(X1 <= B1 and X2 <= B2) with an absolute error 
012 %  less than 1e-15. 
013 %      
014 %  This function is based on the method described by 
015 %  drezner, z and g.o. wesolowsky, (1989), 
016 %  on the computation of the bivariate normal integral, 
017 %  journal of statist. comput. simul. 35, pp. 101-107, 
018 %  with major modifications for double precision, and for |r| close to 1. 
019 % 
020 % Example  
021 %  x = linspace(-5,5,20);   
022 %  [B1,B2] = meshgrid(x); 
023 %  r  = 0.3; 
024 %  F = bvnormcdf(B1,B2,r);   
025 %  surf(x,x,F) 
026 % 
027 % See also  wnormcdf 
028    
029 % Reference 
030 %  Drezner, z and g.o. Wesolowsky, (1989), 
031 %  "on the computation of the bivariate normal integral", 
032 %  journal of statist. comput. simul. 35, pp. 101-107, 
033    
034 %History   
035 % revised pab 19.01.2001 
036 %  -vectorized the code to handle multiple integration limits   
037 %  -renamed from bvnd to bvnormcdf   
038 %  -replaced call to erf with erfc   
039 %    
040 % by  alan genz 
041 %     department of mathematics 
042 %     washington state university 
043 %     pullman, wa 99164-3113 
044 %     email : alangenz@wsu.edu 
045 error(nargchk(3,3,nargin)) 
046 [errorcode,h,k,r] = comnsize(-b1,-b2,r); 
047 if (errorcode > 0) 
048   error ('b1, b2 and r must be of common size or scalar'); 
049 end 
050  
051 bvn  = zeros(size(h)); 
052  
053 zero = 0.d0; 
054 one  = 1.d0; 
055 two  = 2.d0; 
056 twopi = 6.283185307179586d0; 
057 %     gauss legendre points and weights, n = 6 
058 w6 = [ 0.1713244923791705D+00, 0.3607615730481384D+00, 0.4679139345726904D+00]; 
059 x6 = [-0.9324695142031522D+00,-0.6612093864662647D+00,-0.2386191860831970D+00]; 
060 %     gauss legendre points and weights, n = 12 
061 w12 = [ 0.4717533638651177d-01, 0.1069393259953183d+00, 0.1600783285433464d+00,... 
062        0.2031674267230659d+00, 0.2334925365383547d+00, 0.2491470458134029d+00]; 
063 x12=[ -0.9815606342467191d+00,-0.9041172563704750d+00,-0.7699026741943050d+00, ... 
064     -0.5873179542866171d+00,-0.3678314989981802d+00,-0.1252334085114692d+00]; 
065 %     gauss legendre points and weights, n = 20 
066 w20 = [ 0.1761400713915212d-01, 0.4060142980038694d-01, ... 
067       0.6267204833410906d-01,  0.8327674157670475d-01, ... 
068       0.1019301198172404d+00, 0.1181945319615184d+00,... 
069       0.1316886384491766d+00, 0.1420961093183821d+00,... 
070       0.1491729864726037d+00, 0.1527533871307259d+00]; 
071 x20=[ -0.9931285991850949d+00, -0.9639719272779138d+00, ... 
072       -0.9122344282513259d+00, -0.8391169718222188d+00, ... 
073        -0.7463319064601508d+00,-0.6360536807265150d+00,... 
074        -0.5108670019508271d+00,-0.3737060887154196d+00, ... 
075    -0.2277858511416451d+00, -0.7652652113349733d-01]; 
076  
077 hk  = h.*k; 
078  
079 k0 = find(abs(r) < 0.925d0 ); 
080 if (k0)  
081   hs = ( h(k0).^2 + k(k0).^2 )/two; 
082   asr = asin(r(k0)); 
083   k1 = find(r(k0)>0.75); 
084   if any(k1) 
085     k01 = k0(k1); 
086     for i = 1:10 
087       for is = -1:2:1,  
088     sn = sin( asr(k1).*(is.*x20(i) +1 )/2 ); 
089     bvn(k01) = bvn(k01) + w20(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) ); 
090       end  
091     end  
092   end 
093    k1 = find(0.3 <= r(k0) & r(k0)<0.75); 
094   if any(k1) 
095     k01 = k0(k1); 
096     for i = 1:6 
097       for is = -1:2:1,  
098     sn = sin( asr(k1).*(is.*x12(i) +1 )/2 ); 
099     bvn(k01) = bvn(k01) + w12(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) ); 
100       end  
101     end  
102   end 
103    k1 = find( r(k0)<0.3); 
104   if any(k1) 
105     k01 = k0(k1); 
106     for i = 1:3 
107       for is = -1:2:1,  
108     sn = sin( asr(k1).*(is.*x6(i) +1 )/2 ); 
109     bvn(k01) = bvn(k01) + w6(i)*exp( ( sn.*hk(k01) - hs(k1) )./(1 - sn.*sn ) ); 
110       end  
111     end  
112   end 
113   bvn(k0) = bvn(k0).*asr/( two*twopi ) + fi(-h(k0)).*fi(-k(k0)); 
114 end 
115  
116 k1 = find(0.925<=abs(r) & abs(r)<=1 ); 
117 if any(k1) 
118   k2 = find(r(k1) < 0); 
119   if any(k2 )  
120     k12 = k1(k2); 
121     k(k12)  = -k(k12); 
122     hk(k12) = -hk(k12); 
123   end 
124   k3 = find( abs(r(k1)) < 1); 
125   if (k3) 
126     k13 = k1(k3); 
127     as = (1 - r(k13) ).*(1 + r(k13) ); 
128     a  = sqrt(as); 
129     b  = abs( h(k13) - k(k13) ); 
130     bs = b.*b; 
131     c  = ( 4.d0 - hk(k13) )/8.d0; 
132     d  = ( 12.d0 - hk(k13) )/16.d0; 
133     asr = -( bs./as + hk(k13) )/2.d0; 
134     k4 = find(asr > -100.d0); 
135     if (k4)  
136       bvn(k13(k4)) = a(k4).*exp(asr(k4)).*(1 - c(k4).* .... 
137       ( bs(k4) - as(k4)).*(1 - d(k4).*bs(k4)/5 )/3 ... 
138       + c(k4).*d(k4).*as(k4).^2/5 ); 
139     end 
140     k5 = find(hk(k13) < 100.d0); 
141     if ( k5 ) 
142       %               b = sqrt(bs); 
143       k135 = k13(k5); 
144       bvn(k135) = bvn(k135) - exp( -hk(k135)/2 ).*sqrt(twopi).*fi(-b(k5)./a(k5)).*b(k5)... 
145       .*(1- c(k5).*bs(k5).*(1 - d(k5).*bs(k5)/5 )/3 ); 
146     end 
147     a = a/two; 
148     for i = 1:10 
149       for is = -1:2:1, 
150     xs  = ( a.*( is*x20(i) + 1 ) ).^2; 
151     rs  = sqrt(1 - xs ); 
152     asr = -( bs./xs + hk(k13) )/2; 
153     k6 = find( asr > -100.d0 ) ; 
154     if any(k6)  
155       k136 = k13(k6); 
156       bvn(k136) = bvn(k136) + a(k6).*w20(i).*exp( asr(k6) )... 
157           .*( exp(-hk(k136).*(1-rs(k6))./(2*(1+rs(k6))))./rs(k6)... 
158           - (1 + c(k6).*xs(k6).*(1 + d(k6).*xs(k6) ) ) ); 
159     end  
160       end  
161     end  
162     bvn(k3) = -bvn(k3)/twopi; 
163     %bvn = -bvn/twopi; 
164   end 
165   k7 = find(r(k1)>0); 
166   if any(k7 ), 
167     k17 = k1(k7);  
168     bvn(k17) = bvn(k17) + fi( -max( h(k17), k(k17)) ); 
169   end 
170   k8 = find(r(k1)<0); 
171   if any(k8),  
172     k18 = k1(k8); 
173     bvn(k18) = -bvn(k18) + max(0, fi(-h(k18)) - fi(-k(k18)) ); 
174   end 
175 end 
176  
177 k9 = find(abs(r)>1); 
178 if any(k9) 
179   bvn(k9) = NaN; 
180 end 
181 val = bvn; 
182 return 
183  
184 function F = fi(x) 
185 F = 0.5.*(erfc((-x)./sqrt(2))); 
186    
187

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