program harness c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c A. Overview c This section contains a series of FORTRAN procedures for c simulation of multivariate normal rectangle probabilities, and the c derivatives of these probabilities with respect to the mean and c covariances of the normal distribution. Table 2 indexes the c procedures in this section. c c TABLE 2. FORTRAN PROCEDURES c c SECTION MNEUMONIC* DESCRIPTION c c H CFS Crude Frequency Simulator c I NIS Normal Importance Sampling Simulator c Procs NISE and NIST use exponential and c truncated normal comparison c distributions, respectively c J KFS Polynomial Kernel Smoothed Frequency c Simulator. The proc is PKFS c K SDS Stern Decomposition Simulator c L GHK Geweke-Hajivassiliou-Keane Simulator c M PCF Parabolic Cylinder Function Simulator c N DCS Deak chi-square Simulator c O ARS Acceptance/Rejection Simulator c Procs ARSE and ARSR use exponential and c recursive truncated normal comparison c distributions, respectively c P GSS Gibbs Sampler Simulator c Q AUS Approximate (OR Sequential) Unbiased c Simulator c R SUS Sequential Unbiased Simulator c S LIB_A Statistical Functions c T LIB_B Spherical Transformations and c Antithetics c U LIB_C Miscellaneous Utilities c c * The procs are located in libraries named *.G, where '*' denotes the c mneumonic. c c Each of the simulator procedures requires the following standard c inputs; the interpretation of some inputs may vary from routine to c routine, and not all are used in some routines: c c M Dimension of the multivariate normal c MU Mean of multivariate normal, a M x 1 vector c W Covariance matrix of multivariate normal, a M x M c array c WI Inverse of W c C Lower triangular Cholesky factor of W, a M x M array c A Lower bound of rectangle, a M x 1 vector [When the c lower bound is -infinity, set A = (-1.0E10)*ONES(M,1)] c B Upper bound of rectangle, a M x 1 vector [When the c upper bound is +infinity, set B = (+1.0E10)*ONES(M,1)] c R Number of repetitions c U Random variates, a M x R array c PARM Parameters and constants for the simulation c c The simulators all return {P,HU,HC}, where P is the scalar rectangle c probability, HU is the (M+1) x (M+1) array of unconditional partial c moments (6), and HC is the (M+1) x (M+1) array of conditional moments c (7). Parts of the output not provided by a simulator are set to c -999. c c Interpretations of input parameters c ----------------------------------- c ARSR,ARSE c --------- c These procs assume that PARM[.,1] contains seeds for the uniform c random number generator RNDUS called by the routines; there is one c seed for each repetition R. PARM[R+1,1] is assumed to contain a c censoring limit for rejection. c ARSR and ARSE use U = UUBIG; (uniforms) c c AUS c --- c The proc assumes that PARM(.,1) contains a vector of R seeds for the c normal random number generator RNDNS, PARM(1,2) contains the number c of initial GHK simulators of P used in the simulation of 1/P, c PARM(2,2) contains the maximum number of trials using the crude c frequency simulator that are made before the process is censored. If c PARM(2,2) is made large (say 1.e50), then this is computationally the c same as the Sequential Unbiased Simulator (SUS). PARM(3,2) contains c a seed for RNDUS. c AUS uses U = UUBIG[.,1:R]; (uniforms) c c CFS c --- c This proc does not use PARM. c CFS uses U = UNBIG[.,1:R] (standard normals) c c DCS c --- c The routine assumes that PARM(1,1) = DET(W). c DCS uses U = USBIG[.,1:R]; (random normal grid on unit sphere) c c GHK c --- c This proc does not use PARM. c GHK uses U = UUBIG[.,1:R]; (uniforms) c c GSS c --- c It assumes that PARM[1,1] gives the number of rounds to be used in c the construction of the random draws. c GSS uses U = UUBIG[.,1:R]; (uniforms) c c NISE,NIST c --------- c These procs do not use PARM. c NISE and NIST use U = UUBIG[.,1:R]; (uniforms) c c PCF c --- c The routine assumes that PARM(1,1) = DET(W). c PCF uses U = USSBIG[.,1:R]; (random normal grid on unit sphere and c orthant) c c PKF c --- c This proc does not use W. c It assumes that PARM[1,1] contains a window width parameter. c KFS uses U = UNBIG[.,1:R]; (standard normals) c c SDS c --- c The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such that c W-lI is positive definite. c SDS uses U = UNBIG[.,1:R]; (standard normals) c Non-standard feature: c -------------------- c SDS assumes that C is a Cholesky factor of W-lI, not of W. c c SUS c --- c This procedure is the same the Asymptotically Unbiased Simulator c (AUS) with PARM(2,2) made large (1.e50). c The proc assumes that PARM(.,1) contains a vector of R seeds for the c normal random number generator RNDNS, PARM(1,2) contains the number c of initial GHK simulators of P used in the simulation of 1/P, c PARM(2,2) contains the maximum number of trials using the crude c frequency simulator that are made before the process is censored. c This is made very large. c PARM(3,2) contains a seed for RNDUS. c SUS uses U = UUBIG[.,1:R]; (uniforms) c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c Test Harness c ------------ c The following program provides a simple harness for testing the c procs above. A one-factor model is used to generate probabilities, c with gaussian quadrature used for evaluation. The accuracy of this c quadrature can be checked against symmetric cases where the c probabilities are known. c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) double precision nrmpdf character*2 resfil character*130 results character*12 cfilCFS,cfilNISE,cfilNIST,cfilKFS,cfilSDS, / cfilGHK,cfilPCF,cfilDCS,cfilARSE,cfilARSR, / cfilGSS,cfilAUS,cfilSUS c include 'hdim.inc' c parameter( * larse = mmaxim*irmaxim+mmaxim*mmaxim+8*mmaxim, * larsr = mmaxim*irmaxim+irmaxim+2*mmaxim*mmaxim+3*mmaxim, * laus = 4*mmaxim*irmaxim+2*irmaxim+2*mmaxim*mmaxim-2*mmaxim+3, * lcfs = 3*mmaxim*irmaxim+irmaxim+mmaxim*mmaxim+mmaxim, * ldcs = 3*mmaxim*irmaxim+3*irmaxim+mmaxim*mmaxim+3*mmaxim+2, * lghk = 3*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+2*mmaxim, * lgss = mmaxim*irmaxim+3*mmaxim*mmaxim+3*mmaxim-2, * lnise = 6*mmaxim*irmaxim+9*mmaxim+2*irmaxim+mmaxim*mmaxim, * lnist = 5*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+3*mmaxim, * lpcf = 4*mmaxim*irmaxim+11*irmaxim+mmaxim*mmaxim+4*mmaxim+2, * lpkf = 5*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+mmaxim, * lsds = 10*mmaxim*irmaxim+2*mmaxim*mmaxim+2*mmaxim, * lsus = 4*mmaxim*irmaxim+2*irmaxim+2*mmaxim*mmaxim-2*mmaxim+3) parameter( * lcom=lsds) cLaheyBug * lcom = max(larse,larsr,laus ,lcfs ,ldcs ,lghk ,lgss , cLaheyBug * lnise,lnist,lpcf ,lpkf ,lsds ,lsus)) common // dummy(lcom) common /times/ hseclim c dimension hu(mmaxim+1,mmaxim+1),hc(mmaxim+1,mmaxim+1), / rsim(13),parm(irmaxim+1,2), / xmu(mmaxim),a(mmaxim),b(mmaxim),xl(mmaxim),d(mmaxim), / s(mmaxim),p(1),pvec(mmaxim), / w(mmaxim,mmaxim),wi(mmaxim,mmaxim),chw(mmaxim,mmaxim), / tempmxm(mmaxim,mmaxim),tempmx1(mmaxim), / temp1xm(1,mmaxim), / t(mmaxim),tmm1(mmaxim-1), / uuBIG(mmaxim,irmaxim), / us(mmaxim,mmaxim),uss(mmaxim,2*mmaxim+4), / unBIG(mmaxim,irmaxim),unBIG2(mmaxim,irmaxim), / unzsum(irmaxim),unzsumt(1,irmaxim), / usBIG(mmaxim,irmaxim),ussBIG(mmaxim,irmaxim), / aw(mmaxim,mmaxim), / p40(40),t40(40),tt40(40), / pp0(16),dm0(16),xlm0(16),dl0(16),xll0(16),rep0(16), / el0(16),aa(16,7),res(7), / qw(20,2),qw2(40,2),el(16),rep(16),pp(16),dm(16), / xlm(16),dl(16),xll(16),vec40(40),vec40zl(40) dimension index(mmaxim) data qw/ * 0.048613160848618d0,0.146006911993027d0,0.243907555937767d0, * 0.342668563127518d0,0.442667961120606d0,0.544320702552795d0, * 0.648093938827515d0,0.754526972770691d0,0.864258110523224d0, * 0.978062629699707d0,1.096907615661621d0,1.222035169601440d0, * 1.355093598365784d0,1.498357176780701d0,1.655115008354187d0, * 1.830419778823853d0,2.032696485519409d0,2.277773618698120d0, * 2.601793527603149d0,3.127621173858643d0, * 0.038752973079681d0,0.038519907742739d0,0.038055181503296d0, * 0.037361584603786d0,0.036443289369345d0,0.035305824130774d0, * 0.033956021070480d0,0.032402005046606d0,0.030653120949864d0, * 0.028719885274768d0,0.026613922789693d0,0.024347903206944d0, * 0.021935453638434d0,0.019391084089875d0,0.016730098053813d0, * 0.013968503102660d0,0.011122924275696d0,0.008210529573262d0, * 0.005249142181128d0,0.002260638633743d0/ c write(6,'(a/)') *' *-----------------------------------------------------------*', *' | Simulation Algorithms |', *' | for Orthant Normal Probabilities and Derivatives: |', *' | A Harness Test Program Run |', ccc *' | version 1.0, March 22, 1992 |', *' | version 1.1, July 11, 1992 |', *' | Written by Vassilis Hajivassiliou |', *' | with the aid of Yoosoon Chang |', *' | This Program is Part of the Paper: |', *' | Simulation of Normal Orthant Probabilities: |', *' | Methods and Programs |', *' | by Hajivassiliou, McFadden, and Ruud |', ccc *' | February 1992 |', *' | July 1992 |', *' *-----------------------------------------------------------*' call machine(results) write(6,'(/a)') *' --------------------------' write(6,'(a)') *' Computer Used for this Run: '//results(1:ltrim(results)) write(6,'(a/)') *' --------------------------' call when(6) c c The following allows for the first line in the input file c required by the gauss program, giving the number of numeric c parameters in the input file. (See the GAUSS program.) print *, 'Hit when ready:' read(*,'(a)') results print *, ' ' c 5 print *, 'enter input dimension m : ' read(*,*) temp m=temp print *,'***',m if (m.lt.2) then print *, 'm should be >= 2 ' goto 5 elseif (m.gt.mmaxim) then print *, 'm cannot exceed ',mmaxim print *, 'which was the limit at compilation time' goto 5 endif print *, 'enter seed for random number generator: ' read(*,*) temp iseed1=temp print *,'***',iseed1 iseed = iseed1 c print *,'Enter total number of Monte Carlo repetitions: ' read(*,*) temp mcrepet=temp print *,'***',mcrepet c write(6,1) 1 format(' enter 1 to calculate ENDOGENOUSLY number of', / ' simulations for each method,' / / ' enter 0 to specify number for each method yourself: ') read(*,*) temp iendog=temp print *,'***',iendog c if ( iendog .eq. 1 ) then 10 print *, 'Num. of simul. for CFS ', / '(will calculate endogenously for other methods): ' read(*,*) temp irsimCFS=temp print *,'***',irsimCFS if (irsimCFS .lt. m+2 ) then print *,'Num. of simul. for CFS should be >= m+2 ' goto 10 endif print *, 'Num. of simul. for endogenous timing of methods ', / '(irsimTRY): ' read(*,*) temp irsimTRY=temp print *,'***',irsimTRY if ( irsimTRY .gt. irsimCFS ) then print *, 'ERROR!! irsimTRY cannot exceed irsimCFS' goto 10 endif irsimBIG = irsimCFS rsimBIG = irsimBIG else 15 print *, 'Number of repetitions for CFS: ' read(*,*) temp irsimCFS=temp print *,'***',irsimCFS rsim(1) = irsimCFS print *, 'Number of repetitions for NIS_exponential: ' read(*,*) temp irsimNSE=temp print *,'***',irsimNSE rsim(2) = irsimNSE print *, 'Number of repetitions for NIS_truncated_normal: ' read(*,*) temp irsimNST=temp print *,'***',irsimNST rsim(3) = irsimNST print *, 'Number of repetitions for KFS: ' read(*,*) temp irsimKFS=temp print *,'***',irsimKFS rsim(4) = irsimKFS print *, 'Number of repetitions for SDS: ' read(*,*) temp irsimSDS=temp print *,'***',irsimSDS rsim(5) = irsimSDS print *, 'Number of repetitions for GHK: ' read(*,*) temp irsimGHK=temp print *,'***',irsimGHK rsim(6) = irsimGHK print *, 'Number of repetitions for PCF: ' read(*,*) temp irsimPCF=temp print *,'***',irsimPCF rsim(7) = irsimPCF print *, 'Number of repetitions for DCS: ' read(*,*) temp irsimDCS=temp print *,'***',irsimDCS rsim(8) = irsimDCS print *, 'Number of repetitions for ARS_exponential: ' read(*,*) temp irsimARE=temp print *,'***',irsimARE rsim(9) = irsimARE print *, 'Number of repetitions for ', / 'ARS_recursive_truncated_normal: ' read(*,*) temp irsimARR=temp print *,'***',irsimARR rsim(10) = irsimARR print *, 'Number of repetitions for GSS: ' read(*,*) temp irsimGSS=temp print *,'***',irsimGSS rsim(11) = irsimGSS print *, 'Number of repetitions for AUS: ' read(*,*) temp irsimAUS=temp print *,'***',irsimAUS rsim(12) = irsimAUS print *, 'Number of repetitions for SUS: ' read(*,*) temp irsimSUS=temp print *,'***',irsimSUS rsim(13) = irsimSUS print *, ' ' irsimBIG = fmax(rsim,13) if (irsimBIG .lt. m+2 ) then print *, 'For at least one method, num. of repetition ', * 'should be >= m+2 ' goto 15 elseif (irsimBIG .gt. irmaxim) then print *, * '*********************************************', * '* For at least one method, you have requested', * '* TOO LARGE a number of simulations (',irsimBIG,')', * '* which exceeds the LIMIT at compilation (',irmaxim,')', * '*********************************************' goto 15 endif rsimBIG = irsimBIG endif print *, 'Time limit (in seconds) for single run of any method: ' read(*,*) seconds print '(1x,a,f10.3)','***',seconds hseclim = seconds * 100.d0 c write(6,2) 2 format(' Do you want a special case with p analytic? ' / / 3x,'enter 0 = No' / / 3x,' 1 = Factor loading on m only' / / 3x,' 2 = Orthant prob. with full symmetry'/ / 3x,' 3 = Arbitrary omega') read(*,*) temp icase=temp print *,'***',icase c if ( icase .eq. 2 ) then do 40 i = 1,m xmu(i) = 0.0d0 b(i) = 0.0d0 a(i) = -25.0d0 xl(i) = 1.0d0 d(i) = 1.0d0 40 continue else write(6,3) 3 format(' Do you want automatic generation of data? ' / / ' Enter 1=Yes, 0=No') read(*,*) temp iautogen=temp print *,'***',iautogen if ( iautogen .eq. 1 ) then call mrnorm(xmu,mmaxim,m,1,iseed) call fillm(b,mmaxim,1,0.0d0) do 20 i = 1,m a(i) = -2.0d0 - dabs(xmu(i)) 20 continue if ( icase .eq. 0 ) then call mrnorm(xl,mmaxim,m,1,iseed) call mdscal(xl,mmaxim,xl,mmaxim,m,1,-3.0d0) else if (icase .eq. 1) then call fillm(xl,mmaxim,1,0.0d0) xl(m) = -3.0d0 * ranorm(1.0d0,0.0d0,iseed) else if (icase .eq. 3) then call fillm(xl,mmaxim,1,-999.0d0) endif if (icase .eq. 3) then call fillm(d,mmaxim,1,-999.0d0) else call mrnorm(d,mmaxim,m,1,iseed) do 30 i = 1,m d(i) = 0.5d0 + 0.3d0*dabs(d(i)) 30 continue endif else print *, 'Enter means: ' read(*,*) (xmu(i), i=1,m) print *, 'Enter lower bounds: (A)' read(*,*) (a(i), i=1,m) print *, 'Enter upper bounds: (B)' read(*,*) (b(i), i=1,m) if (icase.eq.0 .or. icase.eq.1) then print *, 'Enter indep. std. dev.: ' read(*,*) (d(i), i=1,m) endif call hadcomp(b,mmaxim,a,mmaxim,tempmx1,mmaxim,m,1,'gt') do 70 i = 1,m b(i) = a(i) + 1.0d0 + (b(i)-1.0d0-a(i))*tempmx1(i) 70 continue if (icase .ne. 3) then call fillm(tempmx1,mmaxim,1,0.0d0) call hadcomp(d,mmaxim,tempmx1,mmaxim, * tempmx1,mmaxim,m,1,'gt') do 80 i = 1,m d(i) = 1.0d0 + (d(i)-1.0d0)*tempmx1(i) 80 continue else call fillm(d,mmaxim,1,-999.0d0) endif if ( icase .eq. 0 ) then print *, 'Enter Factor Loadings: ' read(*,*) (xl(i), i=1,m ) else if (icase .eq. 1) then print *, 'Enter Factor Loading Component m: ' read(*,*) temp m=temp call fillm(xl,mmaxim,1,1.0d0) xl(m) = m else if (icase .eq. 3) then call fillm(xl,mmaxim,1,-999.0d0) print *,'Enter lower triange of OMEGA, ', / 'one row at a time:' do 81 i=1,m read(*,*) (w(i,j), j=1,i) print *, (w(i,j), j=1,i) 81 continue endif endif endif c print *, 'Enter 1 to print FULL output on the screen' print *, ' 0 to print SUMMARY output' read(*,*) temp ifullout=temp print *,'***',ifullout c c The following allows for the line in the input file c required by the gauss program, giving the number of string c parameters in the input file. (See the GAUSS program.) print *, 'Hit when ready:' read(*,'(a)') results print *, ' ' c c print *, ' ' print *, 'Basic Name of file to save results in: ', / '(MUST be 2-letters long)' results(3:3)=' ' read(*,'(a)') results if(results(3:3).eq.' ') then resfil=results(1:2) print *, ' '//resfil print *, ' ' else print *,'Error!!! Name must be EXACTLY 2-letters long' goto 999 endif c--assign a name to output file to each methods. cfilCFS = resfil//'CFS.asd' cfilNISE = resfil//'NISE.asd' cfilNIST = resfil//'NIST.asd' cfilKFS = resfil//'KFS.asd' cfilSDS = resfil//'SDS.asd' cfilGHK = resfil//'GHK.asd' cfilPCF = resfil//'PCF.asd' cfilDCS = resfil//'DCS.asd' cfilARSE = resfil//'ARSE.asd' cfilARSR = resfil//'ARSR.asd' cfilGSS = resfil//'GSS.asd' cfilAUS = resfil//'AUS.asd' cfilSUS = resfil//'SUS.asd' 700 continue call fsize(cfilSUS,itemp) if(mcrepet.gt.0) then if(itemp.gt.0) then c *.asd files must NOT exist already print '(1x,a/1x,a)', * 'Error!!!', * 'When MCrepet is Positive, *.asd files must NOT exist.' goto 999 else imcrepet=1 endif elseif(mcrepet.lt.0) then if(itemp.le.0) then c *.asd files MUST exist already print '(1x,a/1x,a)', * 'Error!!!', * 'When MCrepet is Negative, *.asd files MUST exist.' goto 999 else c must find actual completed repetition in files open(10,file=cfilSUS, * access='direct',form='formatted',recl=130) irec=1 777 continue ccc read(10,'(a)',rec=irec,end=778) results read(10,'(a)',rec=irec,err=778) results irec=irec+1 goto 777 778 continue close(10) imcrepet=irec-2 print *,'******************************************' print *,'*** Will proceed from existing results ***' print *,'*** Completed repetitions: ',imcrepet print *,'******************************************' imcrepet=imcrepet+1 mcrepet=-mcrepet endif else print *,'The number of MC repetitions cannot be 0' print *,'Will carry out 1 MC repetition' mcrepet=1 goto 700 endif write(6,84) 84 format(' The inputs for this exercise are:',/1x/ * 5x,'MEANS',10x,'A',14x,'B',9x,'SD',4x,'FACTOR LOADINGS'/ * ' -------------------------------------------------------', * ' ------------') do 85 i = 1,m write(6,'(1x,5g13.5)') xmu(i),a(i),b(i),d(i),xl(i) 85 continue print *,' ' print *,'Initial Seed: ',iseed1 do 90 i = 1,m s(i) = dsqrt( d(i)**2 + xl(i)**2 ) pvec(i) = / pheev((b(i)-xmu(i))/s(i)) - pheev((a(i)-xmu(i))/s(i)) 90 continue index(1)=1 call srtrnk(pvec,m,index) call matcopy(tempmx1,mmaxim,d ,mmaxim,1,m,1,1) do 701 i=1,m d(i)=tempmx1(index(i)) 701 continue call matcopy(tempmx1,mmaxim,xl ,mmaxim,1,m,1,1) do 702 i=1,m xl(i)=tempmx1(index(i)) 702 continue call matcopy(tempmx1,mmaxim,a ,mmaxim,1,m,1,1) do 703 i=1,m a(i)=tempmx1(index(i)) 703 continue call matcopy(tempmx1,mmaxim,b ,mmaxim,1,m,1,1) do 704 i=1,m b(i)=tempmx1(index(i)) 704 continue call matcopy(tempmx1,mmaxim,xmu,mmaxim,1,m,1,1) do 705 i=1,m xmu(i)=tempmx1(index(i)) 705 continue c write(6,706) 706 format(1x/ * ' After sorting to put low probability in first component:'/2x/ * 5x,'Means',10x,'A',12x,'B',11x,'Sd',4x,'Factor loadings'/ * ' --------------------------------', * '---------------------------------') do 100 i = 1,m write(6,'(1x,5g13.5)') xmu(i),a(i),b(i),d(i),xl(i) 100 continue if (icase .ne. 3) then do 110 i = 1,m tempmx1(i) = d(i)**2 110 continue call fillm(w,mmaxim,m,0.0d0) call diagrv(w,mmaxim,m,m,tempmx1,mmaxim) call transps(xl,mmaxim,m,1,temp1xm,1) call matpv(xl,mmaxim,temp1xm,1,tempmxm,mmaxim,m,1,m) call hadoper(w,mmaxim,tempmxm,mmaxim,w,mmaxim,m,m,1) else do 111 i=2,m do 112 j=1,i-1 w(m-i+1,m-j+1) = w(i,j) 112 continue 111 continue endif c c calculate chw=CHOL(W) call matcopy(tempmxm,mmaxim,w,mmaxim,1,m,1,m) call chol(tempmxm,mmaxim,m,chw,mmaxim) c c calculate wi=INV(W) call matcopy(wi,mmaxim,w,mmaxim,1,m,1,m) call invv(wi,tempmxm,mmaxim,mmaxim,tempmx1,m, * npr,nzr,nnr,condnum,*888) c c calculate detw=DET(W) call prodc(tempmx1,mmaxim,m,1,t(1),1) detw=t(1) c print *,' ' print *,'Implied Variance-Covariance Matrix' print *,'----------------------------------' call mprint(w,mmaxim,m,m,'omega',6,6) print '(/1x,a/1x,a,e15.7,a,e15.7)', * 'Var-Cov matrix inverted and its Cholesky factor '// * 'calculated.','Condition Number = ', condnum, * ' Determinant = ',detw c do 120 i = 1,20 qw2(i,1) = qw(i,1) qw2(i,2) = qw(i,2) qw2(i+20,1) = -qw(i,1) qw2(i+20,2) = qw(i,2) 120 continue c ccc---initialize to "Not available" c gresampl = -999.d0 c do 125 i = 1,16 pp(i) = -999.0d0 dm(i) = -999.0d0 xlm(i) = -999.0d0 dl(i) = -999.0d0 xll(i) = -999.0d0 125 continue c do 130 i = 1,2 el(i) = 1.0d0 rep(i) = 1.d0 130 continue c--if omega is arbitrary, skip the quadratures if (icase.eq.3) goto 1000 c DRAWRAND c ccc======================= quadrature 1 ============================ccc c mmaximm1 = mmaxim - 1 mp1 = m + 1 mm1 = m - 1 c if (icase.eq.1) then call diag(w,mmaxim,m,s,mmaxim) do 140 i = 1,m s(i) = dsqrt(s(i)) xhs = pheev( (b(i)-xmu(i)) / s(i) ) rhs = pheev( (a(i)-xmu(i)) / s(i) ) t(i) = xhs - rhs 140 continue call prodc(t,mmaxim,m,1,tempmx1(1),1) pp(1) = tempmx1(1) bm = (b(m)-xmu(m)) / s(m) am = (a(m)-xmu(m)) / s(m) q = nrmpdf(bm,0.0d0,1.0d0) - nrmpdf(am,0.0d0,1.0d0) qq = bm*nrmpdf(bm,0.0d0,1.0d0) - am*nrmpdf(am,0.0d0,1.0d0) call matcopy(tmm1,mmaximm1,t,mmaxim,1,mm1,1,1) call prodc(tmm1,mmaximm1,mm1,1,dlp,1) dl(1) = -dlp dm(1) = dl(1)*q / s(m) xlm(1) = dm(1) / pp(1) dl(1) = dl(1)*qq*xl(m) / s(m)**3 xll(1) = dl(1) / pp(1) endif c if (icase.eq.2) then pp(1) = 1.0d0 / dble(float(mp1)) do 150 i = 1,40 cdf = pheev(qw2(i,1)**mm1) pdf = nrmpdf(qw2(i,1),0.0d0,1.0d0) vec40(i) = cdf*pdf*qw2(i,2) vec40zl(i) = cdf*pdf*qw2(i,1)*qw2(i,2) 150 continue call sumc(vec40,40,40,1,dmp,1) dm(1) = -dmp xlm(1) = dm(1) / pp(1) call sumc(vec40zl,40,40,1,dlp,1) dl(1) = dlp xll(1) = dl(1) / pp(1) endif c ccc======================== quadrature 2 =========================ccc c if (m.lt.3) then call exaf(mmaxim,m,a,b,xmu,w,xl,pzscal,tzscal,ttzscal) pp(2) = pzscal dm(2) = tzscal dl(2) = ttzscal xlm(2) = tzscal / pzscal xll(2) = ttzscal / pzscal else pp(2) = -999.0d0 dm(2) = -999.0d0 dl(2) = -999.0d0 xlm(2) = -999.0d0 xll(2) = -999.0d0 endif c 1000 print *, ' ' c So that for a given iseed1, the same iseed will be used c at a given imcrepet: iseed=iseednew(iseed1,imcrepet) ccc print *, 'DRAWRAND' ccc print *, 'Current seed: ', iseed curseed = iseed call mrunif(uuBIG,mmaxim,m,irsimBIG,iseed) call rndbase(mmaxim,m,iseed,us,t) it = idnint(dble( float(irsimBIG-m-1) / float(m-1) )) call antith(mmaxim,m,irsimBIG,it,us,uss) c-----note: us[1=m,1=m], uss[1=m,1=2*m+4] do 161 j = 1, 2*m+4 do 162 i = 1,m uss(i,j) = -dabs(uss(i,j)) 162 continue 161 continue call mrnorm(unBIG,mmaxim,m,irsimBIG,iseed) call hadoper(unBIG,mmaxim,unBIG,mmaxim,unBIG2,mmaxim,m,irsimBIG,3) call sumc(unBIG2,mmaxim,m,irsimBIG,unzsum,irmaxim) call transps(unzsum,irmaxim,irsimBIG,1,unzsumt,1) do 171 j = 1,irsimBIG do 172 i = 1,m usBIG(i,j) = unBIG(i,j) / dsqrt(unzsumt(1,j)) ussBIG(i,j) = -dabs(usBIG(i,j)) if (ussBIG(i,j) .eq. 0.0d0) then ussBIG(i,j) = ussBIG(i,j) + 1.d-10 endif 172 continue 171 continue c vah if(imcrepet.gt.1) goto 9995 c ENDTIM c ccc=========================== Exact ==============================ccc c ccc obtain exact results only for first mcrepetitions ccc Use as "true" GHK results with 10*irsimBIG repetitions. if (icase.eq.3) then do 174 i=1,2 pp(i)=-999.0d0 dm(i)=-999.0d0 xlm(i)=-999.0d0 dl(i)=-999.0d0 xll(i)=-999.0d0 rep(i)=-999.0d0 el(i)=-999.0d0 174 continue c print *,'*** As TRUE, will use GHK results with 10*irsimBIG ', / 'simulations ***' print *,' ' c ccc------GHK c c vah e1 = hsec1990() e1top = hsec1990() ir = irsimBIG spp = 0 sdm = 0 sxlm = 0 sdl = 0 sxll = 0 do 173 ipart=1,10 print *,'Part ',ipart,' of 10' call mrunif(uuBIG,mmaxim,m,ir,iseed) print *,'back from mrunif' call ghk(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) print *,'back from ghk: p=',p call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) print *,'back from sim: pploc:',pploc spp = spp + pploc sdm = sdm + dmloc sxlm = sxlm + xlmloc sdl = sdl + dlloc sxll = sxll + xllloc 173 continue pp(3) = spp/10.d0 dm(3) = sdm/10.d0 xlm(3) = sxlm/10.d0 dl(3) = sdl/10.d0 xll(3) = sxll/10.d0 rep(3) = 1.0d0 el(3) = hsec1990() - e1top goto 7000 c HAVEXACT endif c e1 = hsec1990() do 177 i = 1,40 vec40(i) = qw2(i,1) vec40zl(i) = qw2(i,2) 177 continue c ccc---note: calling ffact here is different from that in gauss ccc program.(the gauss call does not include xmu as an argument) c call ffact(mmaxim,m,vec40,a,b,xmu,d,xl,p40,t40,tt40) call matptv(p40,40,vec40zl,40,pp(3),1,40,1,1) call matptv(t40,40,vec40zl,40,dm(3),1,40,1,1) call matptv(tt40,40,vec40zl,40,dl(3),1,40,1,1) if(pp(3).gt.0.d0) then xlm(3) = dm(3) / pp(3) xll(3) = dl(3) / pp(3) else xlm(3)=-999.d0 xll(3)=-999.d0 endif if (dl(1) .eq. 0.0d0) then dl(1) = dl(3) endif if (xll(1) .eq. 0.0d0) then xll(1) = xll(3) endif c vah 22 sep 92: added to guard against meanigless derivative zeros c returned from ffact: if(dm(3).eq.0.d0.and.dl(3).eq.0.d0) then dm(3)=-999.d0 xlm(3)=-999.d0 dl(3)=-999.d0 xll(3)=-999.d0 endif write(6,181) 181 format(1x/'Quadrature1, Quadrature2, and Exact Values'/ * ' ---------------------------------------', * '----------------------------------------------') write(6,186) 186 format(20x,'p',8x,'dp/dmu(m)',4x,'dln(p)/dmu(m)',5x,'dp/dl(m)', * 5x,'dln(p)/dl(m)'/ * ' ---------------------------------------', * '----------------------------------------------') write(6,182) pp(1),dm(1),xlm(1),dl(1),xll(1) 182 format(' Quadrature1',(5g15.5)) write(6,183) pp(2),dm(2),xlm(2),dl(2),xll(2) 183 format(' Quadrature2',(5g15.5)) write(6,184) pp(3),dm(3),xlm(3),dl(3),xll(3) 184 format(' Exact',6x,(5g15.5)) write(6,187) 187 format(' ---------------------------------------', * '----------------------------------------------'/) rep(3) = 1.0d0 e2 = hsec1990() ec = e2 - e1 el(3) = ec e1 = e2 C vah HAVEXACT 7000 continue ccc---save exact results c call matcopy(pp0,16,pp,16,1,16,1,1) call matcopy(dm0,16,dm,16,1,16,1,1) call matcopy(xlm0,16,xlm,16,1,16,1,1) call matcopy(dl0,16,dl,16,1,16,1,1) call matcopy(xll0,16,xll,16,1,16,1,1) call matcopy(el0,16,el,16,1,16,1,1) call matcopy(rep0,16,rep,16,1,16,1,1) c ccc---write "true" values, using exact/quadrature average c do 185 i = 1,16 aa(i,1) = pp(i) aa(i,2) = dm(i) aa(i,3) = xlm(i) aa(i,4) = dl(i) aa(i,5) = xll(i) aa(i,6) = rep(i) aa(i,7) = el(i) 185 continue c ccc---check each row of aa, if any element is -999.0d0, then delete c the row. c call packr(aa,16,3,7,-999.0d0,icount) if(icount.gt.0) then call meanc(aa,16,icount,7,res,7) else call fillv(res,7,-999.d0) endif write(6,188) (res(i),i=1,7) 188 format(' Will use the following as TRUE values:'/ * ' ---------------------------------------', * '----------------------------------------', * '----------------------'/ * 11x,'p',8x,'dp/dmu(m)',4x,'dln(p)/dmu(m)',5x,'dp/dl(m)', * 5x,'dln(p)/dl(m)',5x,'nrepet',5x,'time(sec)'/ * ' ---------------------------------------', * '----------------------------------------', * '----------------------'/ * (7g15.5)/ * ' ---------------------------------------', * '----------------------------------------', * '----------------------'/) c ccc put exact values to the first results of each output files c write(results,'(8(1x,g15.9))') (res(i),i=1,7),curseed irec=1 call reswrite(cfilCFS ,results,irec) call reswrite(cfilNISE,results,irec) call reswrite(cfilNIST,results,irec) call reswrite(cfilKFS ,results,irec) call reswrite(cfilSDS ,results,irec) call reswrite(cfilGHK ,results,irec) call reswrite(cfilPCF ,results,irec) call reswrite(cfilDCS ,results,irec) call reswrite(cfilARSE,results,irec) call reswrite(cfilARSR,results,irec) call reswrite(cfilGSS ,results,irec) call reswrite(cfilAUS ,results,irec) call reswrite(cfilSUS ,results,irec) if (iendog .eq. 0) goto 9995 c ENDTIM c c*********************************************************************** c* a simple timing test to set number of simulations endogenously * c*********************************************************************** c print *, 'BEGINTIM' print *,' ' c ccc---CFS e1 = hsec1990() ir = irsimCFS call cfs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc) timeCFS = dmax1(1.0d0, hsec1990()-e1) print *, 'Time for ',ir,' simulations using CFS: ',timeCFS print *, ' ' c write(6,3000) 3000 format(' ---------------------------------------', / '-------------------------------------------------'/ / 'Method',24x,'irsimBIG',7x,'# within',7x,'# equivalent', / 5x,'Number of'/ / 45x,'time limit',11x,'to CFS',5x,'repetitions'/ / '---------------------------------------', / '-------------------------------------------------') c ccc---NISE e1 = hsec1990() ir = irsimTRY call nise(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimNSE = dmax1(10.0d0,out3) write(6,3001) irsimBIG,out1,out2,irsimNSE 3001 format(' NIS_exponential:',15x,i5,2g20.5,i9) c ccc---NIST e1 = hsec1990() ir = irsimTRY call nist(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimNST = dmax1(10.0d0,out3) write(6,3002) irsimBIG,out1,out2,irsimNST 3002 format(' NIS_truncated_normal:',10x,i5,2g20.5,i9) c ccc---KFS e1 = hsec1990() ir = irsimTRY do 200 i = 1, m vec40(i) = b(i) - a(i) 200 continue do 210 i = m+1, 2*m vec40(i) = dsqrt(d(i-m)**2 + xl(i-m)**2) 210 continue val = 0.1d0 * fmin(vec40,2*m) call fillm(parm,irmaxim+1,2,val) call pkfs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimKFS = dmax1(10.0d0,out3) write(6,3003) irsimBIG,out1,out2,irsimKFS 3003 format(' KFS:',27x,i5,2g20.5,i9) c ccc---SDS e1 = hsec1990() ir = irsimTRY call matcopy(aw,mmaxim,w,mmaxim,1,m,1,m) c c-----invv returns eigen values of w in tempmx1 in descending order, c so, sorts is used to make it in ascending order. c call invv(aw,tempmxm,mmaxim,mmaxim,tempmx1,m, * npr,nzr,nnr,rcond,*888) call sorts(tempmx1,m) val = (0.99d0 * tempmx1(1))**2 call fillm(parm,irmaxim+1,2,val) do 220 j = 1,m do 221 i = 1,m if (i.ne.j) then tempmxm(i,j) = w(i,j) else tempmxm(i,j) = w(i,j) - val endif 221 continue 220 continue call chol(tempmxm,mmaxim,m,tempmxm,mmaxim) call sds(mmaxim,irmaxim, / m,xmu,w,wi,tempmxm,a,b,ir,unBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimSDS = dmax1(10.0d0,out3) write(6,3004) irsimBIG,out1,out2,irsimSDS 3004 format(' SDS:',27x,i5,2g20.5,i9) c ccc---GHK e1 = hsec1990() ir = irsimTRY call ghk(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimGHK = dmax1(10.0d0,out3) write(6,3005) irsimBIG,out1,out2,irsimGHK 3005 format(' GHK:',27x,i5,2g20.5,i9) c ccc---PCF e1 = hsec1990() ir = irsimTRY parm(1,1) = detw call pcf(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,ussBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimPCF = dmax1(10.0d0,out3) write(6,3006) irsimBIG,out1,out2,irsimPCF 3006 format(' PCF:',27x,i5,2g20.5,i9) c ccc---DCS e1 = hsec1990() ir = irsimTRY parm(1,1) = detw call dcs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,usBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimDCS = dmax1(10.0d0,out3) write(6,3007) irsimBIG,out1,out2,irsimDCS 3007 format(' DCS:',27x,i5,2g20.5,i9) c ccc---ARSE e1 = hsec1990() ir = irsimTRY issd = idnint(uuBIG(1,1)*1.0d5) do 230 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(issd) ) 230 continue parm(ir+1,1) = ir call arse(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimARE = dmax1(10.0d0,out3) write(6,3008) irsimBIG,out1,out2,irsimARE,dummy(1) 3008 format(' ARS_exponential:',15x,i5,2g20.5,i9, *' % not censored:',f10.7) c ccc---ARSR e1 = hsec1990() ir = irsimTRY issd = idnint(uuBIG(1,1)*1.0d5) do 231 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(issd) ) 231 continue parm(ir+1,1) = ir call arsr(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimARR = dmax1(10.0d0,out3) write(6,3009) irsimBIG,out1,out2,irsimARR,dummy(1) 3009 format(' ARR_recursive_truncated_normal:',i5,2g20.5,i9, *' % not censored:',f10.7) c ccc---GSS e1 = hsec1990() ir = irsimTRY gresampl = idnint( 5.0d0 + dble(float(ir)/20.0d0) ) cvah added so as to make sure at least 10 gibbs resamplings are used: if(val.lt.10.d0) gresampl=10.d0 call fillm(parm,irmaxim+1,2,gresampl) call gss(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimGSS = dmax1(10.0d0,out3) write(6,3010) irsimBIG,out1,out2,irsimGSS 3010 format(' GSS:',27x,i5,2g20.5,i9) c ccc---AUS e1 = hsec1990() ir = irsimTRY do 240 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(iseed) ) parm(i,2) = 1.0d0 240 continue parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) ) parm(2,2) = idnint( 10.0d0 + dble(rsimBIG/10.0d0) ) parm(3,2) = idnint( 1.0d6*uuBIG(2,1) ) call aus(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimAUS = dmax1(10.0d0,out3) write(6,3011) irsimBIG,out1,out2,irsimAUS,dummy(1) 3011 format(' AUS:',27x,i5,2g20.5,i9,' Max. # of CFS draws:',f10.1) c ccc---SUS e1 = hsec1990() ir = irsimTRY do 241 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(iseed) ) parm(i,2) = 1.0d0 241 continue parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) ) parm(3,2) = idnint( 1.0d6*uuBIG(2,1) ) call sus(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) etime = dmax1(1.0d0, hsec1990()-e1) out1 = hseclim / (etime/dble(float(ir))) out2 = dble(float(ir))*timecfs / etime out3 = dmin1(rsimBIG,out1,out2) irsimSUS = dmax1(10.0d0,out3) write(6,3012) irsimBIG,out1,out2,irsimSUS,dummy(1) 3012 format(' SUS:',27x,i5,2g20.5,i9,' Max. # of CFS draws:',f10.1/ * ' ---------------------------------------', * '------------------------------------------------'/) c c c write(6,3013) 3013 format(' Number of Repetitions'/ * ' ------------------------------------------', * '-------------------------------------------------------') write(6,3014) irsimCFS,irsimNSE,irsimNST,irsimKFS,irsimSDS, * irsimGHK,irsimPCF,irsimDCS,irsimARE,irsimARR, * irsimGSS,irsimAUS,irsimSUS 3014 format(1x, * 'CFS =',i5,' NIS_E =',i5,' NIS_T_N =',i5, * ' KFS =',i5,' SDS =',i5/1x, * 'GHK =',i5,' PCF =',i5,' DCS =',i5, * ' ARS_E =',i5,' ARS_R_T_N =',i5/1x, * 'GSS =',i5,' AUS =',i5,' SUS =',i5/ * 1x,'---------------------------------------------', * '---------------------------------------------------') ccc---insert (5) c*********************************************************************** c* a simple timing test to set number of simulations endogenously * c*********************************************************************** c c ENDTIM 9995 continue c c BEGINSIM c c-----asssign ID number to each of 13 methods c methCFS = 4 methNISE = 5 methNIST = 6 methKFS = 7 methSDS = 8 methGHK = 9 methPCF = 10 methDCS = 11 methARSE = 12 methARSR = 13 methGSS = 14 methAUS = 15 methSUS = 16 c c c-----initialize matrices to exact results c call matcopy(pp,16,pp0,16,1,16,1,1) call matcopy(dm,16,dm0,16,1,16,1,1) call matcopy(xlm,16,xlm0,16,1,16,1,1) call matcopy(dl,16,dl0,16,1,16,1,1) call matcopy(xll,16,xll0,16,1,16,1,1) call matcopy(rep,16,rep0,16,1,16,1,1) call matcopy(el,16,el0,16,1,16,1,1) c c-----start of simulators print '(/a/a,i5,a,f13.1/a/)', *' ********************************', *' * Start of MC Repetition ',imcrepet,' * Current seed:',curseed, *' ********************************' irec=imcrepet+1 c ccc---CFS e1 = hsec1990() ir = irsimCFS call cfs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c c-----notice that first three elements of the following vectors are c for the exact results c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','CFS: P =',pploc,' Time =',ec c rep(methCFS) = ir el(methCFS) = ec pp(methCFS) = pploc dm(methCFS) = dmloc xlm(methCFS) = xlmloc dl(methCFS) = dlloc xll(methCFS) = xllloc c c-----time per replication ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methCFS),dm(methCFS),xlm(methCFS),dl(methCFS), / xll(methCFS),rep(methCFS),el(methCFS),curseed call reswrite(cfilCFS,results,irec) c ccc---NISE e1 = hsec1990() ir = irsimNSE call nise(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','NISE: P =',pploc,' Time =',ec c rep(methNISE) = ir el(methNISE) = ec pp(methNISE) = pploc dm(methNISE) = dmloc xlm(methNISE) = xlmloc dl(methNISE) = dlloc xll(methNISE) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methNISE),dm(methNISE),xlm(methNISE),dl(methNISE), / xll(methNISE),rep(methNISE),el(methNISE),curseed call reswrite(cfilNISE,results,irec) c ccc---NIST e1 = hsec1990() ir = irsimNST call nist(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','NIST: P =',pploc,' Time =',ec c rep(methNIST) = ir el(methNIST) = ec pp(methNIST) = pploc dm(methNIST) = dmloc xlm(methNIST) = xlmloc dl(methNIST) = dlloc xll(methNIST) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methNIST),dm(methNIST),xlm(methNIST),dl(methNIST), / xll(methNIST),rep(methNIST),el(methNIST),curseed call reswrite(cfilNIST,results,irec) c ccc---KFS e1 = hsec1990() ir = irsimKFS do 250 i = 1, m 250 vec40(i) = b(i) - a(i) continue do 251 i = m+1, 2*m 251 vec40(i) = dsqrt(d(i-m)**2 + xl(i-m)**2) continue val = 0.1d0 * fmin(vec40,2*m) call fillm(parm,irmaxim+1,2,val) call pkfs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5,a,g13.5)','KFS: P =',pploc, * ' Time =',ec,' Window =',parm(methKFS,1) c rep(methKFS) = ir el(methKFS) = ec pp(methKFS) = pploc dm(methKFS) = dmloc xlm(methKFS) = xlmloc dl(methKFS) = dlloc xll(methKFS) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methKFS),dm(methKFS),xlm(methKFS),dl(methKFS), / xll(methKFS),rep(methKFS),el(methKFS),curseed call reswrite(cfilKFS,results,irec) c ccc---SDS e1 = hsec1990() ir = irsimSDS call matcopy(aw,mmaxim,w,mmaxim,1,m,1,m) call invv(aw,tempmxm,mmaxim,mmaxim,tempmx1,m, * npr,nzr,nnr,rcond,*888) call sorts(tempmx1,m) val = (0.99d0 * tempmx1(1))**2 call fillm(parm,irmaxim+1,2,val) do 260 j = 1,m do 261 i = 1,m if (i.ne.j) then tempmxm(i,j) = w(i,j) else tempmxm(i,j) = w(i,j) - val endif 261 continue 260 continue call chol(tempmxm,mmaxim,m,tempmxm,mmaxim) call sds(mmaxim,irmaxim, / m,xmu,w,wi,tempmxm,a,b,ir,unBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5,a,g13.5)','SDS: P =',pploc, * ' Time =',ec,' Lambda =',parm(methsds,1) rep(methSDS) = ir el(methSDS) = ec pp(methSDS) = pploc dm(methSDS) = dmloc xlm(methSDS) = xlmloc dl(methSDS) = dlloc xll(methSDS) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methSDS),dm(methSDS),xlm(methSDS),dl(methSDS), / xll(methSDS),rep(methSDS),el(methSDS),curseed call reswrite(cfilSDS,results,irec) c ccc---GHK e1 = hsec1990() ir = irsimGHK call ghk(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','KFS: P =',pploc,' Time =',ec c rep(methGHK) = ir el(methGHK) = ec pp(methGHK) = pploc dm(methGHK) = dmloc xlm(methGHK) = xlmloc dl(methGHK) = dlloc xll(methGHK) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methGHK),dm(methGHK),xlm(methGHK),dl(methGHK), / xll(methGHK),rep(methGHK),el(methGHK),curseed call reswrite(cfilGHK,results,irec) c ccc---PCF e1 = hsec1990() ir = irsimPCF parm(1,1) = detw call pcf(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,ussBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','GHK: P =',pploc,' Time =',ec c rep(methPCF) = ir el(methPCF) = ec pp(methPCF) = pploc dm(methPCF) = dmloc xlm(methPCF) = xlmloc dl(methPCF) = dlloc xll(methPCF) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methPCF),dm(methPCF),xlm(methPCF),dl(methPCF), / xll(methPCF),rep(methPCF),el(methPCF),curseed call reswrite(cfilPCF,results,irec) c ccc---DCS e1 = hsec1990() ir = irsimDCS parm(1,1) = detw call dcs(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,usBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5)','DCS: P =',pploc,' Time =',ec c rep(methDCS) = ir el(methDCS) = ec pp(methDCS) = pploc dm(methDCS) = dmloc xlm(methDCS) = xlmloc dl(methDCS) = dlloc xll(methDCS) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methDCS),dm(methDCS),xlm(methDCS),dl(methDCS), / xll(methDCS),rep(methDCS),el(methDCS),curseed call reswrite(cfilDCS,results,irec) c ccc---ARSE e1 = hsec1990() ir = irsimARE issd = idnint(uuBIG(1,1)*1.0d5) do 232 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(issd) ) 232 continue parm(ir+1,1) = ir call arse(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(4(a,g13.5))', * ' ARSE: P =',pploc, * ' Time =',ec,' Parm =',parm(methARSE,1), * ' % not censored:',dummy(1) c rep(methARSE) = ir el(methARSE) = ec pp(methARSE) = pploc dm(methARSE) = dmloc xlm(methARSE) = xlmloc dl(methARSE) = dlloc xll(methARSE) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methARSE),dm(methARSE),xlm(methARSE),dl(methARSE), / xll(methARSE),rep(methARSE),el(methARSE),curseed call reswrite(cfilARSE,results,irec) c ccc---ARSR e1 = hsec1990() ir = irsimARR issd = idnint(uuBIG(1,1)*1.0d5) do 233 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(issd) ) 233 continue parm(ir+1,1) = ir call arsr(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(4(a,g13.5))', * ' ARSR: P =',pploc, * ' Time =',ec,' Parm =',parm(methARSR,1), * ' % not censored:',dummy(1) c rep(methARSR) = ir el(methARSR) = ec pp(methARSR) = pploc dm(methARSR) = dmloc xlm(methARSR) = xlmloc dl(methARSR) = dlloc xll(methARSR) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methARSR),dm(methARSR),xlm(methARSR),dl(methARSR), / xll(methARSR),rep(methARSR),el(methARSR),curseed call reswrite(cfilARSR,results,irec) c ccc---GSS e1 = hsec1990() ir = irsimGSS if (gresampl.lt.0.d0) then c if < 0, this is not an endogenous timing run, so define gresampl: gresampl = idnint( 5.0d0 + dble(float(ir)/20.0d0) ) cvah added so as to make sure at least 10 gibbs resamplings are used: if(val.lt.10.d0) gresampl=10.d0 endif call fillm(parm,irmaxim+1,2,gresampl) call gss(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5,a,g13.5)','GSS: P =',pploc, * ' Time =',ec,' Parm =',parm(methGSS,1) c rep(methGSS) = ir el(methGSS) = ec pp(methGSS) = pploc dm(methGSS) = dmloc xlm(methGSS) = xlmloc dl(methGSS) = dlloc xll(methGSS) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methGSS),dm(methGSS),xlm(methGSS),dl(methGSS), / xll(methGSS),rep(methGSS),el(methGSS),curseed call reswrite(cfilGSS,results,irec) c ccc---AUS e1 = hsec1990() ir = irsimAUS do 242 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(iseed) ) parm(i,2) = 1.0d0 242 continue parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) ) parm(2,2) = idnint( 10.0d0 + dble(rsimBIG/10.0d0) ) parm(3,2) = idnint( 1.0d6*uuBIG(2,1) ) call aus(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5,a,g13.5,a,f10.1)','AUS: P =',pploc, * ' Time =',ec,' Parm =',parm(methAUS,1), * ' Max. # of CFS draws:',dummy(1) c rep(methAUS) = ir el(methAUS) = ec pp(methAUS) = pploc dm(methAUS) = dmloc xlm(methAUS) = xlmloc dl(methAUS) = dlloc xll(methAUS) = xllloc c ec = ec / dble(float(ir)) c write(results,'(8(1x,g15.9))') / pp(methAUS),dm(methAUS),xlm(methAUS),dl(methAUS), / xll(methAUS),rep(methAUS),el(methAUS),curseed call reswrite(cfilAUS,results,irec) c ccc---SUS e1 = hsec1990() ir = irsimSUS do 243 i = 1,ir parm(i,1) = idnint( 1.0d5*randu(iseed) ) parm(i,2) = 1.0d0 243 continue parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) ) parm(3,2) = idnint( 1.0d6*uuBIG(2,1) ) call sus(mmaxim,irmaxim, / m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc) call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc, / dlloc,xllloc) c if (ifullout.gt.0) *print '(1x,a,g13.5,a,g13.5,a,g13.5,a,f10.1)','SUS: P =',pploc, * ' Time =',ec,' Parm =',parm(methSUS,1), * ' Max. # of CFS draws:',dummy(1) c rep(methSUS) = ir el(methSUS) = ec pp(methSUS) = pploc dm(methSUS) = dmloc xlm(methSUS) = xlmloc dl(methSUS) = dlloc xll(methSUS) = xllloc ec = ec / float(ir) c write(results,'(8(1x,g15.9))') / pp(methSUS),dm(methSUS),xlm(methSUS),dl(methSUS), / xll(methSUS),rep(methSUS),el(methSUS),curseed call reswrite(cfilSUS,results,irec) c if (ifullout.gt.0) then c print out full results write(6,1320) 1320 format(1x/' Collected Results'/ * 1x,'-----------------------------------------------', * '------------------------------------------------'/ * 1x,' p dp/dmu(m) dln(p)/dmu(m) ', * ' dp/dl(m) dln(p)/dl(m) nrepet time(sec)'/ * 1x,'-----------------------------------------------', * '------------------------------------------------') print '(a,7g13.5)',' Quad.1:',pp(1),dm(1),xlm(1),dl(1),xll(1), * rep(1),el(1) print '(a,7g13.5)',' Quad.2:',pp(2),dm(2),xlm(2),dl(2),xll(2), * rep(2),el(2) print '(a,7g13.5)',' Exact: ',pp(3),dm(3),xlm(3),dl(3),xll(3), * rep(3),el(3) print '(a,7g13.5)',' CFS: ',pp(methcfs),dm(methcfs), * xlm(methcfs),dl(methcfs),xll(methcfs), * rep(methcfs),el(methcfs) print '(a,7g13.5)',' NISE: ',pp(methnise),dm(methnise), * xlm(methnise),dl(methnise),xll(methnise), * rep(methnise),el(methnise) print '(a,7g13.5)',' NIST: ',pp(methnist),dm(methnist), * xlm(methnist),dl(methnist),xll(methnist), * rep(methnist),el(methnist) print '(a,7g13.5)',' KFS: ',pp(methkfs),dm(methkfs), * xlm(methkfs),dl(methkfs),xll(methkfs), * rep(methkfs),el(methkfs) print '(a,7g13.5)',' SDS: ',pp(methsds),dm(methsds), * xlm(methsds),dl(methsds),xll(methsds), * rep(methsds),el(methsds) print '(a,7g13.5)',' GHK: ',pp(methghk),dm(methghk), * xlm(methghk),dl(methghk),xll(methghk), * rep(methghk),el(methghk) print '(a,7g13.5)',' PCF: ',pp(methpcf),dm(methpcf), * xlm(methpcf),dl(methpcf),xll(methpcf), * rep(methpcf),el(methpcf) print '(a,7g13.5)',' DCS: ',pp(methdcs),dm(methdcs), * xlm(methdcs),dl(methdcs),xll(methdcs), * rep(methdcs),el(methdcs) print '(a,7g13.5)',' ARSE: ',pp(metharse),dm(metharse), * xlm(metharse),dl(metharse),xll(metharse), * rep(metharse),el(metharse) print '(a,7g13.5)',' ARSR: ',pp(metharsr),dm(metharsr), * xlm(metharsr),dl(metharsr),xll(metharsr), * rep(metharsr),el(metharsr) print '(a,7g13.5)',' GSS: ',pp(methgss),dm(methgss), * xlm(methgss),dl(methgss),xll(methgss), * rep(methgss),el(methgss) print '(a,7g13.5)',' AUS: ',pp(methaus),dm(methaus), * xlm(methaus),dl(methaus),xll(methaus), * rep(methaus),el(methaus) print '(a,7g13.5)',' SUS: ',pp(methsus),dm(methsus), * xlm(methsus),dl(methsus),xll(methsus), * rep(methsus),el(methsus) write(6,3021) 3021 format(' -----------------------------------------------', * '------------------------------------------------'/) endif c c c undocumented feature: c Will stop at the end of a Monte Carlo round, if stop.//resfil c exists and has MORE than 0 bytes -- careful!. call fsize('stop.'//resfil,itemp) if (itemp.gt.0) then print *,'---------' print *,'Will stop because file stop.'//resfil//' exists.' print *,'---------' goto 999 endif if (imcrepet.lt.mcrepet) then imcrepet = imcrepet + 1 goto 1000 c DRAWRAND else goto 999 c FINISH endif 888 continue print *, 'error in inversion routine' c end of the program 999 continue c FINISH end