function [ pl, pu, el, eu, en, eq ] = MGP( S, u, n, r, Tf, ua, ns, Ti ) % % Computes approximate upper and lower bounds for the % maximum of a stationary Gaussian process P( M_T > u ). % Inputs: % S = sample size for MCQMC integration method; % u = bound level; % n = number of discretization points for upper bound, default 16; % r = r(t) the covariance function for the process, % vectorized, with default exp(-t.^/2); % Tf, Ti : t limits for the process, default Ti = 0, Tf = 1. % ua (default ua = 0 ), if ua ~= 0, then P( |M_T| > u ) is computed; % ns (default ns = 0 ), if ns ~= 0, then r is scaled so that r''(0) = -1; % Outputs: % pl = approximate lower bound, with % el = error estimate (MCQMC); % pu = approximate upper bound, with % eu = estimate for total error; % en = estimate for discretization error, and % eq = estimate for MCQMC error; % % Example: % [ pl pu el eu ] = MGP( 10000, 1, 12, @(t)exp(-t.^2/2) ); % disp([ pl pu; el eu ]) % if nargin < 3, n = 16; end if nargin < 4, r = inline('exp(-t.^2/2)','t'); end if nargin < 5, Tf = 1; end, if nargin < 6, ua = 0; end if nargin < 7, ns = 0; end, if nargin < 8, Ti = 0; end [ sg, cv, cvl, rs, rps, nc, ncl ] = chlmxr( n, r, Ti, Tf, ns ); urs = u/rs; [ pv, ev ] = qsltmvnevsv( S, n, cv, nc, urs, ua ); h = rps*( Tf - Ti )/n; pz = phi(urs)*phi(0); po = Phi(-urs); pe = po + h*( pz + 2*sum(pv(2:2:n-2)) + pv(n) ); pf = po + h*( pz/2 + sum(pv(1:n-1)) + pv(n)/2 ); pu = ( 4*pf - pe )/3; en = abs( pf - pe )/3; eq = h*( sum(ev(1:n-1)) + ev(n)/2 ); eu = eq + en; if ua ~= 0, pu = 2*pu; en = 2*en; eu = 2*eu; end % Lower Bound [ pl, el ] = lwltmvnsv( S, cvl, ncl, urs, ua ); pl = 1 - pl; % % end MGP % function [ sg, cv, cvl, rs, rps, nc, ncl ] = chlmxr( n, rf, Ti, Tf, ns ) % % Computes Cholesky factors for covariance matrix % using covariance function rf % if nargin < 4, Ti = 0; Tf = 1; end, if nargin < 5, ns = 0; end if ns ~= 0, h = 1e-5; % determine scaling constants sqrt(-rf''(0)/rf(0)) rz = rf(0); rs = sqrt(rz); rps = sqrt( (2*rz-rf(h)-rf(-h))/(rz*h^2) ); else, rz = 1; rps = 1; rs = 1; % assume correct scaling end, % determine rf/rf(0) values at scaled times dl = ( Tf - Ti )*[1:n]/n; r = rf(dl)/rz; % determine rf'/rf(0) values at scaled times h = 1e-6; rp = ( rf(dl+h/2) - rf(dl-h/2) )/( rps*rz*h ); sg = zeros(n+2); sg(1,2:n+2) = [0 r]; sg(2,3:n+2) = rp; for i = 3 : n+1, sg(i,i+1:n+2) = r(1:n-i+2); end sg = eye(n+2) + sg + sg'; [ cv , nc ] = chlsvd( sg ); [ cvl, ncl ] = chlsvd( sg(3:n+2,3:n+2) ); % % end chlmax % function [ ch, nr, e ] = chlsvd( sg ), ep = 5e-9; % % Computes lower Cholesky factor for sg, which may be singular % using svd and qr, returning factor ch, and rank nr. % [n n] = size(sg); [ U, D ] = svd(sg); d = sqrt(diag(D)); nr = n; if d(n) < ep, nr = min( find( d < ep ) ); end [ Q, R ] = qr( diag(d(1:nr))*U(:,1:nr)' ); ch = ( diag(sign(diag(R)))*R )'; e = norm( ch*ch' - sg ); % % end chlsvd % function [ pv, ev ] = qsltmvnevsv( m, n, ch, nc, u, ua ) % % [ P E ] = qsltmvnevsv( m, n, ch, nc, u ) % uses a randomized QMC rule with m points to estimate an % MVN expected value, given n x nc Cholesky factor ch of % +semidefinite nxn covariance matrix. % Probability p is output with error estimate e. % Example use: % >> [ p e ] = qsltmvnevsv( 5000, n, ch, nc, u); disp([ p e ]) % % This function uses a modification of an algorithm in the paper % "Numerical Computation of Multivariate Normal Probabilities", % A. Genz, J. of Comp. Graph. Stat., 1(1992), pp. 141-149. % % The primary references for the numerical integration are % "Randomization of Number Theoretic Methods for Multiple Integration", % R. Cranley & 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 & 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 % ns = 12; pr = fix( m/ns ); if nargin < 6, ua = 0; end if nc > 64, q = sqrt(primes(5*nc*log(nc)/4))'; % Richtmyer generators else, [ q pr ] = fstrnk( pr, nc-2 ); q = q/pr; % Lattice rule end, y = zeros(nc-2,pr); uv = u*ones(1,pr); pu = phi(u)/2; pv = 0; ev = 0; % % Randomization loop for ns samples % for S = 1 : ns, c = 0; % % Compute randomized MVN integrand with tent periodization % % y(1,:) = Phinv( ( rand(1,pr) + 1 )/2 ); % MC y(1,:) = Phinv( ( abs(2*mod(q(1)*[1:pr]+rand,1) - 1 ) + 1 )/2 ); % MCQMC fi = pu*y(1,:); for i = 2 : nc-2, ys = ch(i+1,1:i)*[ uv; y(1:i-1,:) ]; if ua ~=0, c = Phi( max( -( u + ys )/ch(i+1,i+1), -9 ) ); end dc = Phi( min( ( u - ys )/ch(i+1,i+1), 9 ) ) - c; fi = fi.*dc; df(i-1) = mean(fi); % y(i,:) = Phinv( c + rand(1,pr).*dc ); % MC y(i,:) = Phinv( c + abs( 2*mod(q(i)*[1:pr]+rand,1) - 1 ).*dc ); % MCQMC end, bi = 9; ai = -9; for i = nc : n+2, ys = ch(i,1:nc-1)*[ uv; y(1:nc-2,:) ]; if ch(i,nc) > 0 if ua ~= 0, ai = max( ( -uv - ys )/ch(i,nc), ai ); end bi = min( ( uv - ys )/ch(i,nc), bi ); else ai = max( ( uv - ys )/ch(i,nc), ai ); if ua ~= 0, bi = min( ( -uv - ys )/ch(i,nc), bi ); end end end, df(nc-2) = mean( fi.* ( Phi(bi) - Phi( min(ai,bi) ) ) ); df = ( df - pv )/S; pv = pv + df; ev = ( S - 2 )*ev/S + df.^2; end, ev = 3*sqrt(ev); % error estimate is 3 x std error with ns samples. if nc < n + 2, pv(nc-1:n) = pv(nc-2); ev(nc-1:n) = ev(nc-2); end % % end qsltmvnvevsv % function [ p, e ] = lwltmvnsv( m, ch, rc, u, ua ) % % [ p, e ] = lwltmvnsv( M, ch, rc, u ) % uses a randomized QMC rule with m points to estimate an % MVN probability for Cholesky factor of +semidef covariance matrix, % with rank rc, lower limits -inf (or -u if ua ~= 0), upper limits u. % Probability p is output with error estimate e. % Example use: % >> ch = [1 0 0 0;.1 .5 0 0;.1 .2 .4 0;.3 0 0 0]; % >> [ p e ] = lwltmvnsv( 5000, ch, rc, u ); disp([ p e ]) % % This function uses a modification of an algorithm in the paper % "Numerical Computation of Multivariate Normal Probabilities", % A. Genz, J. of Comp. Graph. Stat., 1(1992), pp. 141-149. % % The primary references for the numerical integration are % "Randomization of Number Theoretic Methods for Multiple Integration", % R. Cranley & 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 & 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 % ns = 12; pr = fix(m/ns); [n,n] = size(ch); ci = 0; % disp(ch) if ua ~= 0, ci = Phi(-u/ch(1,1)); end, dci = Phi(u/ch(1,1)) - ci; if rc > 64, q = sqrt(primes(5*rc*log(rc+1)/4))'; % Richtmyer generators else, [ q pr ] = fstrnk( pr, rc-1 ); q = q/pr; % Lattice rule end, y = zeros(rc-1,pr); on = ones(1,pr); uv = u*on; p = 0; e = 0; % % Randomization loop for ns samples % for S = 1 : ns, c = ci; dc = dci; vp = dc; bi = 9; ai = -9; % % Compute randomized MVN integrand with tent periodization % for i = 2 : rc-1 y(i-1,:) = Phinv( c + abs( 2*mod(q(i-1)*[1:pr]+rand,1) - 1 ).*dc ); ys = ch(i,1:i-1)*y(1:i-1,:); if ua ~= 0; c = Phi( max( -( uv + ys )/ch(i,i), -9 )); end dc = Phi( min( ( uv - ys )/ch(i,i), 9 ) ) - c; vp = vp.*dc; end, y(rc-1,:) = Phinv( c + abs( 2*mod(q(rc-1)*[1:pr]+rand,1) - 1 ).*dc ); for i = rc : n, ys = ch(i,1:rc-1)*y(1:rc-1,:); if ch(i,rc) > 0 if ua ~= 0, ai = max( -( uv + ys )/ch(i,rc), ai ); end bi = min( ( uv - ys )/ch(i,rc), bi ); else ai = max( ( uv - ys )/ch(i,rc), ai ); if ua ~= 0, bi = min( -( uv + ys )/ch(i,rc), bi ); end end end, df = ( mean( vp.*( Phi(bi) - Phi(min(ai,bi)) ) ) - p )/S; p = p + df; e = ( S - 2 )*e/S + df^2; end, e = 3*sqrt(e); % error estimate is 3 x standard error with ns samples. % % end lwltmvnsv % function p = phi(x), p = exp(-x.^2/2)/sqrt(2*pi); % Normal pdf % end phi function p = Phi(z), p = erfc(-z/sqrt(2))/2; % Normal cdf % end Phi function z = Phinv(w), z = -sqrt(2)*erfcinv(2*w); % Inverse Normal cdf % 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 = fstrnk( 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