function [ p, e ] = qsclatmvnv( m, r, a, cn, b ) % % [ P E ] = QSCLATMVNV( M, R, A, CN, B ) % uses a randomized quasi-random rule with m points to estimate an % MVN probability for positive semi-definite covariance matrix r, % with constraints a < cn*x < b. If r is nxn and cn is kxn, then % a and b must be column k-vectors. % Probability p is output with error estimate e. % Example use: % >> r = [ 4 3 2 1; 3 5 -1 1; 2 -1 4 2; 1 1 2 5 ]; % >> a = [ -inf 1 -5 ]'; b = [ 3 inf 4 ]'; % >> cn = [ 1 2 3 -2; 2 4 1 2; -2 3 4 1 ]; % >> [ p e ] = qsclatmvnv( 5000, r, a, cn, b ); disp([ p e ]) % % This function uses an algorithm given in the paper by Alan Genz: % "Numerical Computation of Multivariate Normal Probabilities", in % J. of Computational and Graphical Stat., 1(1992), 141-149. % The primary references for the numerical integration are % "Randomization of Number Theoretic Methods for Multiple Integration" % R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13(1976), pp. 904-14. % and % "Fast Component-by-Component Construction, a Reprise for Different % Kernels", D. Nuyens and R. Cools. In H. Niederreiter and D. Talay, % editors, Monte-Carlo and Quasi-Monte Carlo Methods 2004, % Springer-Verlag, 2006, 371-385. % % Alan Genz is the author of this function and following Matlab functions. % Alan Genz, WSU Math, PO Box 643113, Pullman, WA 99164-3113 % Email : alangenz@wsu.edu % % Initialization % [ as ch bs clg n ] = chlsrt( r, a, cn, b ); ns = 10; [ q pr ] = fstrnk( m/ns, n-1 ); q = q/pr; ai = min( max( as(1), -9 ) , 9 ); c = phi(ai); bi = min( max( ai, max( bs(1), -9 ) ), 9 ); d = phi(bi); ci = c; dci = d - ci; y = zeros(n-1,pr); p = 0; e = 0; % % Randomization loop for ns samples % for S = 1 : ns, c = ci; dc = dci; vp = dc; li = 2; lf = 1; % % Compute randomized MVN integrand with tent periodization % for i = 2 : n, lf = lf + clg(i); y(i-1,:) = phinv( c + abs( 2*mod(q(i-1)*[1:pr]+rand,1) - 1 ).*dc ); if lf < li, c = 0; dc = 1; else, ai = -9; bi = 9; for j = li : lf, s = ch(j,1:i-1)*y(1:i-1,:); ai = max( as(j) - s, ai ); bi = min( bs(j) - s, bi ); end, ai = min( ai, 9 ); bi = max( ai, max( bi, -9 ) ); c = phi(ai); dc = phi(bi) - c; vp = vp.*dc; end, li = li + clg(i); end, d = ( mean(vp) - p )/S; p = p + d; if abs(d) > 0, e = abs(d)*sqrt( 1 + ( e/d )^2*( S - 2 )/S ); else, if S > 1, e = e*sqrt( ( S - 2 )/S ); end end end, e = 3*e; % error estimate is 3 x standard error with ns samples. % % end qsclatmvnv % function [ ap, ch, bp, clg, np ] = chlsrt( r, a, cn, b ) % % Computes permuted lower Cholesky factor ch for covariance r which % may be singular, combined with contraints a < cn*x < b, to % form revised lower triangular constraint set ap < ch*x < bp; % clg contains information about structure of ch: clg(1) rows for % ch with 1 nonzero, ..., clg(np) rows with np nonzeros. % ep = 1e-10; % singularity tolerance; % [n n] = size(r); c = r; [m n] = size(cn); ch = cn; dc = sqrt(max(diag(c),0)); for i = 1 : n, d = dc(i); if d > 0, c(:,i) = c(:,i)/d; c(i,:) = c(i,:)/d; ch(:,i) = ch(:,i)*d; end end, ap = a; bp = b; np = 0; % % determine (with pivoting) Cholesky factor for r % and form revised constraint matrix ch % for i = 1 : n, clg(i) = 0; epi = ep*i; j = i; for l = i+1 : n, if c(l,l) > c(j,j), j = l; end, end if j > i, t = c(i,i); c(i,i) = c(j,j); c(j,j) = t; t = c(i,1:i-1); c(i,1:i-1) = c(j,1:i-1); c(j,1:i-1) = t; t = c(i+1:j-1,i); c(i+1:j-1,i) = c(j,i+1:j-1)'; c(j,i+1:j-1) = t'; t = c(j+1:n,i); c(j+1:n,i) = c(j+1:n,j); c(j+1:n,j) = t; t = ch(:,i); ch(:,i) = ch(:,j); ch(:,j) = t; end if c(i,i) < epi, break, end, cvd = sqrt( c(i,i) ); c(i,i) = cvd; for l = i+1 : n c(l,i) = c(l,i)/cvd; c(l,i+1:l) = c(l,i+1:l) - c(l,i)*c(i+1:l,i)'; end, ch(:,i) = ch(:,i:n)*c(i:n,i); np = np + 1; end, y = zeros(n,1); sqtp = sqrt(2*pi); % % use right reflectors to reduce ch to lower triangular % for i = 1 : min( np-1, m ), epi = ep*i; vm = 1 + epi; lm = i; % % permute rows so that smallest variance variables are first. % for l = i : m, v = ch(l,1:np); s = v(1:i-1)*y(1:i-1); ss = max( sqrt( sum( v(i:np).^2 ) ), epi ); al = ( ap(l) - s )/ss; bl = ( bp(l) - s )/ss; dna = 0; dsa = 0; dnb = 0; dsb = 1; if al > -9, dna = exp(-al*al/2)/sqtp; dsa = phi(al); end if bl < 9, dnb = exp(-bl*bl/2)/sqtp; dsb = phi(bl); end, p = dsb - dsa; if p > epi, mn = dna - dnb; vr = al*dna - bl*dnb; if al <= -9, mn = -dnb; vr = -bl*dnb; elseif bl >= 9, mn = dna; vr = al*dna; end, mn = mn/p; vr = 1 + vr/p - mn^2; else, vr = 0; mn = ( al + bl )/2; if al < -9, mn = bl; elseif bl > 9, mn = al; end end, if vr <= vm, lm = l; vm = vr; y(i) = mn; end end, v = ch(lm,1:np); if lm > i, ch(lm,1:np) = ch(i,1:np); ch(i,1:np) = v; tl = ap(i); ap(i) = ap(lm); ap(lm) = tl; tl = bp(i); bp(i) = bp(lm); bp(lm) = tl; end, ch(i,i+1:np) = 0; ss = sum( v(i+1:np).^2 ); if ss > epi, ss = sqrt( ss + v(i)^2 ); if v(i) < 0, ss = -ss; end ch(i,i) = -ss; v(i) = v(i) + ss; vt = v(i:np)'/( ss*v(i) ); ch(i+1:m,i:np) = ch(i+1:m,i:np) - ch(i+1:m,i:np)*vt*v(i:np); end end % % scale and sort constraints % for i = 1 : m, v = ch(i,1:np); clm(i) = min(i,np); jm = 1; for j = 1 : clm(i), if abs(v(j)) > ep*j, jm = j; end, end if jm < np, v(jm+1:np) = 0; end, clg(jm) = clg(jm) + 1; at = ap(i); bt = bp(i); j = i; for l = i-1 : -1 : 1 if jm >= clm(l), break, end ch(l+1,1:np) = ch(l,1:np); j = l; ap(l+1) = ap(l); bp(l+1) = bp(l); clm(l+1) = clm(l); end, clm(j) = jm; vjm = v(jm); ch(j,1:np) = v/vjm; ap(j) = at/vjm; bp(j) = bt/vjm; if vjm < 0, tl = ap(j); ap(j) = bp(j); bp(j) = tl; end end, j = 0; for i = 1 : np, if clg(i) > 0, j = i; end, end, np = j; % % combine constraints for first variable % if clg(1) > 1 ap(1) = max( ap(1:clg(1)) ); bp(1) = max( ap(1), min( bp(1:clg(1)) ) ); ap(2:m-clg(1)+1) = ap(clg(1)+1:m); bp(2:m-clg(1)+1) = bp(clg(1)+1:m); ch(2:m-clg(1)+1,:) = ch(clg(1)+1:m,:); clg(1) = 1; end % % end chlsrt % function p = phi(z) % % Standard statistical normal distribution cdf % p = erfc( -z/sqrt(2) )/2; return % % end phi % function z = phinv(w) % % Standard statistical inverse normal distribution % %z = -sqrt(2)*erfcinv( 2*w ); z = norminv(w); return % % end phinv % function [ z, n ] = fstrnk( ni, sm, om, gm, bt ) % % Reference: % "Fast Component-by-Component Construction, a Reprise for Different % Kernels", D. Nuyens and R. Cools. In H. Niederreiter and D. Talay, % editors, Monte-Carlo and Quasi-Monte Carlo Methods 2004, % Springer-Verlag, 2006, 371-385. % Modifications to original by A. Genz, 05/07 % Typical Use: % om = inline('x.^2-x+1/6'); n = 99991; s = 100; gam = 0.9.^[1:s]; % z = fastrank( n, s, om, gam, 1 + gam/3 ); disp([z'; e]) % n = fix(ni); if ~isprime(n), pt = primes(n); n = pt(length(pt)); end if nargin < 3, om = inline('x.^2-x+1/6'); bt = ones(1,sm); gm = [ 1 (4/5).^[0:sm-2] ]; end, q = 1; w = 1; z = [1:sm]'; m = ( n - 1 )/2; g = prmrot(n); perm = [1:m]; for j = 1 : m-1, perm(j+1) = mod( g*perm(j), n ); end perm = min( n - perm, perm ); c = om(perm/n); fc = fft(c); for s = 2 : sm, q = q.*( bt(s-1) + gm(s-1)*c([w:-1:1 m:-1:w+1]) ); [ es w ] = min( real( ifft( fc.*fft(q) ) ) ); z(s) = perm(w); end % % end fstrnk % function r = prmrot(pin) % % find primitive root for prime p, might fail for large primes (p > 32e7) % p = pin; if ~isprime(p), pt = primes(p); p = pt(length(pt)); end pm = p - 1; fp = unique(factor(pm)); n = length(fp); r = 2; k = 1; while k <= n; d = pm/fp(k); rd = r; for i = 2 : d, rd = mod( rd*r, p ); end % computes r^d mod p if rd == 1, r = r + 1; k = 0; end, k = k + 1; end % % prmrot %