function [ po, e, lw, up ] = tribnd( c, ai, bi, t ) % % [ p e lw up ] = tribnd( c, ai, bi, t ) % computes modified Bonferroni and Hunter-Worsley bounds, with % trivariate corrections, for MVN probabilities. % Input is covariance matrix c, and lower and upper % integration limit vectors ai and bi. % An optional third input parameter t is used to select for % t =< 1 : second order bounds only (default) % t > 1 : hybrid second order bounds % t > 2 : third order bounds. % Output is probability p, with error bound e, % and best lower and upper bounds; % Example use: % c = [ 1 .3 .2 .1; .3 1 .2 .2; .2 .2 1 .3; .1 .2 .3 1]; % [p e lw up] = tribnd(c,-inf*ones(4,1),[3 5 2 2]'); disp([lw p up e]) % % Second Order bounds references: % Hsu, Jason C. (1996), Multiple Comparisons, % Chapman and Hall, London, p. 230; % Dawson, D. and Sankoff, A. (1967), An Inequality for Probabilities, % Proc. AMS 18, pp. 504-507. % Hybrid Second-Third Order bounds References: % Bukszar, J. and Prekopa, A. (2001), Probability Bounds with Cherry % Trees, Math. Oper. Res. 26, pp. 174-192; % Tomescu, I. (1986), Hypertrees and Bonferroni Inequalities, % Journal of Combinatorial Theory 41, pp. 209-217. % d = 1./sqrt(diag(c)); a = d.*ai; b = d.*bi; d = diag(d); r = d*c*d; n = length(a); s1 = 0; s2 = 0; s3 = 0; mxd = 0*[1:n]; for i = 1 : n, p(i) = phid(-b(i)) + phid(a(i)); s1 = s1 + p(i); q(i,i) = 0; for j = 1 : i - 1, ab = [a(i) a(j)]; bb = [b(i) b(j)]; q(i,j) = bvnc( r(i,j), ab, bb ); q(j,i) = q(i,j); s2 = s2 + q(i,j); end, q(i,i) = -1; end, blw = 1 - s1; bup = 1 - min(p); % disp([blw bup]); % Bonferroni bounds k = 1 + floor( 2*s2/s1 ); % Modified Bonferroni upper bound upb = min( 1 - 2*( s1 - s2/k )/( k + 1 ), 1 ); up = upb; % % Determine maximum spanning tree % lwb = 1 - s1; ne = ones(1,n); mxd = q(:,1)'; j = 2; l = 2; for k = 1 : n - 1 for i = 2 : n, if mxd(i) > mxd(j), j = i; end, end lwb = lwb + mxd(j); mxd(j) = -1; if q(j,ne(j)) > q(l,ne(l)), l = j; end for i = 2 : n if mxd(i) >= 0 & mxd(i) < q(i,j); mxd(i) = q(i,j); ne(i) = j; end end end, lw = max( 0, lwb ); % Hunter-Worsley bound if nargin > 3 & t > 1 % % Improve Hunter-Worsley lower bound using Cherry Tree. % lwbt = lwb; ne(1) = 0; mxd(1) = -1; mxd(l) = 1; mxd(ne(l)) = 1; for l = 3 : n, w = -1; for k = 1 : n if mxd(k) < 0 for j = 1 : n if mxd(j) > 0 & ( j == ne(k) | k == ne(j) ) for i = 1 : n if mxd(i) > 0 & ( i == ne(j) | j == ne(i) ) if q(i,k) > w, w = q(i,k); im = i; jm = j; km = k; end end end end end end end i = im; j = jm; k = km; mxd(k) = 1; at = [ a(i) a(j) a(k) ]; bt = [ b(i) b(j) b(k)]; lwbt = lwbt + q(i,k) - tvnc( [r(i,j) r(i,k) r(j,k)], at, bt ); end % % Compute Tomescu Bound % upbt = 1 - s1 + s2; for i = 2 : n for j = 1 : i-1 if j ~= ne(i) & i ~= ne(j); q(i,j) = -1; q(j,i) = -1; end end end for l = 3 : n, w = -1; for k = 1 : n, e = 0; for i = 1 : n, if q(i,k) >= 0, e = e + 1; j = i; end, end if e == 1 & q(j,k) > w, w = q(j,k); v = k; end end, k = v; for i = 1 : n, if q(i,k) >= 0, j = i; end, end, q(k,j) = -1; q(j,k) = -1; for i = 2 : n, j = ne(i); if q(i,j) >= 0, at = [ a(i) a(j) a(k) ]; bt = [ b(i) b(j) b(k)]; upbt = upbt - tvnc( [r(i,j) r(i,k) r(j,k)], at, bt ); end end end lw = max( [ 0 lw lwbt] ); up = min( [ upbt up 1 ] ); end if nargin > 3 & t > 2 for i = 1 : n for j = 1 : i - 1 for k = 1 : j - 1, at = [ a(i) a(j) a(k) ]; bt = [ b(i) b(j) b(k)]; s3 = s3 + tvnc( [r(i,j) r(i,k) r(j,k)], at, bt ); end end end j = 2 + floor( 3*s3/s2 ); % Third Order Bounds lwt = 1 - s1 + 2*( (2*j-1)*s2 - 3*s3 )/(j*(j+1)); i = 1 + floor( 2*( (n-2)*s2 - 3*s3 )/( (n-1)*s1 - 2*s2 ) ); upt = 1 - ( (i+2*n-1)*s1-2*((2*i+n-2)*s2-3*s3)/i )/(n*(i+1)); lw = max( [ 0 lw lwt] ); up = min( [ upt up 1 ] ); % disp([lwb lwt lwbt upbt upt upb]) end po = ( lw + up )/2; e = ( up - lw )/2; % function p = bvnc( r, a, b ) % % Computes complementary bvn value for % lower and upper limits a and b % correlation coefficient r % % x [>b(1) >b(2)] x [b(2)] p = bvnl( -b(1), -b(2), r ) + phid(a(1)) - bvnl( a(1), b(2), r ); % x [>b(1) b(1) >b(2) >b(3)] p = tvnl( -b, r ); % x [b(2) >b(3)] p = p + bvnl( -b(2), -b(3), r(3) ) - tvnl( -[a(1) b(2) b(3)], r ); % x [>b(1) b(3)] p = p + bvnl( -b(1), -b(3), r(2) ) - tvnl( -[b(1) a(2) b(3)], r ); % x [b(3)] p = p + bvnl( a(1), a(2), r(1) ) - tvnl( [a(1) a(2) b(3)], r ); % x [>b(1) >b(2) b(2) b(1) abs(r13) ) h2 = h3; h3 = h(2); r12 = r13; r13 = r(1); end if ( abs(r13) > abs(r23) ) h1 = h2; h2 = h(1); r23 = r13; r13 = r(3); end tvn = 0; if ( abs(h1) + abs(h2) + abs(h3) < epst ) tvn = ( 1 + ( asin(r12) + asin(r13) + asin(r23) )/pt )/8; elseif ( abs(r12) + abs(r13) < epst ) tvn = phid(h1)*bvnl( h2, h3, r23 ); elseif ( abs(r13) + abs(r23) < epst ) tvn = phid(h3)*bvnl( h1, h2, r12 ); elseif ( abs(r12) + abs(r23) < epst ) tvn = phid(h2)*bvnl( h1, h3, r13 ); elseif ( 1 - r23 < epst ) tvn = bvnl( h1, min( [h2 h3] ), r12 ); elseif ( r23 + 1 < epst ) if ( h2 > -h3 ) tvn = bvnl( h1, h2, r12 ) - bvnl( h1, -h3, r12 ); end else % % Compute singular TVN value % tvn = bvnl( h2, h3, r23 )*phid(h1); % % Use numerical integration to compute probability; % % TVNH1 = h1; TVNH2 = h2; TVNH3 = h3; TVNR23 = r23; TVNRUA = asin( r12 ); TVNRUB = asin( r13 ); tvn = tvn + adonet( 'tvnmfn', 0, 1, epst )/( 2*pi ); end tvn = max( 0, min( tvn, 1 ) ) ; % % end tvnl % function f = tvnmfn( x ) % % Computes Plackett formula integrands % global TVNH1 TVNH2 TVNH3 TVNR23 TVNRUA TVNRUB f = 0; [ r12 rr2 ] = sincs( TVNRUA*x ); [ r13 rr3 ] = sincs( TVNRUB*x ); if ( abs(TVNRUA) > 0 ) f = f + TVNRUA*pntgnd( TVNH1,TVNH2,TVNH3,r13,TVNR23,r12,rr2 ); end if ( abs(TVNRUB) > 0 ) f = f + TVNRUB*pntgnd( TVNH1,TVNH3,TVNH2,r12,TVNR23,r13,rr3 ); end % % end tvnmfn % function [sx, cs] = sincs( x ) % % Computes sin(x), cos(x)^2, with series approx. for |x| near pi/2 % ee = ( pi/2 - abs(x) )^2; if ( ee < 5e-5 ) sx = ( 1 - ee*( 1 - ee/12 )/2 )*sign( x ); cs = ee*( 1 - ee*( 1 - 2*ee/15 )/3 ); else sx = sin(x); cs = 1 - sx*sx; end % % end sincs % function f = pntgnd( ba, bb, bc, ra, rb, r, rr ) % % Computes Plackett formula integrand; % f = 0; dt = rr*( rr - ( ra - rb )^2 - 2*ra*rb*( 1 - r ) ); if ( dt > 0 ) bt = ( bc*rr + ba*( r*rb - ra ) + bb*( r*ra -rb ) )/sqrt(dt) ; ft = ( ba - r*bb )^2/rr + bb*bb; if ( bt > -10 & ft < 100 ) f = exp( -ft/2 ); if ( bt < 10 ), f = f*phid(bt); end end end % % end pntgnd % function p = phid(z) p = erfc( -z/sqrt(2) )/2; % % end phid % function p = bvnl( dh, dk, r ) % % A function for computing bivariate normal probabilities. % bvnl calculates the probability that x < dh and y < dk. % parameters % dh 1st upper integration limit % dk 2nd upper integration limit % r correlation coefficient % % 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 - last modifications 7/98. % p = bvnu( -dh, -dk, r ); % % end bvnl % function p = bvnu( dh, dk, r ) % % 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 % % 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 - last modifications 7/98. % Note: to compute the probability that x < dh and y < dk, use % bvnu( -dh, -dk, r ). % 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; twopi = 2*pi; if abs(h) > 10, if h < -10, bvn = phid(-k); end elseif abs(k) > 10, if k < -10, bvn = phid(-h); end else 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 ) + phid(-h)*phid(-k) ; else 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*( is*x(i,ng) + 1 ) )^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 end p = max( 0, min( 1, bvn ) ); % % end bvnu % function fin = adonet( f, a, b, tol ) % % 1-dimensional adaptive integration % nl = 100; ai(1) = a; bi(1) = b; ip = 1; im = 1; err = 1; while ( 4*err > tol & im < nl ) im = im + 1; bi(im) = bi(ip); ai(im) = ( ai(ip) + bi(ip) )/2; bi(ip) = ai(im); [ fi(ip) ei(ip) ] = krnrdt( ai(ip), bi(ip), f ); [ fi(im) ei(im) ] = krnrdt( ai(im), bi(im), f ); fin = sum( fi(1:im) ); err = sqrt( sum( ei(1:im).^2 ) ); for i = 1 : im, if ei(i) > ei(ip), ip = i; end, end end % function [ resk, err ] = krnrdt( a, b, f ) % % kronrod rule % wg0 = 0.2729250867779007e+00; wg(1:3) = [ .05566856711617449 0.1255803694649048 0.1862902109277352 ]; wg(4:5) = [ 0.2331937645919914 0.2628045445102478 ]; % xgk(1:3) = [ 0.9963696138895427 0.9782286581460570 0.9416771085780681 ]; xgk(4:6) = [ 0.8870625997680953 0.8160574566562211 0.7301520055740492 ]; xgk(7:9) = [ 0.6305995201619651 0.5190961292068118 0.3979441409523776 ]; xgk(10:11) = [ 0.2695431559523450 0.1361130007993617 ]; % wgk0 = 0.1365777947111183e+00; wgk(1:3) = [ .00976544104596129 .02715655468210443 .04582937856442671 ]; wgk(4:6) = [ .06309742475037484 .07866457193222764 .09295309859690074 ]; wgk(7:9) = [ 0.1058720744813894 0.1167395024610472 0.1251587991003195 ]; wgk(10:11) = [ 0.1312806842298057 0.1351935727998845 ]; wid = ( b - a )/2; cen = ( b + a )/2; fc = feval( f, cen ); resg = fc*wg0; resk = fc*wgk0; for j = 1 : 5 t = wid*xgk(2*j-1); fc = feval( f, cen - t ) + feval( f, cen + t ); resk = resk + wgk(2*j-1)*fc; t = wid*xgk(2*j); fc = feval( f, cen - t ) + feval( f, cen + t ); resk = resk + wgk(2*j)*fc; resg = resg + wg(j)*fc; end t = wid*xgk(11); fc = feval( f, cen - t ) + feval( f, cen + t ); resk = wid*( resk + wgk(11)*fc ); err = abs( resk - wid*resg ); % % end krnrdt %