/* 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="mlout"; /* output file name */
simout="mlsim";  /* name of the simulation output file if nwrite=1 */

      /* do not write anything below */
let LD=100; let LDK=100; let LSET=9;
dum=0;
if nburn < 0;
   nburn=0;
endif;

if nwrite;
   output file=^simout reset; 
endif;

load in[]=^hfile;
niter=in[1];
npar=in[2];
if not (rows(in) == (niter*(4+npar) + 2));
   print "
          Error in simulation input file line 1 
          or unexpected end of file reached at line " (niter*npar +2) "
          or some of the data is not numeric \n";
   end;
endif;
par=zeros(1,npar);    
if nburn > niter;
   print nburn "burn in requested but only " niter "available \n";
endif;
nrange=niter;
call repar0;
nuse=niter-nburn;

/* First pass */
wtsum=0;
tmu=zeros(1,npar);
tvar=zeros(npar,npar);
base=0;
jter=0; wtlog=0; pri=0; pdat=0; theta=0; nk=0;
k=(4+npar)*nburn + 3;
for iter (nburn+1,niter,1);
   call readrec;
   call repar;
   base=base+pri+pdat+nk/2;
   wt=exp(wtlog);
   wtsum=wtsum+wt;
   for ik (1, nk, 1);
      aa=wt*theta[ik];
      tmu[ik]=tmu[ik]+aa;
      for jk (1, nk, 1); 
            tvar[ik,jk]=tvar[ik,jk]+aa*theta[jk];
      endfor;
   endfor;
endfor; /* iter */
base=base/nuse;
tmu=tmu[1:nk];
tvar=tvar[1:nk,1:nk];
tmu=tmu/wtsum;
tvar=tvar/wtsum-tmu'*tmu;
tpre=invpd(tvar);

/* prepare normal density constants */
set1=LSET+1;
dk=nk;
adj=((seqa(LSET,-1,LSET))^(-1)*set1)';
cut=((0.5)*cdfchii(adj^(-1),dk+zeros(1,LSET)));
lall=0;
ldummy=lrange(nk,theta);  
theta1=zeros(1,nk); theta2=zeros(1,nk);
if not lall;
  a1=zeros(nk,1); v1=zeros(nk,1);
  call gau(nk,tmu,tpre);
  top=zeros(1,LSET);
  den=zeros(1,LSET);
  for irange (1,nrange,1);
     call gaua(nk,tmu,tpre);
     size=(0.5)*(theta1-tmu)*tpre*(theta1-tmu)';
     add=0;
     if lrange(nk,theta1);
        add=.5;
     endif;
     if lrange(nk,theta2);
        add=add+.5;
     endif;
/* can next 4 lines be moved outside the loop */
     den=den+(cut .>= size);
     top=top+(cut .>= size)*add;
  endfor; /* irange */
  top=-(top .<= 0) +top.*(top.>0);
  adj=adj.*den./top;
   
endif;
/* section 2b to be inserted here */

/*construct Gelfand-Dey posterior moment and print */
top=zeros(1,LSET);
den=zeros(1,LSET);
if nwrite; 
   format /rdn 16,0;
   print nuse;; print LSET;
endif; 
wtsum=0;
k=(4+npar)*nburn + 3;
for iter (nburn+1,niter,1);
   call readrec;
   call repar;
   wt=exp(wtlog);
   wtsum=wtsum+wt;
   size=(0.5)*(theta-tmu)*tpre*(theta-tmu)';
   ratio=exp(base-size-pri-pdat);
   g=(cut .>= size)*ratio .* adj;
   top=top+(cut .>= size)*wt .* g;
   den=den+(cut .>= size)*wt;
   if nwrite; 
      call writerec(iter); 
   endif;
endfor; /* iter */
aa=base+(0.5)*nk*ln(2*pi)-(0.5)*ln(abs(det(tpre)));

output file=^outfile reset;
/* Format the Output header */
print "Marginal Likelihood Analysis";
format /rdn 10, 0;
print "Total iterations available:" niter;
print "        Burn-in iterations:" nburn;
print "       Analysis iterations:" nuse;
print;
if not lall;
   print "Adjustment factors for coverage ";
   format /re 16,7; print adj;
   print;
endif;
print "      Marginal Likelihood (ML) Computations ";
print "Version    Prob     Weighted Fractions Included    Log Ml";
for iset (1,LSET,1);
   format /rdn 3,0; print iset;;
   format /rd 9,6; print ((LSET+1-iset)/set1);;
   format /rd 25,6; print den[iset]/wtsum;;
   format /re 25,7; print (ln(wtsum/top[iset])+aa);
endfor;
print;
print "Inverse entries written to file must be scaled by";
print "exp(";; 
format /re 10,7; print -aa;;
print ") to recover correct inverse ML";
end;