(* enter the following *)
nburn=100;  (* burn in iterations *)
nwrite=1;   (* 1 if a simulation file is to be written, 0 otherwise *)
hfile="sim"; (* input data file *)
outfile="mout"; (* output file name *)  
simout="mlsim"; (* name of the simulation output file if nwrite=1 *)

(* DO NOT MODIFY ANYTHING BELOW *)

Off[General::spell];
Off[General::spell1];
LD=100; LDK=100; LSET=9;
dum=0;
<<Statistics`ContinuousDistributions`; 
<<readrec.m; <<pars.m; <<dscal.m; <<dfull.m; << gau.m; <<for.m; << writerec.m;

If[nburn<0, nburn=0];
OpenRead::noopen="Cannot open input file";
Check[h=OpenRead[hfile], Quit[1]];
out=OpenWrite[outfile];                               
lwrite=True;
If[nwrite==0,lwrite=False];                               
If[lwrite,outs=OpenWrite[simout]];                               

SetStreamPosition[hfile,0];
Check[niter=Read[h, Number],
      WriteString[out, "read error in simulation input file, line 1 \n"];
      Quit[1]];
If[niter==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];
Check[npar=Read[h,Number],
      WriteString[out,"read error in simulation input file, line 1 \n"];
      Quit[1]];
If[npar==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];    

If[nburn > niter,
   WriteString[out,nburn," burn in requested but only ",
               niter, " available\n"];
   Quit[1];
  ];
nrange=niter;
repar0[];

nuse=niter-nburn;
(* Format the output header *)
WriteString [out, "Marginal Likelihood Analysis\n"];
WriteString [out, "          Total iterations: ", niter, "\n"];
WriteString [out, "        Burn-in iterations: ", nburn, "\n"];
WriteString [out, "       Analysis iterations: ", nuse, "\n\n"];

(* First pass *)
wtsum=0;
tmu=Table[0,{i,npar}];
tvar=Table[0,{i,npar},{j,npar}];
base=0;
For[iter=1, iter <= niter, iter++, 
   readrec[h,npar]; 
   If [iter<=nburn, Continue[]];
   (repar[];
    base=base+pri+pdat+nk/2;
    wt=Exp[wtlog];
    wtsum=wtsum+wt;
    For[ik=1, ik<=nk, ik++,
        aa=wt*theta[[ik]];
        tmu[[ik]]=tmu[[ik]]+aa;
        For[jk=ik, jk<=nk, jk++,
            tvar[[ik,jk]]=tvar[[ik,jk]]+aa*theta[[jk]];
           ];
       ];
   ) (* if iter> nburn *)
   ];  (* for iter=1 to niter *)      
base=base/nuse;
tmu=dscal[nk,1/wtsum,tmu];
tvar=Table[ If[ik<=jk, tvar[[ik,jk]]/wtsum-tmu[[ik]]*tmu[[jk]], tvar[[ik,jk]] ],
            {ik,nk}, {jk,nk}];
tvar=dfull[nk,tvar];
LinearSolve::nosol="Singular parameter vector variance";
tpre=Inverse[tvar];

(* prepare normal density constants *)
set1=LSET+1;
dk=nk;
adj=Table[set1/(LSET+1-iset),{iset,LSET}];
cut=Table[(0.5)*Quantile[ChiSquareDistribution[dk],1/adj[[iset]]],{iset,LSET}];
ldummy=lrange[nk,theta];
If[!lall,(
  gau[nk,tmu,tpre];
  den=Table[0,{i,LSET}]; top=Table[0,{i,LSET}];
  For [irange=1, irange <= nrange, irange++,
       gaua[nk,tmu,tpre];
       size=-gauk[nk,tmu,tpre,theta1];
       add=0;
       If [lrange[nk,theta1], add=0.5];
       If [lrange[nk,theta2], add=add+0.5];
       For [iset=1, iset <=LSET, iset++,
            If [size <= cut[[iset]], 
                  den[[iset]]=den[[iset]]+1;
                  top[[iset]]=top[[iset]]+add;
               ];
       ];  (* For iset *)
  ];  (* For irange *)
  For [iset=1, iset <=LSET, iset++,
       If [top[[iset]]<= 0, top[[iset]]=-1];
  ];
  WriteString[out,"Adjustment factors for coverage\n"];
  For[i=1, i<=LSET, i++,
             WriteString[out, e[adj[[i]],7,6]];
             WriteString[out,If [Mod[i,5]==0 || i==LSET, "\n", "     "]];
     ];   
)];

(* Construct Gelfand-Dey Posetrior moment and print *)
SetStreamPosition[hfile,0];
niter=Read[h, Number];
npar=Read[h, Number];
top=Table[0,{iset,LSET}];
den=Table[0,{iset,LSET}];
If [lwrite,WriteString[outs, nuse, " ",  LSET, "\n"]];
wtsum=0;
For[iter=1, iter <= niter, iter++, 
   readrec[h,npar]; 
   If [iter<=nburn, Continue[]];
   (repar[];
    wt=Exp[wtlog];
    wtsum=wtsum+wt;
    size=-gauk[nk,tmu,tpre,theta];
    ratio=Exp[base-size-pri-pdat];
    g=Table[0,{iset,LSET}];
    For[iset=1, iset <= LSET, iset++, 
        If [size>cut[[iset]], Continue[]];
        (g[[iset]]=adj[[iset]]*ratio;
         top[[iset]]=top[[iset]]+wt*g[[iset]];
         den[[iset]]=den[[iset]]+wt;)
       ]; 
    If [lwrite, writerec[] ];
   ) (* if iter> nburn *)
];  (* for iter=1 to niter *)      
aa=base-gaun[nk,tpre];
WriteString[out, "            Marginal Likelihood (ML) Computations\n"];
WriteString[out, "Version    Prob  
   Weighted Fractions Included    Log Ml \n"]; 
For[iset=1, iset <= LSET, iset++, 
    WriteString[out, iset, "       ", 
     PaddedForm[N[(LSET+1-iset)/set1],6], "         ", 
    PaddedForm[den[[iset]]/wtsum,6], "               ", 
    e[Log[wtsum/top[[iset]] ] + aa,8,7] , "\n"];
   ];
WriteString[out, "Inverse entries written to file must be scaled by\n"];
WriteString[out,"exp(",e[-aa,8,7],") to recover correct inverse ML\n"];

Quit[1];
Close[hfile];