/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ This code is part of the paper "Simulation of Multivariate Normal Orthant Probabilities: Methods and Programs" by V.A. Hajivassiliou, Yale University, D.L. McFadden, University of California at Berkeley, and P.A. Ruud, University of California at Berkeley July 1992 A. Overview This section contains a series of GAUSS procedures for simulation of multivariate normal rectangle probabilities, and the derivatives of these probabilities with respect to the mean and covariances of the normal distri- bution. Table 2 indexes the procedures in this section. TABLE 2. GAUSS PROCEDURES SECTION MNEUMONIC* DESCRIPTION I CFS Crude Frequency Simulator J NIS Normal Importance Sampling Simulator Procs NISE and NIST use exponential and truncated normal comparison distributions, respectively K KFS Polynomial Kernel Smoothed Frequency Simulator The proc is PKFS L SDS Stern Decomposition Simulator M GHK Geweke-Hajivassiliou-Keane Simulator N PCF Parabolic Cylinder Function Simulator O DCS Deak chi-square Simulator P ARS Acceptance/Rejection Simulator Procs ARSE and ARSR use exponential and recursive truncated normal comparison distributions, respectively Q GSS Gibbs Sampler Simulator R AUS Approximate (OR Sequential) Unbiased Simulator S SUS Sequential Unbiased Simulator T LIB_A Statistical Functions U LIB_B Spherical Transformations and Antithetics V LIB_C Miscellaneous Utilities * The procs are located in libraries named *.G, where '*' denotes the mneumonic. Each of the simulator procedures requires the following standard inputs; the interpretation of some inputs may vary from routine to routine, and not all are used in some routines: M Dimension of the multivariate normal MU Mean of multivariate normal, a M x 1 vector W Covariance matrix of multivariate normal, a M x M array WI Inverse of W C Lower triangular Cholesky factor of W, a M x M array A Lower bound of rectangle, a M x 1 vector [When the lower bound is -infinity, set A = (-1.0E10)*ONES(M,1)] B Upper bound of rectangle, a M x 1 vector [When the upper bound is +infinity, set B = (+1.0E10)*ONES(M,1)] R Number of repetitions U Random variates, a M x R array PARM Parameters and constants for the simulation The simulators all return {P,HU,HC}, where P is the scalar rectangle probabi- lity, HU is the (M+1) x (M+1) array of unconditional partial moments (6), and HC is the (M+1) x (M+1) array of conditional moments (7). Parts of the out- put not provided by a simulator are set to -999. Interpretations of input parameters ----------------------------------- ARSR,ARSE --------- These procs assume that PARM[.,1] contains seeds for the uniform random number generator RNDUS called by the routines; there is one seed for each repetition R. PARM[R+1,1] is assumed to contain a censoring limit for rejection. ARSR and ARSE use U = UUBIG; (uniforms) AUS --- The proc assumes that PARM(.,1) contains a vector of R seeds for the normal random number generator RNDNS, PARM(1,2) contains the number of initial GHK simulators of P used in the simulation of 1/P, PARM(2,2) contains the maximum number of trials using the crude frequency simulator that are made before the process is censored. If PARM(2,2) is made large (say 1.e50), then this is computationally the same as the Sequential Unbiased Simulator (SUS). PARM(3,2) contains a seed for RNDUS. AUS uses U = UUBIG[.,1:R]; (uniforms) CFS --- This proc does not use PARM. CFS uses U = UNBIG[.,1:R] (standard normals) DCS --- This proc assumes that PARM(1,1) = DET(W). DCS uses U = USBIG[.,1:R]; (random normal grid on unit sphere) GHK --- This proc does not use PARM. GHK uses U = UUBIG[.,1:R]; (uniforms) GSS --- It assumes that PARM[1,1] gives the number of rounds to be used in the construction of the random draws, JL. GSS uses U = UUBIG[.,1:JL]; (uniforms) NISE,NIST --------- These procs do not use PARM. NISE and NIST use U = UUBIG[.,1:R]; (uniforms) PCF --- This proc assumes that PARM(1,1) = DET(W). PCF uses U = USSBIG[.,1:R]; (random normal grid on unit sphere and orthant) PKF --- This proc does not use W. It assumes that PARM[1,1] contains a window width parameter. KFS uses U = UNBIG[.,1:R]; (standard normals) SDS --- The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such that W-lI is positive definite. SDS uses U = UNBIG[.,1:R]; (standard normals) Non_standard feature: -------------------- SDS assumes that C is a Cholesky factor of W-lI, not of W. SUS --- This procedure is the same the Asymptotically Unbiased Simulator (AUS) with PARM(2,2) made large (1.e50). The proc assumes that PARM(.,1) contains a vector of R seeds for the normal random number generator RNDNS, PARM(1,2) contains the number of initial GHK simulators of P used in the simulation of 1/P, PARM(2,2) contains the maximum number of trials using the crude frequency simulator that are made before the process is censored. This is made very large. PARM(3,2) contains a seed for RNDUS. SUS uses U = UUBIG[.,1:R]; (uniforms) ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Test Harness The following program provides a simple harness for testing the procs above. A one-factor model is used to generate probabilities, with gaussian quadrature used for evaluation. The accuracy of this quadrature can be checked against symmetric cases where the probabilities are known. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ NEW; BEGIN: "INPUT DIMENSION M";; M = CON(1,1); "SEED FOR RANDOM NUMBER GENERATOR";; SEED = CON(1,1); SEED1 = SEED; "NUMBER OF REPETITIONS FOR NIS, GHK, SDS (CFS, KFS USE 10X)";; RBASE = CON(1,1); RBIG = 30*RBASE; "DO YOU WANT A SPECIAL CASE WITH P ANALYTIC? ENTER 0 = NO, 1 = FACTOR"; " LOADING ON M ONLY, 2 = ORTHANT PROB. WITH FULL SYMMETRY ";; ICASE = CON(1,1); IF ICASE == 2; MU = ZEROS(M,1); B = ZEROS(M,1); A = -25*ONES(M,1); L = ONES(M,1); D = ONES(M,1); ELSE; "DO YOU WANT AUTOMATIC GENERATION OF DATA? ENTER 1 = YES, 0 = NO ";; DDD=CON(1,1); IF DDD == 1; MU = RNDNS(M,1,SEED); B = ZEROS(M,1); A = -2 - ABS(RNDNS(M,1,SEED)); IF ICASE == 0; L = .3*RNDNS(M,1,SEED); ELSE; L = ZEROS(M-1,1)|(.3*RNDNS(1,1,SEED)); ENDIF; D = .5+.3*ABS(RNDNS(M,1,SEED)); ELSE; "ENTER MEANS"; MU = CON(M,1); "ENTER LOWER BOUNDS A"; A = CON(M,1); "ENTER UPPER BOUNDS B"; B = CON(M,1); "ENTER INDEP. ST. DEV."; D = CON(M,1); B = A+1+(B-1-A).*(B.>A); D = 1+(D-1).*(D.>0); IF ICASE == 0; "ENTER FACTOR LOADINGS"; L = CON(M,1); ELSE; "ENTER FACTOR LOADING COMPONENT M";; L = ZEROS(M-1,1)|CON(1,1); ENDIF; ENDIF; ENDIF; OUTPUT FILE=MSM.OUT ON; FORMAT /M1 /RO 12,5; "THE INPUTS FOR THIS EXERCISE ARE"; " MEANS A B SD FACTOR LOADINGS"; MU~A~B~D~L; FORMAT /M1 /RD 12,0; " RBASE SEED"; RBASE~SEED1; " "; FORMAT /M1 /RO 12,5; S = SQRT(D.*D+L.*L); P = CDFN((B-MU)./S)-CDFN((A-MU)./S); I = SORTIND(P); D = D[I]; L = L[I]; A = A[I]; B = B[I]; MU = MU[I]; "AFTER SORTING TO PUT LOW PROB. IN FIRST COMPONENT"; " MEANS A B SD FACTOR LOADINGS"; MU~A~B~D~L; " "; W = VEC(D.*D); W = DIAGRV(EYE(ROWS(W)),W); W = W+(L*(L')); C = CHOL(W)'; WI = INVPD(W); DETW = DET(W); LET ETW[1,4] = {360000 6000 100 1}; LET QW[20,2] = { 0.048613160848618 0.038752973079681, 0.146006911993027 0.038519907742739, 0.243907555937767 0.038055181503296, 0.342668563127518 0.037361584603786, 0.442667961120606 0.036443289369345, 0.544320702552795 0.035305824130774, 0.648093938827515 0.033956021070480, 0.754526972770691 0.032402005046606, 0.864258110523224 0.030653120949864, 0.978062629699707 0.028719885274768, 1.096907615661621 0.026613922789693, 1.222035169601440 0.024347903206944, 1.355093598365784 0.021935453638434, 1.498357176780701 0.019391084089875, 1.655115008354187 0.016730098053813, 1.830419778823853 0.013968503102660, 2.032696485519409 0.011122924275696, 2.277773618698120 0.008210529573262, 2.601793527603149 0.005249142181128, 3.127621173858643 0.002260638633743}; QW = QW|((-QW[.,1])~QW[.,2]); PP = -999; DM = -999; LM = -999; DL = -999; LL = -999; LET EL = {1,1}; LET REP = {1,1}; IF ICASE == 1; S = SQRT(DIAG(W)); T = CDFN((B-MU)./S)-CDFN((A-MU)./S); PP = PRODC(T); BM = (B[M]-MU[M])/S[M]; AM = (A[M]-MU[M])/S[M]; Q = PDFN(BM)-PDFN(AM); QQ = BM*PDFN(BM)-AM*PDFN(AM); DL = -PRODC(T[1:M-1]); DM = DL*Q/S[M]; LM = DM/PP; @ DL = DL*QQ*L[M]/(S[M]*S[M]);@ DL = DL*QQ*L[M]/(S[M]^3); @ vah: ^3 to conform to eq 72 @ LL = DL/PP; ENDIF; IF ICASE == 2; PP = 1/(M+1); DM = -SUMC((CDFN(QW[.,1])^(M-1)).*PDFN(QW[.,1]).*QW[.,2]); LM = DM/PP; DL = SUMC((CDFN(QW[.,1])^(M-1)).*PDFN(QW[.,1]).*QW[.,1].*QW[.,2]); LL = DL/PP; ENDIF; IF M < 3; {P,T,TT} = EXAF(M,A,B,MU,W,L); PP = PP|P; DM = DM|T; DL = DL|TT; LM = LM|(T/P); LL = LL|(TT/P); ELSE; PP = PP|-999; DM = DM|-999; DL = DL|-999; LM = LM|-999; LL = LL|-999; ENDIF; @vah: set to 1 hour limit per method@ hseclim=360000; UUBIG = RNDUS(M,RBIG,SEED); @SET UP UNIFORM RANDOM NUMBERS@ {US,T} = RNDBASE(M,SEED); @RANDOM BASIS@ T = ROUND((RBASE-M-1)/(M-1)); US = ANTITH(M,T,US); @RANDOM GRID IN UNIT SPHERE@ USS = -ABS(US); @SPHERE & ORTHANT@ @ vah:instead of inverting UUBIG to get UNBIG, a new set of normals is drawn @ /*UNBIG = CDFNI(UUBIG); @SET UP NORMAL RANDOM NUMBERS@*/ UNBIG = rndns(M,RBIG,SEED); @SET UP NORMAL RANDOM NUMBERS@ USBIG = UNBIG./(ONES(M,1)*(SQRT(SUMC(UNBIG.*UNBIG))')); USSBIG = -ABS(USBIG); USSBIG = USSBIG-.0000000001*(USSBIG .== 0); E1 = ETW*TIME; {P,T,TT} = FFACT(QW[.,1],A,B,D,L); PP = PP|(P')*QW[.,2]; DM = DM|(T')*QW[.,2]; DL = DL|(TT')*QW[.,2]; LM = LM|DM[3,1]/PP[3,1]; LL = LL|DL[3,1]/PP[3,1]; @ vah: Case 1, M>2: error in quad1, because L[M]=0 after reordering @ @ temporary fix follows @ if dl[1,1] eq 0; dl[1,1]=dl[3,1]; endif; if ll[1,1] eq 0; ll[1,1]=ll[3,1]; endif; "EXACT 1 AND 2 (-999 IF NOT AVAILABLE) FOLLOWED BY QUADRATURE VALUES:"; " P dP/dMU(M) dLN(P)/dMU(M) dP/dL(M) dLN(P)/dL(M)"; FORMAT /M1 /RO 13,6; PP~DM~LM~DL~LL; REP = REP|1; E2 = ETW*TIME; EC = E2-E1; EL = EL|EC; E1 = E2; @CFS@ PARM = 1; @ vah: irrelevant @ {P,HU,HC} = CFS(M,MU,W,WI,C,A,B,RBIG,UNBIG,PARM); EC = SIM(M,MU,W,WI,C,A,B,RBIG,UNBIG,PARM); "CFS: P AND TIME ";; P~EC; @NIS@ PARM = 1; @ vah: irrelevant @ R = 10*RBASE; U = UUBIG[.,1:R]; {P,HU,HC} = NISE(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "NISE: P AND TIME ";; P~EC; @NIST@ PARM = 1; @ vah: irrelevant @ R = 5*RBASE; U = UUBIG[.,1:R]; {P,HU,HC} = NIST(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "NIST: P AND TIME ";; P~EC; @KFS@ PARM = .1*MINC((B-A)|SQRT(D.*D+L.*L)); R = 10*RBASE; U = UNBIG[.,1:R]; {P,HU,HC} = PKFS(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "KFS: P, TIME, WINDOW ";; P~EC~PARM; @SDS@ VA = EIGRS(W); PARM = .99*VA[1,1]; @ vah: parm raised to 2 @ CC = CHOL(W-PARM^2*EYE(M))'; R = RBASE; U = UNBIG[.,1:R]; {P,HU,HC} = SDS(M,MU,W,WI,CC,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,CC,A,B,R,U,PARM); "SDS: P, TIME, LAMBDA ";; P~EC~PARM; @GHK@ PARM = 1; @ vah: irrelevant @ R = 3*RBASE; U = UUBIG[.,1:R]; {P,HU,HC} = GHK(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "GHK: P AND TIME ";; P~EC; @PCF@ PARM = DETW; R = 6*RBASE; U = USSBIG[.,1:R]; {P,HU,HC} = PCF(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "PCF: P AND TIME ";; P~EC; @DCS@ PARM = DETW; R = 2*RBASE; U = USBIG[.,1:R]; {P,HU,HC} = DCS(M,MU,W,WI,C,A,B,R,U,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,U,PARM); "DCS: P AND TIME ";; P~EC; @ARSE@ R = 10+ROUND(RBASE/4); SDD = CINT(UUBIG[1,1]*100000); PARM = CINT(100000*RNDUS(R,1,SDD))|R; {P,HU,HC} = ARSE(M,MU,W,WI,C,A,B,R,UUBIG,PARM); EC = SIM(M,MU,W,WI,C,A,B,R,UUBIG,PARM); "ARSE: P, TIME ";; P~EC; @ARSR@ R = 10+ROUND(RBASE/8); SDD = CINT(UUBIG[1,1]*100000); PARM = CINT(100000*RNDUS(R,1,SDD))|R; {P,HU,HC} = ARSR(M,MU,W,WI,C,A,B,R,UUBIG,PARM); EC = SIM (M,MU,W,WI,C,A,B,R,UUBIG,PARM); "ARSR: P, TIME ";; P~EC; @GSS@ PARM = ROUND(5+RBASE/20); R = PARM; @watch dimension@ {P,HU,HC} = GSS(M,MU,W,WI,C,A,B,R,UUBIG[.,1:R],PARM); EC = SIM(M,MU,W,WI,C,A,B,R,UUBIG[.,1:R],PARM); "GSS: P AND TIME ";; P~EC; @AUS@ R = ROUND(5+RBASE/20); PARM = CINT(100000*RNDU(R,1))~ONES(R,1); PARM[1,2] = ROUND(5+RBASE/20); PARM[2,2] = ROUND(10+RBASE/10); PARM[3,2] = ROUND(1000000*UUBIG[2,1]); {P,HU,HC} = AUS(M,MU,W,WI,C,A,B,R,UUBIG[.,1:R],PARM); EC = SIM(M,MU,W,WI,C,A,B,R,UUBIG[.,1:R],PARM); "AUS: P AND TIME ";; P~EC; PAUSE(5); @COLLECT RESULTS@ " "; FORMAT /M1 /RO 13,6; " PROB. DERIV.MU LOG DER. DERIV. LOG DER. "; " MU LAM LAM "; " "; AA = PP~DM~LM~DL~LL; "EXACT1";; AA[1,.]; "EXACT2";; AA[2,.]; "QUAD. ";; AA[3,.]; "CFS ";; AA[4,.]; "NISE ";; AA[5,.]; "NIST ";; AA[6,.]; "KFS ";; AA[7,.]; "SDS ";; AA[8,.]; "GHK ";; AA[9,.]; "PCF ";; AA[10,.]; "DCS ";; AA[11,.]; "ARSE ";; AA[12,.]; "ARSR ";; AA[13,.]; "GSS ";; AA[14,.]; "AUS ";; AA[15,.]; PAUSE(5); " "; FORMAT /M1 /RD 9,0; BB = (REP~EL~(100*REP./EL))'; " QUAD CFS NISE NIST PKF SDS GHK"; "REPET ";; BB[1,3:9]; "TIME ";; BB[2,3:9]; "REP/SEC";; BB[3,3:9]; " PCF DCS ARSE ARSR GSS AUS"; "REPET ";; BB[1,10:15]; "TIME ";; BB[2,10:15]; "REP/SEC";; BB[3,10:15]; END; PROC 1 = SIM(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL EC; REP = REP|R; PP = PP|P; DM = DM|HU[M+1,1]; LM = LM|HC[M+1,1]; DL = DL|FDL(HU,L); LL = LL|FDL(HC,L); E2 = ETW*TIME; EC = E2-E1; EL = EL|EC; E1 = E2; RETP(EC); ENDP; PROC 3 = FFACT(X,A,B,D,L); @INTEGRAND IN ONE-FACTOR MODEL@ LOCAL XXX,BB,AA,DD,T,TT,Y0,Y1,Y2,TTT; XXX = L*(X'); @ARRAY M BY 40@ BB = (B-MU)*ONES(1,40); AA = (A-MU)*ONES(1,40); DD = (1./D)*ONES(1,40); T = CDFN((BB-XXX).*DD)-CDFN((AA-XXX).*DD); @M BY 40@ Y0 = PRODC(T); @COL VECTOR LENGTH 40@ TT = PRODC(T[1:(M-1),.]); @LENGTH 40@ TTT = PDFN((BB[M,.]-XXX[M,.]).*DD[M,.])-PDFN((AA[M,.]-XXX[M,.]).*DD[M,.]); Y1 = (-1/D[M,1])*(TT.*(TTT')); Y2 = Y1.*X; RETP(Y0,Y1,Y2); ENDP; PROC 1 = FDL(H,L); RETP(H[M+1,2:(M+1)]*L*2); ENDP; PROC 1 = PF(M,AA,BB,MU,W); LOCAL S,SS,R,A,B,P,R,R1,R2,R3,A1,A2,A3,B1,B2,B3; S = SQRT(DIAG(W)); SS = DIAGRV(EYE(M),(1./S)); R = SS*W*SS; A = (AA-MU)./S; B = (BB-MU)./S; IF M == 2; R = R[1,2]; A1 = A[1,1]; A2 = A[2,1]; B1 = B[1,1]; B2 = B[2,1]; P = CDFBVN(B1,B2,R)-CDFBVN(A1,B2,R)-CDFBVN(B1,A2,R)+CDFBVN(A1,A2,R); ELSE; IF M == 3; R1 = R[1,2]; R2 = R[2,3]; R3 = R[3,1]; A1 = A[1,1]; A2 = A[2,1]; A3 = A[3,1]; B1 = B[1,1]; B2 = B[2,1]; B3 = B[3,1]; P = CDFTVN(B1,B2,B3,R1,R2,R3)-CDFTVN(A1,B2,B3,R1,R2,R3); P = P-CDFTVN(B1,A2,B3,R1,R2,R3)-CDFTVN(B1,B2,A3,R1,R2,R3); P = P+CDFTVN(A1,A2,B3,R1,R2,R3)+CDFTVN(B1,A2,A3,R1,R2,R3); P = P+CDFTVN(A1,B2,A3,R1,R2,R3)-CDFTVN(A1,A2,A3,R1,R2,R3); ELSE; P = -999; ENDIF; ENDIF; RETP(P); ENDP; PROC 3 = EXAF(M,A,B,MU,W,L); LOCAL MUP,LP,WP,P0,p1m,p1w,DM,DL; @ vah: -999 check @ MUP = MU[1:M-1,1]|(MU[M,1]+.001); LP = L[1:M-1,1]|(L[M,1]+.001); WP = W - L*(L') + LP*(LP'); P0 = PF(M,A,B,MU,W); p1m=PF(M,A,B,MUP,W); p1w=PF(M,A,B,MU,WP); if p0 ne -999 and p1m ne -999; DM = (p1m-P0)/.001; else; dm = -999; endif; if p0 ne -999 and p1w ne -999; DL = (p1w-P0)/.001; else; dl = -999; endif; RETP(P0,DM,DL); ENDP; proc hsec1990(); /* ** HSEC1990.G ** ** x=hsec1990(); ** ** vah ** ** Returns number of hundredth's of a second since Jan 1, 1990 ** */ local temp,year,month,dom,hsecs,daysinm,doy; temp=date; year=temp[1,1]; month=temp[2,1]; dom=temp[3,1]; hsecs=temp[4,1]; daysinm=31|28|31|30|31|30|31|31|30|31|30|31; @J F M A M J J A S O N D@ if year%4 eq 0; daysinm[2,1]=29; endif; doy=dom; if month gt 1; doy=doy+sumc(daysinm[1:(month-1),1]); endif; retp(((year-1990)*365+doy)*24*3600*100 + hsecs); endp; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ I. Crude Frequency Simulator (CFS) This simulator assumes that U is an array of standard normal variates, independent within each column, and either independent or antithetic across columns. It does not use PARM. It returns both HU and HC. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = CFS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL AA,BB,VD,FF,P,HU,HC; VD = C*U; @ Latent Utility Vectors @ AA = (A-MU)*ONES(1,R); @ mean adjusted lower bound @ BB = (B-MU)*ONES(1,R); @ mean adjusted upper bound @ FF = (VD .< BB).*(VD .> AA); @ Indicators for in rectangles @ FF = MINC(FF); @ Calculate 1(v in B) @ P = SUMC(FF)/R; @ Calculate E1(v in B) @ /* Calculate the matrix HU = E1(v in B)*h(v), formula (6) */ VD = WI*(VD.*(ONES(M,1)*(FF'))); @ Expand into an MxR matrix @ HU = VD*(VD')/R-P*WI; @SE Corner of HU@ VD = SUMC(VD')/R; @SW Corner of HU@ HU = (P|VD)~((VD')|(HU/2)); /* Calculate the matrix HC = HU / P, formula (7) */ IF P > 0; HC = HU/P; ELSE; HC = -999*ONES(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ J. Normal Importance Sampling Simulator (NIS) This simulator provides HC = HU/P. It provides two alternative simu- lators, corresponding to the exponential (NISE) and independent truncated normal (NIST) comparison distributions. These simulators assume U is an array of uniform random (0,1) variates. These procs call the routine NDEN that returns a vector of values of the multivariate normal density. Procedure NISE does not use C or PARM. It calls the bilateral exponential CDF BEXP and its inverse BEXPI from LIB_A.G. It guarantees P positive. Provides both HU and HC. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = NISE(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL S,T1,T2,T3,T4,AA,BB,LAM,PA,PB,VD,D,WGT,VD,P,HU,HC; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ S = SQRT(DIAG(W)); T1 = (B-MU).<(-3*S); T2 = (A-MU).>(3*S); T3 = (B-MU).<(3*S); T4 = (A-MU).>(-3*S); AA = T1.*(B-MU-2*S)+T2.*(A-MU)+(1-T1-T2).*(-3*S+T4.*(A-MU+3*S)); BB = T1.*(B-MU)+T2.*(A-MU+2*S)+(1-T1-T2).*(3*S+T3.*(B-MU-3*S)); S = 1-(BB .< 0)-(AA .> 0); @vah: 1.e-100 added to prevent zero elements in LAM@ LAM = 1.e-100+(BB-AA).*(0.5-S/4)./(DIAG(W)+1.e-100); PA = BEXP(A-MU,LAM)*ONES(1,R); PB = BEXP(B-MU,LAM)*ONES(1,R); {VD,D} = BEXPI(U.*PA+(1-U).*PB,LAM,PA,PB); @EXPONENTIAL VARIATES@ @"lam " lam;"pa " pa;"pb " pb;"prodc(d)" prodc(d);@ WGT = NDEN(M,R,VD,WI,C)./(PRODC(D)+1.e-100); @RATIO OF DENSITIES = @ @ the integrand @ @ Prodc(D) is @ @ formula (23) @ P = SUMC(WGT)/R; @SIMULATE PROB@ S = (ONES(M,1)*(WGT')).*VD; HU = WI*S*(VD')*WI/R - P*WI; @SE CORNER OF H@ VD = WI*SUMC(S')/R; @SW CORNER OF H@ HU = (P|VD)~((VD')|(HU/2)); IF P > 0; HC = HU/P; ELSE; @ HC = ZEROS(M+1,M+1);@ HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /* Procedure NIST guarantees P positive. It calls the proc RNDT for samp- ling from univariate truncated normals.*/ PROC 3 = NIST(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL VD,Y,WGT,VD,P,HU,HC; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ {VD,Y} = RNDT(M,R,U,W,A-MU,B-MU); @TRUNCATED NORMAL VARIATES@ WGT = NDEN(M,R,VD,WI,C)./(Y+1.e-100); @RATIO OF DENSITIES@ @ Y is formula (24) @ P = SUMC(WGT)/R; @SIMULATE PROB@ Y = (ONES(M,1)*(WGT')).*VD; HU = WI*Y*(VD')*WI/R - P*WI; @SE CORNER OF H@ VD = WI*SUMC(Y')/R; @SW CORNER OF H@ HU = (P|VD)~((VD')|HU/2); IF P > 0; HC = HU/P; ELSE; @ HC = ZEROS(M+1,M+1);@ HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ K. Polynomial Kernel-Smoothed Frequency Simulator (PKF) The KFS simulator uses the kernel (31), defined by the function PKER. It assumes that U is an array of standard normal variates, independent within rows. It does not use W. PARM[1,1] contains a window width parameter. This simulator does not return a value for HC. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = PKFS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL AA,BB,VD,F,P,HU,HC; if parm[1,1] le 0; "Error in calling PKFS: smoothing parameter PARM[1,1] MUST be positive"; stop; endif; VD = C*U; @LATENT VECTORS@ AA = (A-MU)*ONES(1,R); BB = (B-MU)*ONES(1,R); F = PKER((VD-BB)/PARM[1,1]) - PKER((VD-AA)/PARM[1,1]); @SMOOTHED INDICATORS@ P = SUMC(F)/R; @FORM SIMULATORS@ AA = WI*VD.*(ONES(M,1)*(F')); HU = WI*VD*(AA')/R-P*WI; @SE CORNER OF H@ VD = SUMC(AA')/R; @SW CORNER OF H@ HU = (P|VD)~((VD')|(HU/2)); IF P > 0; HC = HU/P; ELSE; @ HC = ZEROS(M+1,M+1);@ HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; PROC 1 = PKER(Z); @POLYNOMIAL KERNEL FUNCTION@ LOCAL WGT; WGT = 0.5*(1-Z.*(2-Z)).*(Z.<1)-Z.*Z.*(Z.<0)+0.5*(1+Z.*(2+Z)).*(Z.<-1.); RETP(PRODC(WGT)); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ L. Stern Decomposition Simulator (SDS) The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such that W-lI is positive definite, and assumes that C is a Cholesky factor of W-lI. It assumes that U is an array of standard normal variates, independent within each column. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = SDS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL V,TA,TB,TK,S,J,HT,QQ,DD,TT,P,HU,HC,FB0,FB1,FB2,qi1,qi2; @"SDS: parm[1,1]" parm[1,1];@ if parm[1,1] le 0; "Error in calling SDS: STD parameter PARM[1,1] MUST be positive"; stop; endif; V = C*U; @DRAWS FROM OUTER INTEGRAL@ TB = ((B-MU)*ONES(1,R)-V)/PARM[1,1]; @INNER BOUNDS@ TA = ((A-MU)*ONES(1,R)-V)/PARM[1,1]; tk = -v/parm[1,1]; {FB0,FB1,FB2} = FI(TA,TB,TK); @PARTIAL MOMENTS@ FB0=FB0+.000000001*(FB0 .== 0); qi1 = PARM[1,1].*fb1; qi2 = PARM[1,1]^2 .* fb2; S = qi1./FB0; @EQ (35)@ J = 0; HU = ZEROS(M+1,M+1); DO WHILE J < R; J = J+1; @"j" j "fb0[.,j]" fb0[.,j];@ @"j" j "qi1[.,j]" qi1[.,j];@ @"j" j "qi2[.,j]" qi2[.,j];@ @"j" j "s[.,j]" s[.,j];@ QQ = PRODC(FB0[.,J]); @EQ (36)@ DD = (qi2[.,J]./FB0[.,J]) - (qi1[.,j]./fb0[.,j])^2; DD = DIAGRV(EYE(M),DD); @EQ (38)@ TT = S[.,J]; @EQ (39)@ HT = DD-W+TT*TT'; HT = QQ*((1|TT)~(TT'|(HT/2))); HU = HU+HT; ENDO; S = (1|ZEROS(M,1))~(ZEROS(1,M)|WI); HU = S*HU*S/R; P = HU[1,1]; IF P > 0; HC = HU/P; ELSE; HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ M. Geweke-Hajivassiliou-Keane Simulator (GHK) The GHK procedure assumes that U is an array of uniform [0,1] variates, independent within columns, and in general independent between columns. It does not use PARM. It returns both HU and HC. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = GHK(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL J,II,TA,TB,TT,WGT,V,HU,P,HC; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ J = 1; II = 1; TA = CDFN((A[1,1]-MU[1,1])/(C[1,1]+1.e-100))*ONES(1,R); @1ST COMPONENT BOUNDS@ TB = CDFN((B[1,1]-MU[1,1])/(C[1,1]+1.e-100))*ONES(1,R); /*TT = CDFNI(U[1,.].*TA+(1-U[1,.]).*TB);*/ TT = CDFINVN(U[1,.].*TA+(1-U[1,.]).*TB); WGT = TB-TA; DO WHILE J < M; @LOOP OVER COMPONENTS@ J = J+1; @COMPONENT J BOUNDS@ TA = CDFN(((A[J,1]-MU[J,1])*ONES(1,R)-C[J,II]*TT)/(C[J,J]+1.e-100)); TB = CDFN(((B[J,1]-MU[J,1])*ONES(1,R)-C[J,II]*TT)/(C[J,J]+1.e-100)); /* TT = TT|CDFNI(U[J,.].*TA+(1-U[J,.]).*TB); @COMPONENT J DRAW@*/ TT = TT|CDFINVN(U[J,.].*TA+(1-U[J,.]).*TB); @COMPONENT J DRAW@ II = II|J; WGT = WGT.*(TB-TA); @ ACCUMULATE WEIGHT@ ENDO; V = C*TT; @TRANSFORM TO V-MU@ TT = (ONES(M,1)*WGT).*V; P = SUMC(WGT')/R; HU = WI*V*(TT')*WI/R-P*WI; V = WI*SUMC(TT')/R; HU = (P|V)~((V')|(HU/2)); IF P > 0; HC = HU/P; ELSE; HC = -999*ONES(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ N. Parabolic Cylinder Function Simulator (PCF) The PCF proc assumes that U is an array whose column vectors are points in the intersection of the unit sphere and the negative orthant, and that PARM[1,1] = DET(W). ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = PCF(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL KP,BB,AA,VWV,BWV,S,RO,Y0,Y1,Y2,TT,HU,P,HC; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ /*vah: KP=1/(GAMMA(M/2)*(2^(M/2-1))*SQRT(DET(W))+1.e-100);*/ KP=1/(GAMMA(M/2)*(2^(M/2-1))*SQRT(parm[1,1])+1.e-100); BB = WI*(MU-B); AA = WI*U; VWV = SUMC(U.*AA); BWV = (U')*BB; S = -(MU-B)'*BB/2; /*vah: S = ONES(ROWS(BWV),1)*S;*/ S = ONES(R,1)*S; S = S+BWV.*BWV./(2*VWV+1.e-100); S = KP*EXP(S); RO = MINC(((A-B)*ONES(1,R))./(U+1.e-100)); @UPPER BOUND, EQ (46)@ {Y0,Y1,Y2} = CF(M-1,VWV,BWV,RO); @CYLINDER FUNCTION VALUES@ P = MEANC(Y0.*S); TT = MEANC(U'.*((Y1.*S)*ONES(1,M))); HU = WI*TT*(BB'); TT = WI*TT-P*BB; HU = -HU-(HU')+AA*((AA').*((Y2.*S)*ONES(1,M)))/R; HU = HU+P*(BB*BB'-WI); HU = (P|TT)~(TT'|HU/2); HU = HU/(2^M); @AVERAGE HU@ P = HU[1,1]; IF P > 0; HC = HU/P; ELSE; HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; PROC 3 = CF(K,ALP,BET,GAM); @PARABOLIC CYLINDER FN, EQ (48)@ LOCAL T,TT,Y0,Y1,Y2,Z0,Z1,Z2,GJ,J; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ alp=alp+1.e-100; T = SQRT(ALP); TT = EXP(-0.5*ALP.*(GAM-BET./ALP)^2)./ALP; @FIRST 3 FUNCTIONS@ @ y's are the C's in eq(49-51) @ Z0 = 2.506628274631*(CDFN(BET./T))./T; Y0 = -2.506628274631*(CDFN((BET-ALP.*GAM)./T))./T+Z0; Z1 = (EXP(-BET.*BET./(2*ALP)))./ALP; Y1 = Z1-TT+Y0.*BET./ALP; Z1 = Z1+Z0.*BET./ALP; Z2 = Z1.*BET./ALP+Z0./ALP; Y2 = Y1.*BET./ALP+Y0./ALP-GAM.*TT; IF K > 0; @RECURSE IF NECESSARY@ GJ = GAM; J = 2; DO WHILE J < (K+2); GJ = GJ.*GAM; Y0 = Y1; Y1 = Y2; Y2 = Y1.*BET./ALP+Y0*J./ALP-GJ.*TT; Z0 = Z1; Z1 = Z2; Z2 = Z1.*BET./ALP+Z0*J./ALP; J = J+1; ENDO; ENDIF; RETP(Y0,Y1,Y2); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ O. Deak Chi-Square Simulator (DCS) The DCS assumes that U is an array whose columns are points in the unit sphere, antithetic across columns, and that PARM[1,1] = DET(W). ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = DCS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL K,K1,k2,UU,SN,SP,BB,AA,QL,QU,J,HU,T,TR,TT,D0,D1,D2,HC,HT,P,gammo2; local timelim; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ /* "m" m;"mu" mu';"w";w;"wi";wi;"c";c;"a" a';"b" b';"r" r; call msize("u",u); "u";u;"parm" parm; */ /*vah: K = 1@2@/SQRT(DET(W)+1.e-100); @LEAD CONSTANT@*/ K = 1@2@/SQRT(parm[1,1]+1.e-100); @LEAD CONSTANT@ gammo2=GAMMA(M/2)+1.e-100; K1 = K*GAMMA((M+1)/2)*SQRT(2)/gammo2; k2 = k*2*gamma((m+2)/2)/gammo2; UU = U+(1.E-100)*(U .== 0.); @FUZZ EXACT ZEROS@ SN = U .< 0; @LOCATE NEGATIVES@ SP = U .> 0; BB = (B-MU)*ONES(1,R); AA = (A-MU)*ONES(1,R); QL = MAXC(((AA.*SP+BB.*SN)./(UU+1.e-100))|ZEROS(1,R)); @EQ (58)@ QU = MINC((BB.*SP+AA.*SN)./(UU+1.e-100)); QU = QL+(QU.>QL).*(QU-QL); @CONSISTENT BOUNDS@ J = 0; HU = ZEROS(M+1,M+1); /*timelim=hsec1990();*/ DO WHILE J < R; @ACCUM. HU OVER REPETITIONS@ J = J+1; T = UU[.,J]'*WI*UU[.,J]; TR = 1/(SQRT(T)+1.e-100); TT = TR^M; IF QU[J] > QL[J]; @EQ (57)@ D0 = K*TT*(CDFCHIC(T*QL[J]^2,M)-CDFCHIC(T*QU[J]^2,M)); D1 = K1*TT*TR*(CDFCHIC(T*QL[J]^2,M+1)-CDFCHIC(T*QU[J]^2,M+1)); @ D2 = K*2*(M/2)*TT*TR*TR*(CDFCHIC(T*QL[J]^2,M+2)-CDFCHIC(T*QU[J]^2,M+2));@ D2 = K2*TT*TR*TR*(CDFCHIC(T*QL[J]^2,M+2)-CDFCHIC(T*QU[J]^2,M+2)); ELSE; D0 = 0; D1 = 0; D2 = 0; ENDIF; TR = WI*U[.,J]; @this is v' in (59)@ HT = D2*(TR*(TR'))-D0*WI; @EQ (59)@ HT = (D0|(D1*TR))~((D1*TR')|(HT/2)); HU = HU+HT; /* if (hsec1990()-timelim) gt hseclim;*/ /* @give up@*/ /* r=j;*/ /* j=r;*/ /* endif;*/ ENDO; HU = HU/R; P = HU[1,1]; IF P > 0; HC = HU/P; ELSE; HC = -999*ones(M+1,M+1); ENDIF; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ P. Acceptance/Rejection Simulator (ARS) The ARS simulator is given in a version using the truncated bilateral exponential comparison distribution (ARSE), and in a version using a recursive truncated comparison distribution (ARSR). PARM[.,1] is assumed to contain seeds for the uniform random number generator RNDUS called by the routines; there is one seed for each repetition R. PARM[R+1,1] is assumed to contain a censoring limit for rejection. Should only compute HC. The ARS method provides a mechanism for drawing from a conditional density when transformations from uniform or standard normal variates are not available. We want to sample from f(x/a). We generate x from g(x). In our case g(x) = ARSE. Generate tsi from uniform. Generate unconditional f(x) and bound alpha. We accept x if tsi < f(x)/(g(x)*alpha) x has conditional density f(x/a) = HC. We set x = VD ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = ARSE(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL ALP,JL,Q,SD,I,U,D,V,J,HT,S,T1,T2,T3,T4,AA,BB,LAM,PA,PB,D,WGT,VD,P,HU,HC; local timebeg; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ Q = 0; @ set arse proportion @ @ accepted @ VD = EIGRS(WI); @ eigenroots of covariance @ S = SQRT(DIAG(W)); @ matrix @ T1 = (B-MU).<(-3*S); T2 = (A-MU).>(3*S); T3 = (B-MU).<(3*S); T4 = (A-MU).>(-3*S); AA = T1.*(B-MU-2*S)+T2.*(A-MU)+(1-T1-T2).*(-3*S+T4.*(A-MU+3*S)); BB = T1.*(B-MU)+T2.*(A-MU+2*S)+(1-T1-T2).*(3*S+T3.*(B-MU-3*S)); S = 1-(BB .< 0)-(AA .> 0); LAM = (BB-AA).*(0.5-S/4)./(DIAG(W)+1.e-100); ALP = LAM'*LAM/(2*VD[1,1]+1.e-100); @ bound on density ratio @ @ equal to alpha @ PA = BEXP(A-MU,LAM); @ lower bound for truncation @ @ to get x = VD from ARSE @ PB = BEXP(B-MU,LAM); @ upper bound for truncation @ HC = ZEROS(M+1,M+1); JL = PARM[R+1,1]; @ Censoring limit for rejection @ timebeg=hsec1990(); I = 0; DO WHILE I < R; @LOOP OVER REPETITIONS@ I = I+1; J = 0; SD = PARM[I,1]; @SEEDS@ DO WHILE J < JL; @DRAW UNTIL CENSOR LIMIT@ J = J+1; U = RNDUS(M,1,SD); {VD,D} = BEXPI(U.*PA+(1-U).*PB,LAM,PA,PB); @EXPONENTIAL VARIATES@ @ the x drawn from g @ U = RNDUS(1,1,SD); @ Draw from uniform @ @ so u = tsi @ IF (LN(U)) < (-(VD')*WI*VD/2+(LAM')*ABS(VD)-ALP); @ ACCEPT x = VD ? @ /* LN(tsi) < LN(f(x)/g(x)*alpha) */ @ Note f(x) is normal @ @ and g(x) is BEXP @ Q = Q+1/R; J = JL; ENDIF; if (hsec1990()-timebeg) gt hseclim; @give up@ r=i; jl=j; "ARSE: giving up because time-limit exceeded. Completed R and J=" r j; endif; ENDO; /* Now find HC using formula (7) */ VD = WI*VD; HT = (VD*(VD')-WI)/2; @ACCUMULATE HC@ HC = HC + (1|VD)~(VD'|HT); ENDO; HC = HC/R; P = -999; @HU = ZEROS(M+1,M+1);@ HU = -999*ones(M+1,M+1); "ARSE PROPORTION NON-CENSORED";; Q; RETP(P,HU,HC); ENDP; /* The following uses a recursive truncated */ /* comparison distribution for g(x) */ PROC 3 = ARSR(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL HC,JL,Q,SD,I,J,K,II,TA,TB,UU,HT,WGT,ALP,TT,V,P,HU; local timebeg; HC = ZEROS(M+1,M+1); JL = PARM[R+1,1]; @CENSORING LIMIT@ Q = 0; timebeg=hsec1990(); I = 0; DO WHILE I < R; @LOOP OVER REPETITIONS@ I = I+1; J = 0; SD =PARM[I,1]; @SEEDS@ DO WHILE J < JL; @DRAW UNTIL ACCEPT OR LIMIT@ J = J+1; K = 1; @START RECURSIVE DRAW@ II = 1; TA = CDFN((A[1,1]-MU[1,1])/C[1,1]); @INITIAL PROB BOUNDS@ TB = CDFN((B[1,1]-MU[1,1])/C[1,1]); UU = RNDUS(M+1,JL,SD); /* TT = CDFNI(UU[1,J]*TA+(1-UU[1,J])*TB); @1ST COMPONENT@*/ TT = CDFINVN(UU[1,J]*TA+(1-UU[1,J])*TB); @1ST COMPONENT@ WGT = TB-TA; @DENSITY RATIO@ @ ie f(x)/g(x) @ ALP = WGT; @BOUND ON DENSITY RATIO@ @ equal to alpha @ DO WHILE K < M; K = K+1; TA = CDFN(((A[K,1]-MU[K,1])-C[K,II]*TT)/C[K,K]); TB = CDFN(((B[K,1]-MU[K,1])-C[K,II]*TT)/C[K,K]); /* TT = TT|(CDFNI(UU[K,J].*TA+(1-UU[K,J]).*TB));*/ TT = TT|(CDFINVN(UU[K,J].*TA+(1-UU[K,J]).*TB)); II = II|K; WGT = WGT*(TB-TA); @UPDATE DENSITY RATIO@ ALP = ALP*(2*CDFN((B[K]-A[K])/(2*C[K,K]))-1); @UPDATE BOUND@ ENDO; IF UU[M+1,J] < WGT/ALP; @ accept if tsi < (f(x)/g(x))/alpha @ Q = Q+1/R; J = JL; ENDIF; if (hsec1990()-timebeg) gt hseclim; @give up@ r=i; jl=j; "ARSR: giving up because time-limit exceeded. Completed R and J=" r j; endif; ENDO; /* Find HC using formula (7) */ V = WI*C*TT; @ACCUMULATE HC@ HT = V*V'-WI; HT = (1|V)~(V'|HT/2); HC = HC+HT; ENDO; HC = HC/R; P = -999; @HU = ZEROS(M+1,M+1);@ HU = -999*ones(M+1,M+1); "ARSR PROPORTION NON-CENSORED";; Q; RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Q. Gibbs's Sampler Simulator (GSS) The GSS simulates HC, but does not simulate P or HU. It assumes that PARM[1,1] gives the number of rounds to be used in the construction of the random draws (JL). It assumes that U is an MxJL array of uniform [0,1] random numbers. It returns only HC. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = GSS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL HC,JL,I,J,K,MK,TA,TB,E,F,G,D,S,HT,V,P,HU,UI,UIKJ; local wimkk,wmkmki,wmkk,vbar,vvpbar,condmukf,condsdk,csdk,cmuk,seedgss; local timelim; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ ui=u; condmukf=zeros(m,m-1); condsdk =zeros(m,1); k=1; do until k>m; IF K == 1; @INDICES EXCLUDING K@ MK = SEQA(2,1,M-1); ELSE; IF K == M; MK = SEQA(1,1,M-1); ELSE; MK = SEQA(1,1,K-1)|SEQA(K+1,1,M-K); ENDIF; ENDIF; wimkk = WI[MK,K]; wmkk = W[MK,K]; wmkmki = wi[mk,mk] - (wimkk*wimkk')/(wi[k,k]+1.e-100); @INVERSE OF W EXCL.ROW/COL K@ condsdk[k,1] = sqrt(w[k,k] - (wmkk')*wmkmki*wmkk); @CONDIT. SD OF V[K]@ condmukf[k,.] = (wmkk')*wmkmki; @FACTOR FOR CONDIT. MEAN OF V[K]-MU[K]@ k=k+1; endo; /*begin*/ vbar=0; vvpbar=0; JL = PARM[1,1]; @NUMBER OF GIBBS ROUNDS@ I = 0; V = (B+A)/2-MU; @INITIAL VECTOR@ /*timelim=hsec1990();*/ DO WHILE I < R; @LOOP OVER REPETITIONS@ I = I+1; if i > 1; seedgss=ui[1,1]*1.e+9; ui=rndus(m,jl,seedgss); endif; J = 0; DO WHILE J < JL; @DRAW UNTIL LIMIT@ J = J+1; K = 0; @START RECURSIVE STEPS@ DO WHILE K < M; @STEP OVER K@ K = K+1; IF K == 1; @INDICES EXCLUDING K@ MK = SEQA(2,1,M-1); ELSE; IF K == M; MK = SEQA(1,1,M-1); ELSE; MK = SEQA(1,1,K-1)|SEQA(K+1,1,M-K); ENDIF; ENDIF; csdk=condsdk[k,1]; cmuk=condmukf[k,.]*v[mk,1]; TA = CDFN((A[K]-MU[K]-cmuk)/(csdk+1.e-100)); @BOUNDS@ TB = CDFN((B[K]-MU[K]-cmuk)/(csdk+1.e-100)); UIKJ = UI[K,J]; @"ta tb uikj" ta tb uikj k;@ V[K] = cmuk + csdk*CDFINVN(UIKJ*TA+(1-UIKJ)*TB); @DRAW FROM TRUNCATED N@ ENDO; ENDO; vbar=vbar+v; vvpbar=vvpbar+v*v'; /* if (hsec1990()-timelim) gt hseclim;*/ /* @give up@*/ /* r=j;*/ /* j=r;*/ /* endif;*/ ENDO; vbar=wi*vbar/r; vvpbar=wi*(vvpbar/r - w)*wi/2; hc=(1|vbar)~(vbar'|vvpbar); P = -999; HU = -999*ones(M+1,M+1); RETP(P,HU,HC); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ S. Asymptotically Unbiased Simulator (AUS) The AUS simulator is obtained as the product of an unbiased simulator of H and an asymptotically unbiased simulator of 1/P. The proc is defined with the GHK simulator used for H and for the initial independent simulations of P(B)/P(A), where A is the event that the first component is in the rectangle bounds. The proc assumes that PARM(.,1) contains a vector of R seeds for the normal random number generator RNDNS, PARM(1,2) contains the number of initial GHK simulators of P used in the simulation of 1/P, PARM(2,2) contains the maximum number of trials using the crude frequency simulator that are made before the process is censored. If PARM(2,2) is made large (say 1.e50), then this is computationally the same as the Sequential Unbiased Simulator (SUS). PARM(3,2) contains a seed for RNDUS. Gives HU and HC. This does continous parameter estimation of 1/p. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = AUS(M,MU,W,WI,C,A,B,R,U,PARM); LOCAL F,G,TA,TB,A2,B2,JLIM,SD,SDN,II,HU,PA,PT,PINV,Q,J,PPINV,HC,I,RR,UU,T,TT,P; local timebeg; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ F = W[2:M,1]/(W[1,1]+1.e-100); G = C[2:M,2:M]; @CHOL FACTOR OF W[2:M,2:M]@ TA = CDFN((A[1,1]-MU[1,1])/(C[1,1]+1.e-100)); @BOUNDS FOR 1ST COMPONENT@ TB = CDFN((B[1,1]-MU[1,1])/(C[1,1]+1.e-100)); A2 = A[2:M]-MU[2:M]; @BOUNDS FOR REMAINING COMPONENTS@ B2 = B[2:M]-MU[2:M]; JLIM = PARM[2,2]; @MAX NUMBER OF CFS DRAWS@ /*"in AUS: Max num of CFS draws= " jlim;*/ SD = PARM[3,2]; @SEED FOR RNDUS@ II = 2; {PA,PT,HU} = FGHK(II,M,MU,W,WI,C,A,B,R,SD); @NUMERATOR@ @ i.e. H=HU @ PINV = 1; Q = 1; J = 0; II = 1; DO WHILE J < PARM[1,2]; @INITIAL TERMS FOR DENOMINATOR@ J = J+1; {PA,RR,SDN} = FGHK(II,M,MU,W,WI,C,A,B,R,SD); @P(A)@ Q = Q.*(1-PT./(PA+1.e-100)); PINV = PINV+Q; ENDO; @COMPLETES PIECES COMMON TO ALL REPETITIONS@ I = 0; PPINV = 0; timebeg=hsec1990(); DO WHILE I < R; @LOOP OVER REPETITIONS@ I = I+1; RR = 0; @START FOR DENOM SEQ@ J = 0; SDN = PARM[I,1]; DO WHILE J < JLIM; @SEQUENTIAL DRAWS UNTIL STOP@ J = J+1; UU = RNDUS(1,1,SD); @DRAW FROM A@ /* TT = C[1,1]*CDFNI(UU*TA+(1-UU)*TB); @1ST COMPONENT DRAW@*/ TT = C[1,1]*CDFINVN(UU*TA+(1-UU)*TB); @1ST COMPONENT DRAW@ T = F*TT+G*RNDNS(M-1,1,SDN); @REMAINING COMPONENTS DRAW@ IF MINC(T-A2) > 0 AND MAXC(T-B2) < 0; @ if E1(v not in B) @ J = JLIM; ELSE; RR = RR+1; @ACCUMULATE CFS RUN LENGTH@ ENDIF; @"aus: hsec1990()-timebeg, hseclim, i j" hsec1990()-timebeg hseclim i j;@ if (hsec1990()-timebeg) gt hseclim; @give up@ r=i; jlim=j; "AUS: giving up because time-limit exceeded. Completed R and J=" r j; endif; ENDO; PPINV = PPINV+PINV+Q*RR; @ACCUMULATE EST. OF P(B)/P(A)@ ENDO; PPINV=PPINV/(R*PA+1.e-100); @ = 1/P @ HC = HU*PPINV; @ H * 1/P @ P = HU[1,1]; RETP(P,HU,HC); ENDP; /* Gives unbiased simulator of H=HU */ PROC 3 = FGHK(I,M,MU,W,WI,C,A,B,R,SEED); @VERSION OF GHK SIMULATOR USED IN AUS, SEE GHK FOR COMMENTS@ LOCAL J,II,TA,TB,TT,WGT,PA,V,HU,P,U; J = 1; II = 1; U = RNDUS(M,R,SEED); TA = CDFN((A[1,1]-MU[1,1])/(C[1,1]+1.e-100))*ONES(1,R); @1ST COMPONENT BOUNDS@ TB = CDFN((B[1,1]-MU[1,1])/(C[1,1]+1.e-100))*ONES(1,R); /*TT = CDFNI(U[1,.].*TA+(1-U[1,.]).*TB); @ inverse cumulative @*/ TT = CDFINVN(U[1,.].*TA+(1-U[1,.]).*TB); @ inverse cumulative @ @ normal @ PA = TB-TA; WGT = PA; DO WHILE J < M; @LOOP OVER COMPONENTS@ J = J+1; @COMPONENT J BOUNDS@ TA = CDFN(((A[J,1]-MU[J,1])*ONES(1,R)-C[J,II]*TT)/(C[J,J]+1.e-100)); TB = CDFN(((B[J,1]-MU[J,1])*ONES(1,R)-C[J,II]*TT)/(C[J,J]+1.e-100)); /* TT = TT|CDFNI(U[J,.].*TA+(1-U[J,.]).*TB); @COMPONENT J DRAW@*/ TT = TT|CDFINVN(U[J,.].*TA+(1-U[J,.]).*TB); @COMPONENT J DRAW@ II = II|J; WGT = WGT.*(TB-TA); @ ACCUMULATE WEIGHT@ ENDO; P = MEANC(VEC(WGT)); IF I == 1; HU = 0; ELSE; V = C*TT; TT = (ONES(M,1)*WGT).*V; HU = WI*V*(TT')*WI/R-P*WI; V = WI*SUMC(TT')/R; HU = (P|V)~((V')|(HU/2)); ENDIF; PA = MEANC(VEC(PA)); RETP(PA,P,HU); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ R. Sequentially Unbiased Simulator (SUS) The SUS simulator is obtained as the product of an unbiased simulator of H and an unbiased simulator of 1/P. Computationally, if PARM(2,2) is made large (say 1.e50), SUS is the same as the Asymptotically Unbiased Simulator (AUS). ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 3 = SUS(M,MU,W,WI,C,A,B,R,U,PARM); local m,mu,w,wi,c,a,b,r,u,parm,p,hu,hc; parm[2,2]=1.e50; {p,hu,hc}=aus(M,MU,W,WI,C,A,B,R,U,PARM); retp(p,hu,hc); endp; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ T. Statistical Functions (LIB_A) This paragraph contains a series of procedures that return statistical functions. NDEN is the multivariate normal density, CDFNI and CDFINVN are two different algorithms for calculating the inverse of the cumulative standard univariate normal, RNDT is a proc for drawing an array from truncated independent normals that have the same means and variances as the multivariate density, and FI is a proc that returns partial moments of the standard normal. CDFINVN appears faster and more accurate than CDFNI. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 1 = NDEN(M,R,VD,WI,C); @MULTIVARIATE NORMAL DENSITY, DIMENSION M@ /* WI = INVERSE OF COV MATRIX, C = CHOLESKY FACTOR, VD = M BY R ARRAY OF DEVIATIONS FROM MEAN, EACH COL. A POINT, R = NUMBER OF POINTS PROCESSED SIMULTANEOUSLY */ LOCAL KK,E,I,Q; KK = 0.39894228040143^M; KK = KK/(PRODC(DIAG(C))); @CONSTANT@ E = WI*VD; @PART OF QUADRATIC FORM@ Q = KK*EXP(-0.5*SUMC(VD.*E)); @DENSITY@ RETP(Q); ENDP; /***CDFINVN used instead of CDFNI: PROC 1 = CDFNI(P); @INVERSE CUMULATIVE NORMAL@ LOCAL S,Z,T,X0,X; S = P .> 0.5; @P > .5?@ Z = S.*(1-P)+(1-S).*P; @IF SO, USE 1-P@ T = SQRT(LN(1/(Z.*Z))); @ABRAMOWITZ & STEGUM@ X0 = 2.515517+T.*(0.802853+T*0.010328); @POLYNOMIAL APPROX.@ X0 = X0./(1+T.*(1.432788+T.*(0.189269+T*0.001308))); X0 = X0-T; X = X0 + (Z-CDFN(X0))./PDFN(X0); @NEWTON STEP@ X = (1-2*S).*X; @FIX SIGN OF X@ RETP(X); ENDP; not used***/ /* ** CDFINVN -- computes the inverse normal cdf. The algorithm used is taken ** from Kennedy and Gentle, STATISTICAL COMPUTING, p. 95. ** ** vah 22:57:35.78, Sun, Jul 23, 1989 ** Taken intact from V149b MAXPROC3.ARC ** ** xp=cdfinvn(p); ** ** INPUT ** p -- nx1 vector of probabilities (that is, 0 < p < 1) ** OUTPUT ** xp -- the corresponding vector of values such that: cdfn(xp)=p; ** ** This algorithm seems to be very accurate. It is ** rated to be accurate to 1e-7. However, it appears ** to be substantially more accurate -- approximately ** 1e-10. ** */ proc 1=cdfinvn(p); local lim, p0, p1, p2, p3, p4, q0, q1, q2, q3, q4, maskgt, maskeq, sgn, y, xp, pn; @ Constants @ lim = 1e-20; p0 = - 0.322232431088; q0 = 0.0993484626060; p1 = - 1.0; q1 = 0.588581570495; p2 = - 0.342242088547; q2 = 0.531103462366; p3 = - 0.0204231210245; q3 = 0.103537752850; p4 = - 0.453642210148*1e-4; q4 = 0.38560700634*1e-2; @ Main body of code @ @ Create masks for handling p > .5 and p >= .5 @ maskgt=(p .> 0.5); maskeq=(p .ne 0.5); sgn=missrv(miss(maskgt,0),-1); @ Convert p > .5 --> 1-p @ pn=(maskgt-p).*sgn; clear maskgt; @ Computation of function for p < 0.5 @ y=sqrt(-2*ln(pn)); clear pn; xp=y + ((((y*p4 + p3).*y + p2).*y + p1).*y + p0)./ ((((y*q4 + q3).*y + q2).*y + q1).*y + q0); clear y; @ Convert results for p > .5 and p = .5 @ xp=(xp.*sgn).*maskeq; clear sgn, maskeq; retp(xp); endp; PROC 2 = RNDT(M,R,U,W,A,B); /* RETURNS M BY R ARRAY OF INDEP. TRUNCATED NORMAL DEVIATES (IN V) AND ASSOCIATED R BY 1 VECTOR OF LEVELS OF DENSITY (IN Y), DRAWN FROM INDEP. NORMALS WITH SAME VARIANCES AS THE MULTIVARIATE NORMAL DENSITY. U IS ASSUMED TO BE A M BY R ARRAY OF UNIFORM (0,1) VARIATES. */ LOCAL S,YA,YB,VD,V,Y,KK; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ S = SQRT(DIAG(W))+1.e-100; @STANDARD DEVIATIONS@ YB = CDFN(B./S)*ONES(1,R); @BOUNDS@ YA = CDFN(A./S)*ONES(1,R); /*VD = CDFNI(U.*YB+(1-U).*YA); @DRAWS@*/ VD = CDFINVN(U.*YB+(1-U).*YA); @DRAWS@ V = (S*ONES(1,R)).*VD; @TRUNCATED NORMAL DENSITY@ KK = 0.39894228040143; Y = ((KK ./ S)*ONES(1,R)).*EXP(-0.5*VD.*VD)./(YB-YA+1.e-100); Y=PRODC(Y); RETP(V,Y); ENDP; PROC 3 = FI(A,B,KAP); /*FIRST THREE PARTIAL MOMENTS OF (Y-KAP)^K FROM A TO B FOR STANDARD NORMAL. CALLS CDFN AND PDFN*/ LOCAL FI0,FI1,FI2; FI0 = CDFN(B)-CDFN(A); FI1 = -KAP.*FI0 - PDFN(B) + PDFN(A); FI2 = FI0 - KAP.*FI1 - (B-KAP).*PDFN(B) + (A-KAP).*PDFN(A); RETP(FI0,FI1,FI2); ENDP; PROC 1 = BEXP(V,LAM); @BILATERAL EXPONENTIAL CDF@ LOCAL VV,VVV; VV = - ABS(V); VVV = V .> 0; RETP(VVV+(.5-VVV).*EXP(LAM.*VV)); ENDP; PROC 2 = BEXPI(P,LAM,PA,PB); @INVERSE TRUNCATED BILATERAL@ LOCAL PP,X,D; @EXPONENTIAL CDF@ PP = P+(P.>0.5).*(1-2*P); @vah: prevent numerical problems in difficult cases @ @LN() : add 1.e-100 so as to prevent problems if PP comes out equal to 0@ @./LAM : add 1.e-100 in case LAM has 0 element(s) @ @./(PB-PA): add 1.e-100 in case PB has equal element(s) with PA @ X = (1-2*(P.>0.5)).*LN(2*PP+1.e-100)./(LAM+1.e-100); D = LAM.*PP./(PB-PA+1.e-100); RETP(X,D); ENDP; /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ U. Spherical Transformations and Antithetics (LIB_B) RNDBASE is a routine that draws a random basis in K space, where K is the dimension of the random coefficient vector. ANTITH is a routine that creates a "grid" of antithetic points S on the unit sphere, starting from a random basis provided by RNDBASE. These routines will ordinarily be called for each observation. To avoid "chatter" when parameter values are changed, the seed for the random number generator for each observation should be stored so that the same outputs can be regenerated in iteration to parameter estimates. RNDBASE takes as input a dimension K and returns a K by K array V that is column orthonormal and random, bordered by a column of random lengths of K-dim. normal random vectors. SEED1 is a seed for the random number generator. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ PROC 2 = RNDBASE(K,SEED1); LOCAL S,KK,I,V,T; @vah: guard against possible division by 0, by adding 1.e-100 to denominators@ V = RNDNS(K,K,SEED1); @ GET RANDOM ARRAY @ T = DIAG(V*V'); @ GET RANDOM LENGTHS @ RFIRST: S = V[.,1]'V[.,1]; IF S < 1.0E-16; @ IF NEAR ZERO, REDRAW @ V[.,1] = RNDNS(K,1,SEED1); GOTO RFIRST; ENDIF; V[.,1] = V[.,1]/(SQRT(S)+1.e-100); @ NORMALIZE @ KK = 2; DO UNTIL KK > K; @ LOOP OVER COLUMNS 2 TO K @ RNEXT: I = 0; DO WHILE I < KK; @ ORTHOGONALIZE TO EARLIER @ I = I+1; V[.,KK] = V[.,KK]-V[.,I]*(V[.,KK]'V[.,I]); ENDO; S = V[.,KK]'V[.,KK]; IF S < 1.0E-16; @ IF NEAR ZERO, REDRAW @ V[.,KK] = RNDNS(K,1,SEED1); GOTO RNEXT; ENDIF; V[.,KK] = V[.,KK]/(SQRT(S)+1.e-100); @ NORMALIZE @ KK = KK+1; ENDO; RETP(V,T); @ RETURN RANDOM BASIS AND LENGTHS@ ENDP; /* Starting from a random basis B, obtained from RNDBASE with seed SEED, the ANTITH routine returns an array A with columns of length K that have an antithetic property. R is an integer that determines the number of columns M as follows: If T <= 1, then M = 2K, and A contains the columns of B and their antipodes. If T > 1, then M = 2K(1+(T-1)(K-1)), obtained by taking 4T equally spaced points on each of the K(k-1)/2 great circles through pairs of directions in the random basis.*/ PROC ANTITH(K,T,V0); LOCAL B,J,I,II,JJ,V1,V2,S,A; B = V0[.,1:K]; @ START FROM RANDOM BASIS @ A = B~(-B); @ BASIS AND ANTIPODES @ IF T <= 1; GOTO ALAST; @ STOP IF ENOUGH POINTS @ ELSE; J = 2 + INT((T/2-K)/(2*K*(K-1))); @ DETERMINE # PTS NEEDED @ I = 1; DO UNTIL I == K; @ LOOP ON 1ST CIRCLE PT @ II = I; DO WHILE II < K; @ LOOP ON 2ND CIRCLE PT @ II = II+1; JJ = 1.; DO WHILE JJ < J; @ LOOP ON INTERMED. PTS @ V1 = B[.,I]*COS(1.5707963*JJ/J); V2 = B[.,II]*SIN(1.5707963*(1.-JJ/J)); A = A~(V1+V2)~(V1-V2)~(-V1+V2)~(-V1-V2); @ ADD NEW PTS TO A @ JJ = JJ+1.; ENDO; ENDO; I = I+1; ENDO; ENDIF; ALAST: RETP(A); ENDP;