c ***************************************************************** c This program assumes that there are n Y's and msamp W's, c replicated exactly twice (I keep a term numw around, but c much of the program assumes numw = 2). You have to c set n and msamp. c c The number of Gibbs samples is ngibbs c Kbigest = current number of mixtures. c Kbigsim = maximum number of mixtures, <= 5 = maxmix c c ***************************************************************** implicit real*8 (a-h,o-z) c c ***************************************************************** c These are all of size n, the sample size. The "5"=maxmix c ***************************************************************** real*8 xx(4001),yy(4001,1),pweight(4001,1) real*8 xxest(4001,1),vvnorm(4001),vvunif(4001) real*8 designx(4001,2),xxest1(4001,1) real*8 dummyx(4001) integer iiest(4001,5),ig(4001) real*8 bigsort(4001,3) integer iiselect(4001),inewsort(4001) c c ***************************************************************** c This matrix summarizes the results, and is the last thing c printed c Column #1 = n (#FFQ) c Column #2 = msamp (# in calibration) c Column #3 = irandom (1 = SRS, 0 = Stratified RS) c Column #4 = Kbigtrue (# of mixtures) c Column #5 = isnegexp (Family Type) c Column #6 = ifbeta (0 = mom, 1 = 1.0 for starters) c Column #7 = method (1 = mom, 2 = mle, 3 = Gibbs c Column #8 = ifrho (0 = beta, 1 = rho) c Column #9 = skewnorm (coefficient for skewed normal) c Column #10 = sigma_u^2 c Column #11 = mean c Column #12 = s.d. c Column #13 = median c Column #14 = mad c Column #15 = mae c Column #16 = mse c Column #17 = ase c ***************************************************************** real*8 bigsumer(6,17),smsumer(2,7) c ***************************************************************** c These matrices go with the clustering algorithm of Hartigan. c They are all of size n, the sample size c ***************************************************************** real*8 dclus(4001) integer ic2(4001) c c ***************************************************************** c These are of size msamp, the number of observations c in which replicate W's are observed c ***************************************************************** real*8 ww(397,2),yyw(397,1),dummyw(397,2) c c ***************************************************************** c These are all fixed, of size "5" = maxmix c ***************************************************************** real*8 pmixtrue(5,1),zmixtrue(5,1),smixtrue(5,1) real*8 pmixest(5,1),zmixest(5,1),smixest(5,1) real*8 trueprop(5,1) real*8 psimest(5,5),zsimest(5,5),ssimest(5,5) real*8 anslike(5,1),xbar(5),s(5) integer njmatrix(5,1) c c ***************************************************************** c These matrices are dimensioned by the number of c Gibbs samples (ngibbs). The "5" = maxmix c c The matrices ppgibbs,zmugibbs and ssgibbs are supposed to keep c the results of the Gibbs simulations, which you'll then c average to get estimates of the mixtrue probabilities c The dimension is ppgibbs(ngibbs,maxmix,maxmix). The third c aregument is the value of Kbigest. c c The matrix bigans is ngibbs x 6 x maxmix, and contains the c estimates of beta0,beta1,siguest,sigeest and rhoqt c in the current Gibbs steps c c ***************************************************************** real*8 answer(2016,6),dummyg(2016,1),bigans(2016,6,5) real*8 ppgibbs(2016,5,5),zmugibbs(2016,5,5),ssgibbs(2016,5,5) c c ***************************************************************** c These matrices are of fixed size. It is assumed that there c are to be no more than 1,000 simulated data sets! c c npopsim gives the values of the number of populations c ***************************************************************** real*8 anssim(500,6),dummysim(500),ansmom(500,5) real*8 ansmle(500,6),ansmletr(500,5) real*8 zdtzd(3,3),zdsum(3,1) integer npopsim(500,1) c c c do 53 i=1,6 do 52 j=1,17 bigsumer(i,j) = 0.0d0 52 continue 53 continue open(6,file="Butmle01.out") print *, " ############################################### " write(6,*) " Simulated Data" print *, " " rnseed = 12345 rnseed = 322398. rnseed = 54321 rnseed = 722398. rnseed = 1339003. rnseed = 7339004. rnseed = 88755439. rnseed = 7339003. do 99999 mmmqqq=1,1 print *, "mmmqqq = ",mmmqqq n = 4001 msamp = 397 numw = 2 maxmix = 5 Nbigsim = 250 c c ***************************************************************** c You allow K to be from Kestlow to Kesthigh. If you c are not estimating K, these should be equal and c Knoest should equal the common value, and iestK = 0. c If you are estimating K, set the range and iestK = 1. c ***************************************************************** iestK = 1 Kestlow = 1 Kesthigh = 3 Knoest = 1 print *, "Minimum number of Mixtures = ",Kestlow print *, "Maximum number of mixtures = ",Kesthigh c c ***************************************************************** c if irandom=1, then do simple random sampling c if irandom=0, then make the final sample have the c proportions 0.20, 0.15, 0.10, 0.15, 0.40 c ***************************************************************** irandom = 0 if (irandom.eq.0) print *, "Selection on basis of FFQ" if (irandom.eq.1) print *, "Selection at Random" c c ***************************************************************** c You can also get sampling with X not being c a mixture of normals. To do this, set c Kbigtrue = 1 and isnegexp = 1 c to get a negative exponential, and c Kbigtrue = 1, isnegexp = 2, skewnorm c to get a scaled skewed normal with parameter skewnorm and c Kbigtrue = 1 and isnegexp = 3 c to get a rescaled log(chisquared with 1 degree of freedom). c The default is mixture of normals, in which case isnegexp = 0 c ***************************************************************** skewnorm = 0.0d0 Kbigtrue = 1 isnegexp = 2 if (isnegexp.eq.0) print *, "Mixture of Normals" if (isnegexp.eq.1) print *, "X is Negative Exponential" if (isnegexp.eq.2) print *, "X is Skewed normal, " if (isnegexp.eq.2) print *, " with parameter",skewnorm if (isnegexp.eq.3) print *, "X is log(Chisquared(1))" c ***************************************************************** c You can decide to turn off the Gibbs sampler. c It runs if ifgibbs=1, and doesn't run if ifgibbs=0 c ***************************************************************** ifgibbs=1 c ***************************************************************** c You can decide to turn start the sampler for beta1 either c at the method of moments (ifbeta = 0) or at 1.0 (ifbeta = 1) c ***************************************************************** ifbeta = 1 if (ifbeta.eq.0) print *, "Start beta_1 at MOM" if (ifbeta.eq.1) print *, "Start beta_1 at 1.0" c ***************************************************************** c You have to set the prior variance for beta_1, c which is ssbprior c ***************************************************************** ssbprior = 0.25 ** 2 c c c c do 56 i=1,6 bigsumer(i,1) = dfloat(n) bigsumer(i,2) = dfloat(msamp) bigsumer(i,3) = dfloat(irandom) bigsumer(i,4) = dfloat(Kbigtrue) bigsumer(i,5) = dfloat(isnegexp) bigsumer(i,6) = dfloat(ifbeta) bigsumer(i,8) = 0.0d0 bigsumer(i,9) = skewnorm 56 continue bigsumer(1,7) = 1.0d0 bigsumer(2,7) = 1.0d0 bigsumer(3,7) = 2.0d0 bigsumer(4,7) = 2.0d0 bigsumer(5,7) = 3.0d0 bigsumer(6,7) = 3.0d0 print *, "Initial seed = ",rnseed print *, "Sample size = ",n print *, "# Replicated W's = ",msamp print *, "# of simulations = ",Nbigsim ngibbs = 2016 nuprior = 1 Mprior = 5 cycper = 0.05 aeps = 5. auu = 5. nueps = 5. nuuu = 5. c ************************************************************* c ************************************************************* c ************************************************************* c **** The program runs automatically from here ***** c ************************************************************* c ************************************************************* c ************************************************************* print *, "Number of Gibbs samples = ",ngibbs print *, "Prior value of nu = ",nuprior print *, "Prior value of M = ",Mprior print *, "Cycling percentage = ",cycper print *, "Value of Aeps = ",aeps print *, "Value of Auu = ",auu print *, "Value of nu_eps = ",nueps print *, "Value of nu_uu = ",nuuu print *, "Value of sigma^2_{beta} = ",ssbprior if (iestK.eq.0) then print *, "K not estimated, fixed = ",Knoest endif if (iestK.eq.1) then print *, "K is estimated *********" endif print *, " #################################################### " print *, " " c c c Set up the Anderson Darling Tests anddar15 = 0. anddar10 = 0. anddar05 = 0. call zeros(psimest,maxmix,maxmix) call zeros(zsimest,maxmix,maxmix) call zeros(ssimest,maxmix,maxmix) call zeros(xx,n,1) call zeros(xxest,n,1) call zeros(yy,n,1) call zeros(ww,msamp,numw) call zeros(pmixest,maxmix,1) call zeros(smixest,maxmix,1) call zeros(zmixest,maxmix,1) call zeros(vvunif,n,1) call zeros(dummyx,n,1) call zeros(dummyw,msamp,numw) call zeros(anssim,500,6) call zeros(ansmom,500,5) call zeros(ansmle,500,6) call zeros(ansmletr,500,5) call izeros(npopsim,500,1) c ***************************************************************** c Set the true valuies for the simulation c ***************************************************************** call settrue(smixtrue,zmixtrue,pmixtrue,maxmix, 6 b0true,b1true,sigetrue,sigutrue,njmatrix, 6 Kbigtrue,rhoqtrue) c c ******************************************************************** c Start the simulation Here c ******************************************************************** c print *, " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" do 9800 ksim=1,Nbigsim print *, " Simulation Number = ",ksim call izeros(iiest,n,maxmix) c ******************************************************************** c Generate the data and the sampling weights c ******************************************************************** call gendata(iiest,xxest,yy,ww,n,msamp,numw,rnseed, 6 vvunif,maxmix,pmixtrue,zmixtrue,smixtrue, 6 vvnorm,b0true,b1true,sigutrue,sigetrue, 6 yyw,irandom,trueprop,Kbigtrue,designx,bigsort, 6 iiselect,inewsort,dummyx,pweight,isnegexp, 6 skewnorm) do 57 i=1,n 57 xx(i) = xxest(i,1) c c ******************************************************************** c Get moments estimate of sigma_u c ******************************************************************** siguest = findsigu(ww,msamp,numw) c ***************************************************************** c Get method of moments estimates for beta0, beta1, c and sigma_e while accounting for sampling weights c ***************************************************************** call mom2(yy,ww,siguest,pweight,n,msamp,numw,dummyw, 6 sigeest,beta0,beta1,sigxmom,zmuxmom,ierrcde) beta0mom = beta0 beta1mom = beta1 sigumom = siguest sigemom = sigeest zmuxmle = zmuxmom sigxmle = sigxmom c c ******************************************************************** c Compute Moments estimate of rhoqt c ******************************************************************** varx = sigxmom ** 2 rhoqtmom = (beta1 * varx) 6 / (dsqrt(varx) 6 * dsqrt( (beta1 * beta1 * varx) 6 + (sigeest ** 2))) c c ******************************************************************** c Put the answers for moments into ansmom c ******************************************************************** ansmom(ksim,1) = sigeest ** 2 ansmom(ksim,2) = siguest ** 2 ansmom(ksim,3) = beta0 ansmom(ksim,4) = beta1 ansmom(ksim,5) = rhoqtmom ansmletr(ksim,1) = sigeest ** 2 ansmletr(ksim,2) = siguest ** 2 ansmletr(ksim,3) = beta0 ansmletr(ksim,4) = beta1 ansmletr(ksim,5) = rhoqtmom c print *, "MOM for beta0 = ",beta0 c print *, "MOM for beta1 = ",beta1 c print *, "MOM fo rhoqt = ",rhoqtmom c print *, "MOM for sigma_u = ",siguest c print *, "MOM for sigma_e = ",sigeest c print *, "MOM for sigma_x = ",sigxmom c ********************************************************************* c Compute the actual mle c ********************************************************************* call zeros(zdtzd,3,3) call zeros(zdsum,3,1) do 100 i=1,msamp dummyw(i,1) = ( (ww(i,1) + ww(i,2)) / 2.) zdsum(1,1) = zdsum(1,1) + yy(i,1) zdsum(2,1) = zdsum(2,1) + ww(i,1) zdsum(3,1) = zdsum(3,1) + ww(i,2) zdtzd(1,1) = zdtzd(1,1) + (yy(i,1) * yy(i,1)) zdtzd(1,2) = zdtzd(1,2) + (yy(i,1) * ww(i,1)) zdtzd(1,3) = zdtzd(1,3) + (yy(i,1) * ww(i,2)) zdtzd(2,2) = zdtzd(2,2) + (ww(i,1) * ww(i,1)) zdtzd(2,3) = zdtzd(2,3) + (ww(i,1) * ww(i,2)) zdtzd(3,3) = zdtzd(3,3) + (ww(i,2) * ww(i,2)) 100 continue zdtzd(2,1) = zdtzd(1,2) zdtzd(3,1) = zdtzd(1,3) zdtzd(3,2) = zdtzd(2,3) yotsum = 0. yotysum = 0. if (msamp.lt.n) then do 200 i=(msamp+1),n yotsum = yotsum + (yy(i,1)) yotysum = yotysum + (yy(i,1) ** 2) 200 continue endif th1 = beta0mom + (beta1mom * zmuxmom) th2 = zmuxmom th3 = ((beta1mom * sigxmom) ** 2) + (sigemom ** 2) th4 = beta1mom * (sigxmom ** 2) th6 = (sigxmom ** 2) + (sigumom ** 2) th7 = sigxmom ** 2 ierrcde2 = 0 call fullmlso(th1,th2,th3,th4,th6,th7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,n,msamp,ierrcde2,zmuxmle, 6 sigxmle,sigumle,sigemle,beta0mle,beta1mle,rhoqtmle) c print *, "MLE for mu_x = ",zmuxmle c print *, "MLE for beta0 = ",beta0mle c print *, "MLE for beta1 = ",beta1mle c print *, "MLE fo rhoqt = ",rhoqtmle c print *, "MLE for sigma_u = ",sigumle c print *, "MLE for sigma_e = ",sigemle c print *, "MLE for sigma_x = ",sigxmle c c ******************************************************************** c Put the answers for actual mle into ansmlrtr c ******************************************************************** if (ierrcde2.eq.0) then ansmletr(ksim,1) = (sigemle ** 2) ansmletr(ksim,2) = (sigumle ** 2) ansmletr(ksim,3) = beta0mle ansmletr(ksim,4) = beta1mle ansmletr(ksim,5) = rhoqtmle endif c ***************************************************************** c For Kbigest=1, use the mle. Saves on computing time c ***************************************************************** do 43300 i=1,ngibbs bigans(i,1,1) = ansmletr(ksim,1) bigans(i,2,1) = ansmletr(ksim,2) bigans(i,3,1) = ansmletr(ksim,3) bigans(i,4,1) = ansmletr(ksim,4) bigans(i,5,1) = ansmletr(ksim,5) bigans(i,6,1) = ansmletr(ksim,5) do 43200 j=1,maxmix ppgibbs(i,j,1) = 1. zmugibbs(i,j,1) = zmuxmle ssgibbs(i,j,1) = sigxmle 43200 continue 43300 continue c c ******************************************************************** c Now you start the Bayesian analysis, selecting the c current number of populations Kbigest c ******************************************************************** if (ifgibbs.eq.1) then do 8950 Kbigest = Kestlow,Kesthigh beta0 = beta0mom beta1 = beta1mom if (ifbeta.eq.1) beta1 = 1.0d0 siguest = sigumom sigeest = sigemom c c c ******************************************************************** c Get starting values or iiest, zmixest, smixest, c pmixest, iiest, and xxest. This uses a cluster algorithm c on Y to get the first four, and then the usual c linear regression formulae to get xxest c c ******************************************************************** if (kbigest.eq.1) go to 8949 call clustyy(yy,dummyx,n,Kbigest,pmixest,smixest, 6 zmixest,njmatrix,maxmix,dclus,ig,ic2, 6 iiest,beta0,beta1,sigeest,siguest,ksim, 6 zmuxmom) kklow = 1 kkhigh = n call rnorm(rnseed,n,vvnorm) call findxxxx(yy,ww,sigeest,siguest,n,msamp,numw,beta0,beta1, 6 pmixest,zmixest,smixest,iiest,vvnorm,kbigest, 6 xxest,kklow,kkhigh,maxmix,xxest1) do 4539 kk=1,Kbigest xbar(kk) = zmixest(kk,1) 4539 s(kk) = smixest(kk,1) do 4541 i=1,n 4541 xx(i) = xxest(i,1) kkgibbs = 1 c c ******************************************************************** c Cycle through the Gibbs step once in order to get c the population indicators in line c ******************************************************************** call gibbloop(xx,xxest,ig,smixest,pmixest, 6 zmixest,xbar,s,n,Mprior,nuprior,rnseed, 6 maxmix,Kbigest,njmatrix,iiest, 6 kklow,kkhigh,kkdelta,kkgibbs) c c ******************************************************************** c Start the Gibbs Estimation c ******************************************************************** kklow = 1 kkhigh = n * cycper kkdelta = kkhigh - kklow + 1 do 8900 kkgibbs=1,ngibbs c n1 = kkgibbs / 200 c n2 = (kkgibbs+1) / 200 c if(n2.gt.n1) then c print *, " " c print *, " " c print *, " " c print *, " &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" c print *, "Gibbs step = ",kkgibbs,beta0,beta1,siguest,sigeest c print *, "Kbigest = ",Kbigest c print *, "pmixest = ",(pmixest(j,1),j=1,Kbigest) c print *, "njmatrix = ",(njmatrix(j,1),j=1,Kbigest) c print *, "zmixest = ",(zmixest(j,1),j=1,Kbigest) c print *, "smixest = ",(smixest(j,1),j=1,Kbigest) c print *, "b0,b1,su,se,rho = ",beta0,beta1,siguest, c 6 sigeest,rhoqt c print *, " &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" c endif bigans(kkgibbs,1,Kbigest) = sigeest ** 2 bigans(kkgibbs,2,Kbigest) = siguest ** 2 bigans(kkgibbs,3,Kbigest) = beta0 bigans(kkgibbs,4,Kbigest) = beta1 c if (ksim.ge.25) print *,"Calling findvxex" call findvxex(maxmix,pmixest,zmixest,smixest, 6 Kbigest,expx,varx) rhoqt = (beta1 * varx) 6 / (dsqrt(varx) 6 * dsqrt( (beta1 * beta1 * varx) 6 + (sigeest ** 2))) bigans(kkgibbs,5,Kbigest) = rhoqt sigy = stdc(yy,n,1,1) rhoqtnew = beta1 * varx / (sigy * dsqrt(varx)) bigans(kkgibbs,6,Kbigest) = rhoqtnew c ***************************************************************** c Update the parameters c ***************************************************************** c c ****************************************************************** c Update sigma_u and sigma_e c ****************************************************************** c if (ksim.ge.25) print *,"Calling findsesu" call findsesu(sigeest,siguest,n,nueps,msamp,nuuu,aeps, 6 yy,beta0,beta1,xxest,ww,auu,dummyx,dummyw,numw, 6 rnseed) c c ****************************************************************** c Update beta_0 and beta_1 c ****************************************************************** c if (ksim.ge.25) print *,"Calling findbeta" call findbeta(rnseed,n,xxest,yy,sigeest,designx, 6 beta0,beta1,ssbprior) c ***************************************************************** c Update the mixture estimates c ***************************************************************** ccc if (Kbigest.eq.1) then ccc zmixest(1,1) = zmeanc(xxest,n,1,1) ccc smixest(1,1) = stdc(xxest,n,1,1) ccc endif ccc if (kbigest.ge.2) then c if (ksim.ge.25) print *,"Calling gibbloop" call gibbloop(xx,xxest,ig,smixest,pmixest, 6 zmixest,xbar,s,n,Mprior,nuprior,rnseed, 6 maxmix,Kbigest,njmatrix,iiest, 6 kklow,kkhigh,kkdelta,kkgibbs) ccc endif c c ****************************************************************** c Update the estimated value of x c ****************************************************************** call rnorm(rnseed,n,vvnorm) c if (ksim.ge.25) print *,"Calling findxxxx" call findxxxx(yy,ww,sigeest,siguest,n,msamp,numw,beta0,beta1, 6 pmixest,zmixest,smixest,iiest,vvnorm,kbigest, 6 xxest,kklow,kkhigh,maxmix,xxest1) c if (ksim.ge.25) print *,"Finished findxxxx" do 4549 i=1,n 4549 xx(i) = xxest(i,1) c c ************************************************************ c Put the results of the Gibbs step for the mixture parameters c into the appropriate matrix c ************************************************************ do 3200 j=1,maxmix ppgibbs(kkgibbs,j,Kbigest) = pmixest(j,1) zmugibbs(kkgibbs,j,Kbigest) = zmixest(j,1) ssgibbs(kkgibbs,j,Kbigest) = smixest(j,1) 3200 continue 8900 continue c c ************************************************************ c Finished the Gibbs sampling for this value of Kbigest c ************************************************************ c 8949 continue 8950 continue c print *, "******************************* " c print *, "rnseed = ",rnseed c print *, "******************************* " c c ************************************************************ c You have now finished the Gibbs sampling for c all values of Kbigest, and it is now time to try to c estimate the optimal value c ************************************************************ c c c ************************************************************ c The first step is to get the estimates of the mixture c parameters, and put them into matrices psimest,zsimest c and ssimest. Also, you want to get the estimates c of beta0, beta1, sigma_e^2 and sigma_u^2. c ************************************************************ if (iestK.eq.1) then do 9203 j=1,maxmix do 9201 k=1,maxmix psimest(j,k) = st3mean(ppgibbs,ngibbs,maxmix,maxmix,j,k) zsimest(j,k) = st3mean(zmugibbs,ngibbs,maxmix,maxmix,j,k) ssimest(j,k) = st3mean(ssgibbs,ngibbs,maxmix,maxmix,j,k) 9201 continue 9203 continue do 6933 Kbigest=Kestlow,Kesthigh do 6925 j=1,maxmix pmixest(j,1) = psimest(j,Kbigest) zmixest(j,1) = zsimest(j,Kbigest) smixest(j,1) = ssimest(j,Kbigest) 6925 continue do 6931 i=1,ngibbs do 6929 j=1,4 6929 answer(i,j) = bigans(i,j,Kbigest) 6931 continue sigseest = zmedian(answer,dummyg,ngibbs,6,1) sigsuest = zmedian(answer,dummyg,ngibbs,6,2) beta0 = zmedian(answer,dummyg,ngibbs,6,3) beta1 = zmedian(answer,dummyg,ngibbs,6,4) c ************************************************************ c Compute the likelihood function for this number of c mixtures c ************************************************************ sigu = dsqrt(sigsuest) sige = dsqrt(sigseest) c ********************************************************************* c Compute the loglikelihood of the data and penalize c it for the number of paramaters, which is c dd = 4 + (2 * Kbigest) + (Kbigest - 1) c c The penalty is to subtract dd * dlog(dfloat(n)) c from the loglikelihood. c ********************************************************************* call cloglike(anslike,n,msamp,numw,maxmix, 6 pmixest,smixest,zmixest, 6 Kbigest,yy,ww,beta0,beta1,sigu,sige, 6 dummyx) 6933 continue c c ************************************************************ c You now have the likelihood contributions for all c possible combinations, and it is time to estimate c Kest c c ************************************************************ print *, " " andtest = andar(yy,dummyx,vvunif,vvnorm,n) if (andtest.ge.0.57d0) anddar15=anddar15+(1.0d0/Nbigsim) if (andtest.ge.0.64d0) anddar10=anddar10+(1.0d0/Nbigsim) if (andtest.ge.0.77d0) anddar05=anddar05+(1.0d0/Nbigsim) print *, "A-D test statistic = ",andtest print *, "The # Pops. and loglikelihood = " do 2241 j=1,maxmix 2241 print *, " ",j,anslike(j,1) call rightmix(anslike,Kestlow,Kesthigh,maxmix,Kest) endif if (iestK.eq.0) Kest = Knoest c c ************************************************************ c Record in the simulation the number of subpopulations c you have estimated are correct c ************************************************************ npopsim(ksim,1) = Kest print *, "For ksim = ",ksim," # populations = ",Kest c c ************************************************************ c Compute the Gibbs estimates of the parameters c c The first step is to put into "answer" the Gibbs c steps for the selected number of mixtures, Kest c ************************************************************ do 7831 i=1,ngibbs do 7829 j=1,6 7829 answer(i,j) = bigans(i,j,Kest) 7831 continue sigseest = zmedian(answer,dummyg,ngibbs,6,1) sigsuest = zmedian(answer,dummyg,ngibbs,6,2) beta0 = zmedian(answer,dummyg,ngibbs,6,3) beta1 = zmedian(answer,dummyg,ngibbs,6,4) rhoqt = zmedian(answer,dummyg,ngibbs,6,5) rhoqtnew = zmedian(answer,dummyg,ngibbs,6,6) anssim(ksim,1) = sigseest anssim(ksim,2) = sigsuest anssim(ksim,3) = beta0 anssim(ksim,4) = beta1 anssim(ksim,5) = rhoqt anssim(ksim,6) = rhoqtnew c c ***************************************************************** c If Kestlow=1, then you've computed the ordinary Bayes c estimate for 1 population, so store those results c ***************************************************************** c if (Kestlow.eq.1) then do 8331 i=1,ngibbs do 8329 j=1,6 8329 answer(i,j) = bigans(i,j,1) 8331 continue sigseest = zmedian(answer,dummyg,ngibbs,6,1) sigsuest = zmedian(answer,dummyg,ngibbs,6,2) beta0 = zmedian(answer,dummyg,ngibbs,6,3) beta1 = zmedian(answer,dummyg,ngibbs,6,4) rhoqt = zmedian(answer,dummyg,ngibbs,6,5) rhoqtnew = zmedian(answer,dummyg,ngibbs,6,6) ansmle(ksim,1) = sigseest ansmle(ksim,2) = sigsuest ansmle(ksim,3) = beta0 ansmle(ksim,4) = beta1 ansmle(ksim,5) = rhoqt ansmle(ksim,6) = rhoqtnew endif endif 9800 continue bigsumer(1,10) = sigutrue ** 2 bigsumer(2,10) = sigutrue ** 2 bigsumer(3,10) = sigutrue ** 2 bigsumer(4,10) = sigutrue ** 2 bigsumer(5,10) = sigutrue ** 2 bigsumer(6,10) = sigutrue ** 2 c c ***************************************************************** c Summarize the simulation for method of moments c ***************************************************************** call momsum(ansmom,dummysim,Nbigsim,smsumer) print *, "True value of rho_{QT} = ",rhoqtrue bigsumer(1,8) = 0.0d0 bigsumer(2,8) = 1.0d0 do 9801 i=1,7 bigsumer(1,i+10) = smsumer(1,i) bigsumer(2,i+10) = smsumer(2,i) 9801 continue c c ***************************************************************** c Summarize the simulation for True MLE c ***************************************************************** call mletrsum(ansmletr,dummysim,Nbigsim,smsumer) print *, "True value of rho_{QT} = ",rhoqtrue bigsumer(3,8) = 0.0d0 bigsumer(4,8) = 1.0d0 do 9802 i=1,7 bigsumer(3,i+10) = smsumer(1,i) bigsumer(4,i+10) = smsumer(2,i) 9802 continue if (ifgibbs.eq.1) then c c ***************************************************************** c Summarize the simulation for Bayes c ***************************************************************** call simsum(anssim,dummysim,Nbigsim,smsumer) print *, "True value of rho_{QT} = ",rhoqtrue bigsumer(5,8) = 0.0d0 bigsumer(6,8) = 1.0d0 do 9803 i=1,7 bigsumer(5,i+10) = smsumer(1,i) bigsumer(6,i+10) = smsumer(2,i) 9803 continue c c ***************************************************************** c Summarize the number of populations chosen c ***************************************************************** call simpop(npopsim,Nbigsim) endif print *, " " print *, "Level of A-D nominal .15 test = ",anddar15 print *, "Level of A-D nominal .10 test = ",anddar10 print *, "Level of A-D nominal .05 test = ",anddar05 print *, " " print *, "Final seed = ",rnseed print *, " " print *, " " print *, "Here are the final answers, in terms of bigsumer" print *, "Column #1 = n " print *, "Column #2 = msamp " print *, "Column #3 = irandom " print *, "Column #4 = Kbigtrue " print *, "Column #5 = isnegexp " print *, "Column #6 = ifbeta " print *, "Column #7 = method " print *, "Column #8 = ifrho " print *, "Column #9 = skewnorm " print *, "Column #10 = sigma_u^2 " print *, "Column #11 = mean " print *, "Column #12 = s.d. " print *, "Column #13 = median " print *, "Column #14 = mad " print *, "Column #15 = mae " print *, "Column #16 = mse " print *, "Column #17 = ase " do 99993 i=1,6 write(6,99991) (bigsumer(i,j),j=1,17) 99991 format(2F6.0,6F3.0,9F10.6) 99993 continue 99999 continue stop end c ***************************************************************** c Set the true valuies for the simulation c ***************************************************************** subroutine settrue(smixtrue,zmixtrue,pmixtrue,maxmix, 6 b0true,b1true,sigetrue,sigutrue,njmatrix, 6 Kbigtrue,rhoqtrue) implicit real*8(a-h,o-z) real*8 smixtrue(maxmix,1),zmixtrue(maxmix,1) real*8 pmixtrue(maxmix,1) integer njmatrix(maxmix,1) b0true = 5.95 b1true = 0.83 sigutrue = dsqrt(30.36d0) sigutrue = dsqrt(83.35d0) sigetrue = dsqrt(40.92d0) do 100 j=1,maxmix njmatrix(j,1) = 0 smixtrue(j,1) = 0. zmixtrue(j,1) = 0. 100 pmixtrue(j,1) = 0. if (kbigtrue.eq.1) then pmixtrue(1,1) = 1. smixtrue(1,1) = dsqrt(24.45d0) zmixtrue(1,1) = 38.25 endif if (kbigtrue.eq.2) then pmixtrue(1,1) = .5 pmixtrue(2,1) = .5 smixtrue(1,1) = 5. smixtrue(2,1) = 5. zmixtrue(1,1) = 33.25 zmixtrue(2,1) = 43.25 endif if (kbigtrue.eq.3) then pmixtrue(1,1) = .45 pmixtrue(2,1) = .30 pmixtrue(3,1) = .25 smixtrue(1,1) = dsqrt(16.04d0) smixtrue(2,1) = dsqrt(16.04d0) smixtrue(3,1) = dsqrt(16.04d0) zmixtrue(1,1) = 35.777778 zmixtrue(2,1) = 38. zmixtrue(3,1) = 43. endif do 200 j=1,Kbigest 200 njmatrix(j,1) = n * pmixtrue(j,1) print *, "True value of sigma_e = ",sigetrue print *, "True value of sigma_u = ",sigutrue print *, "True value of beta_0 = ",b0true print *, "True value of beta_1 = ",b1true print *, "Kbigtrue = ",Kbigtrue print *, "zmixtrue = ", 6 (zmixtrue(j,1),j=1,Kbigtrue) print *, "pmixtrue = ", 6 (pmixtrue(j,1),j=1,Kbigtrue) print *, "smixtrue = ", 6 (smixtrue(j,1),j=1,Kbigtrue) call findvxex(maxmix,pmixtrue,zmixtrue,smixtrue, 6 Kbigtrue,expx,varx) sstdx = dsqrt(varx) print *, "Here is varx = ",varx print *, "Here is stdx = ",sstdx print *, "Here is expx = ",expx rhoqtrue = (b1true * varx) 6 / (dsqrt(varx) * 6 dsqrt( (b1true * b1true * varx) 6 + (sigetrue ** 2))) print *, "True value of rho_{QT} = ",rhoqtrue return end