function p = bvtl( nu, dh, dk, r ) %BVTL % p = bvtl( nu, dh, dk, r ) % A function for computing bivariate t probabilities. % bvtl calculates the probability that x < dh and y < dk; % Parameters % nu integer number of degrees of freedom, nu < 1, gives Normal case % dh 1st upper integration limit % dk 2nd upper integration limit % r correlation coefficient % Example: p = bvtl( 6, 3, 4, .35 ) % % This function is based on the method described by % Dunnett, C.W. and M. Sobel, (1954), % A bivariate generalization of Student's t-distribution % with tables for certain special cases, % Biometrika 41, pp. 153-169, % % Alan Genz % Department of Mathematics % Washington State University % Pullman, Wa 99164-3113 % Email : alangenz@wsu.edu % if nu < 1, p = bvnl( dh, dk, r ); elseif dh == -inf | dk == -inf, p = 0; elseif dh == inf, if dk == inf, p = 1; else p = studnt( nu, dk ); end elseif dk == inf, p = studnt( nu, dh ); elseif 1 - r < eps, p = studnt( nu, min([ dh dk ]) ); elseif r + 1 < eps, p = 0; if dh > -dk, p = studnt( nu, dh ) - studnt( nu, -dk ); end else, tpi = 2*pi; ors = 1 - r*r; hrk = dh - r*dk; krh = dk - r*dh; if abs(hrk) + ors > 0 xnhk = hrk^2/( hrk^2 + ors*( nu + dk^2 ) ); xnkh = krh^2/( krh^2 + ors*( nu + dh^2 ) ); else, xnhk = 0; xnkh = 0; end, hs = sign( dh - r*dk ); ks = sign( dk - r*dh ); if mod( nu, 2 ) == 0 bvt = atan2( sqrt(ors), -r )/tpi; gmph = dh/sqrt( 16*( nu + dh^2 ) ); gmpk = dk/sqrt( 16*( nu + dk^2 ) ); btnckh = 2*atan2( sqrt( xnkh ), sqrt( 1 - xnkh ) )/pi; btpdkh = 2*sqrt( xnkh*( 1 - xnkh ) )/pi; btnchk = 2*atan2( sqrt( xnhk ), sqrt( 1 - xnhk ) )/pi; btpdhk = 2*sqrt( xnhk*( 1 - xnhk ) )/pi; for j = 1 : nu/2 bvt = bvt + gmph*( 1 + ks*btnckh ); bvt = bvt + gmpk*( 1 + hs*btnchk ); btnckh = btnckh + btpdkh; btpdkh = 2*j*btpdkh*( 1 - xnkh )/(2*j+1); btnchk = btnchk + btpdhk; btpdhk = 2*j*btpdhk*( 1 - xnhk )/(2*j+1); gmph = gmph*( j - 1/2 )/( j*( 1 + dh^2/nu ) ); gmpk = gmpk*( j - 1/2 )/( j*( 1 + dk^2/nu ) ); end else, qhrk = sqrt( dh^2 + dk^2 - 2*r*dh*dk + nu*ors ); hkrn = dh*dk + r*nu; hkn = dh*dk - nu; hpk = dh + dk; bvt = atan2( -sqrt(nu)*(hkn*qhrk+hpk*hkrn), hkn*hkrn-nu*hpk*qhrk )/tpi; if bvt < -10*eps, bvt = bvt + 1; end gmph = dh/( tpi*sqrt(nu)*( 1 + dh^2/nu ) ); gmpk = dk/( tpi*sqrt(nu)*( 1 + dk^2/nu ) ); btnckh = sqrt( xnkh ); btpdkh = btnckh; btnchk = sqrt( xnhk ); btpdhk = btnchk; for j = 1 : ( nu - 1 )/2 bvt = bvt + gmph*( 1 + ks*btnckh ); bvt = bvt + gmpk*( 1 + hs*btnchk ); btpdkh = (2*j-1)*btpdkh*( 1 - xnkh )/(2*j); btnckh = btnckh + btpdkh; btpdhk = (2*j-1)*btpdhk*( 1 - xnhk )/(2*j); btnchk = btnchk + btpdhk; gmph = gmph*j/( ( j + 1/2 )*( 1 + dh^2/nu ) ); gmpk = gmpk*j/( ( j + 1/2 )*( 1 + dk^2/nu ) ); end end, p = bvt; end % % end bvtl % function st = studnt( nu, t ) % % Student t Distribution Function; % % t % studnt = c i ( 1 + y*y/nu )^( -(nu+1)/2 ) dy % nu -inf % if t == inf, st = 1; elseif t == -inf, st = 0; elseif nu < 1, st = phid(t); elseif nu == 1, st = ( 1 + 2*atan(t)/pi )/2; elseif nu == 2, st = ( 1 + t/sqrt( 2 + t*t ))/2; else, tt = t*t; cssthe = 1/( 1 + tt/nu ); polyn = 1; for j = nu-2 : -2 : 2, polyn = 1 + ( j - 1 )*cssthe*polyn/j; end if mod( nu, 2 ) == 1, rn = nu; ts = t/sqrt(rn); st = ( 1 + 2*( atan(ts) + ts*cssthe*polyn )/pi )/2; else, snthe = t/sqrt( nu + tt ); st = ( 1 + snthe*polyn )/2; end st = max( [ 0 min( [st 1] )] ); end % % end studnt % function p = bvnl( dh, dk, r ) %BVNL p = bvnu( -dh, -dk, r ); % % end bvnl % function p = bvnu( dh, dk, r ) %BVNU % A function for computing bivariate normal probabilities. % bvnu calculates the probability that x > dh and y > dk. % parameters % dh 1st lower integration limit % dk 2nd lower integration limit % r correlation coefficient % Example: p = bvnu( -3, -1, .35 ) % Note: to compute the probability that x < dh and y < dk, % use bvnu( -dh, -dk, r ). % % Author % Alan Genz % Department of Mathematics % Washington State University % Pullman, Wa 99164-3113 % Email : alangenz@wsu.edu % % This function is based on the method described by % Drezner, Z and G.O. Wesolowsky, (1989), % On the computation of the bivariate normal inegral, % Journal of Statist. Comput. Simul. 35, pp. 101-107, % with major modifications for double precision, for |r| close to 1, % and for Matlab by Alan Genz. Minor bug modifications 7/98, 2/10. % if dh == inf | dk == inf, p = 0; elseif dh == -inf, if dk == -inf, p = 1; else p = phid(dk); end elseif dk == -inf, p = phid(dh); else if abs(r) < 0.3, ng = 1; lg = 3; % Gauss Legendre points and weights, n = 6 w(1:3,1) = [0.1713244923791705 0.3607615730481384 0.4679139345726904]'; x(1:3,1) = [0.9324695142031522 0.6612093864662647 0.2386191860831970]'; elseif abs(r) < 0.75, ng = 2; lg = 6; % Gauss Legendre points and weights, n = 12 w(1:3,2) = [.04717533638651177 0.1069393259953183 0.1600783285433464]'; w(4:6,2) = [0.2031674267230659 0.2334925365383547 0.2491470458134029]'; x(1:3,2) = [0.9815606342467191 0.9041172563704750 0.7699026741943050]'; x(4:6,2) = [0.5873179542866171 0.3678314989981802 0.1252334085114692]'; else, ng = 3; lg = 10; % Gauss Legendre points and weights, n = 20 w(1:3,3) = [.01761400713915212 .04060142980038694 .06267204833410906]'; w(4:6,3) = [.08327674157670475 0.1019301198172404 0.1181945319615184]'; w(7:9,3) = [0.1316886384491766 0.1420961093183821 0.1491729864726037]'; w(10,3) = 0.1527533871307259; x(1:3,3) = [0.9931285991850949 0.9639719272779138 0.9122344282513259]'; x(4:6,3) = [0.8391169718222188 0.7463319064601508 0.6360536807265150]'; x(7:9,3) = [0.5108670019508271 0.3737060887154196 0.2277858511416451]'; x(10,3) = 0.07652652113349733; end h = dh; k = dk; hk = h*k; bvn = 0; if abs(r) < 0.925, hs = ( h*h + k*k )/2; asr = asin(r); for i = 1 : lg sn = sin( asr*( 1 - x(i,ng) )/2 ); bvn = bvn + w(i,ng)*exp( ( sn*hk - hs )/( 1 - sn*sn ) ); sn = sin( asr*( 1 + x(i,ng) )/2 ); bvn = bvn + w(i,ng)*exp( ( sn*hk - hs )/( 1 - sn*sn ) ); end, bvn = bvn*asr/( 4*pi ); bvn = bvn + phid(-h)*phid(-k); else, twopi = 2*pi; if r < 0, k = -k; hk = -hk; end if abs(r) < 1 as = ( 1 - r )*( 1 + r ); a = sqrt(as); bs = ( h - k )^2; c = ( 4 - hk )/8 ; d = ( 12 - hk )/16; asr = -( bs/as + hk )/2; if asr > -100 bvn = a*exp(asr)*( 1 - c*(bs-as)*(1-d*bs/5)/3 + c*d*as*as/5 ); end if hk > -100, b = sqrt(bs); sp = sqrt(twopi)*phid(-b/a); bvn = bvn - exp(-hk/2)*sp*b*( 1 - c*bs*( 1 - d*bs/5 )/3 ); end, a = a/2; for i = 1 : lg for is = -1 : 2 : 1, xs = ( a + a*is*x(i,ng) )^2; rs = sqrt( 1 - xs ); asr = -( bs/xs + hk )/2; if asr > -100, sp = ( 1 + c*xs*( 1 + d*xs ) ); ep = exp( -hk*( 1 - rs )/( 2*( 1 + rs ) ) )/rs; bvn = bvn + a*w(i,ng)*exp(asr)*( ep - sp ); end end end, bvn = -bvn/twopi; end if r > 0, bvn = bvn + phid( -max( h, k ) ); end if r < 0, bvn = -bvn + max( 0, phid(-h)-phid(-k) ); end end, p = max( 0, min( 1, bvn ) ); end % % end bvnu % function p = phid(z), p = erfc( -z/sqrt(2) )/2; % Normal cdf % % end phid %