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)

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. erfc Complementary error function. error Display message and abort function. nan Not-a-Number.
This function is called by:
 bvnormprb Bivariate Normal probability

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