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(271),yy(271,1),pweight(271,1)
        real*8  xxest(271,1),vvnorm(271),vvunif(271)
        real*8  designx(271,2),xxest1(271,1)
        real*8  dummyx(271)
        real*8  bigdata(271,9)
        integer iiest(271,5),ig(271)

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(271)
        integer ic2(271)
c
c       *****************************************************************
c       These are of size msamp, the number of observations
c         in which replicate W's are observed
c       *****************************************************************
        real*8  ww(271,2),yyw(271,1),dummyw(271,2)
c
c       *****************************************************************
c       These are all fixed, of size "5" = maxmix
c       *****************************************************************
        real*8  pmixest(5,1),zmixest(5,1),smixest(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(1016,6),dummyg(1016,1),bigans(1016,6,5)
        real*8  ppgibbs(1016,5,5),zmugibbs(1016,5,5),ssgibbs(1016,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
        open(6,file="wish.pccal01.out")
        print *, " ############################################### "
        write(6,*) "   Wish % Calories Data"
        print *, " "
c       ********************************************************************
c       Set n and msamp. Generate the data and the sampling weights
c       ********************************************************************
        open(3,file="wish.mcpccal.dtt")
        n         = 271
        msamp     = 271
        do 7 i=1,n
        read(3,5) (bigdata(i,j),j=1,9)
5       format(F8.2,8F9.2)
7       continue
        do 23 i=1,n
        yy(i,1) = bigdata(i,1)
        ww(i,1) = (bigdata(i,2) + bigdata(i,3) + bigdata(i,4))/3
        ww(i,2) = (bigdata(i,5) + bigdata(i,6) + bigdata(i,7))/3
        pweight(i,1) = 1.
 23     xxest(i,1) = (ww(i,1) + ww(i,2)) / 2
           do 57 i=1,n
 57        xx(i) = xxest(i,1)
        aa = zmeanc(yy,n,1,1)
        bb = zmeanc(xxest,n,1,1)
        print *, " "
        print *, "*************************************** "
        print *, " Mean of FFQ = ",aa
        print *, " Mean of  FR = ",bb
        print *, "*************************************** "
        print *, " "
        do 99905 kkdat=1,2
        do 99903 kkprior=1,2
        do 99901 kkbeta=1,2
c       ********************************************************************
c       set ifdatyy=1 if you cluster on Y
c       set ifdatyy=0 if you cluster on Wbar
c       ********************************************************************
        ifdatyy = kkdat-1
        if (ifdatyy.eq.0) print *, "Cluster on basis of Wbar"
        if (ifdatyy.eq.1) print *, "Cluster on basis of Y"
c       ********************************************************************
c       Set the prior variance for beta_1
c       ********************************************************************
        if (kkprior.eq.1)        ssbprior   = .25 ** 2
        if (kkprior.eq.2)        ssbprior   = 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 = kkbeta-1
        if (ifbeta.eq.0) print *, "Start beta_1 at MOM"
        if (ifbeta.eq.1) print *, "Start beta_1 at 1.0"
        rnseed = 12345
        rnseed = 322398.
        rnseed = 54321
        rnseed = 722398.
        rnseed = 88755439.
        rnseed = 7339003.
        rnseed = 1339003.
        numw      = 2
        maxmix    = 5
        Nbigsim   = 1
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  = 2
        Knoest    = 1
        print *, "Minimum number of Mixtures = ",Kestlow
        print *, "Maximum number of mixtures = ",Kesthigh
        irandom = 0
c       *****************************************************************
c       You can decide to turn off the Gibbs sampler.
c       if runs if ifgibbs=1, and doesn't run if ifgibbs=0
c       *****************************************************************
        ifgibbs=1
c
c
c
c
        print *, "Initial seed            = ",rnseed
        print *, "Sample size             = ",n
        print *, "# Replicated W's        = ",msamp
        print *, "# of simulations        = ",Nbigsim
        ngibbs    = 1016
        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(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       ********************************************************************
c       Start the simulation Here
c       ********************************************************************
c
        print *, " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
        do 9800 ksim=1,Nbigsim
        print *, "      Simulation Number = ",ksim
c        call izeros(iitrue,n,maxmix)
        call izeros(iiest,n,maxmix)
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
        print *, "MOM for beta0      = ",beta0
        print *, "MOM for beta1      = ",beta1
        print *, "MOM fo  rhoqt      = ",rhoqtmom
        print *, "MOM for sigma_u    = ",siguest
        print *, "MOM for sigma_e    = ",sigeest
        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
         if (ifdatyy.eq.0) then
         beta0 = 0.
         beta1 = 1.
         do 3031 i=1,n
 3031       yyw(i,1) = (ww(i,1) + ww(i,2)) / 2.
            endif
         if (ifdatyy.eq.1) then
         beta0 = beta0mom
         beta1 = beta1mom
         do 3033 i=1,n
 3033       yyw(i,1) = yy(i,1)
            endif
         call clustyy(yyw,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)
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,
c     6             sigeest,zmuxmom,sigxmom
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
         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      ******************************************************************
       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      ******************************************************************
        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
           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)
        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 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
        print *, "******************************* "
        print *, "rnseed = ",rnseed
        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
c
c      *****************************************************************
c      Summarize the simulation for method of moments
c      *****************************************************************
        call momsum(ansmom,dummysim,Nbigsim)
        print *, "True value of rho_{QT}  = ",rhoqtrue
c
c      *****************************************************************
c      Summarize the simulation for True MLE
c      *****************************************************************
        call mletrsum(ansmletr,dummysim,Nbigsim)
        print *, "True value of rho_{QT}  = ",rhoqtrue
       if (ifgibbs.eq.1) then
c
c      *****************************************************************
c      Summarize the simulation for Bayes
c      *****************************************************************
        call simsum(anssim,dummysim,Nbigsim)
        print *, "True value of rho_{QT}  = ",rhoqtrue
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 *, " "
99901  continue
99903  continue
99905  continue
       stop
       end