PROGRAM MLIKE
      PARAMETER(LD=100,LDK=100)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4 CPSEC
      PARAMETER(LDSET=9)
      COMMON/AAML/WT,PRI,PDAT,PAR(250),ITER,NPAR,LWRITE
      COMMON/ATWSIM/LSET
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      LOGICAL LALL,LDUMMY,LRANGE,LWRITE
      DIMENSION THETA(LDK),THETA1(LDK),THETA2(LDK),TMU(LDK)
      DIMENSION TVAR(LDK,LDK),TPRE(LDK,LDK),WPP(3)
      DIMENSION ADJ(LDSET),CUT(LDSET),DEN(LDSET),G(LDSET),TOP(LDSET)
      LSET=LDSET
      CALL NOWPR(6)
      TIME=CPSEC()
      CALL INIT
      READ(5,*)NBURN,NWRITE
      LWRITE=.TRUE.
      IF(NWRITE.EQ.0)LWRITE=.FALSE.
      CALL RDSIM0(1,250,NITER,NPAR)
      IF(LWRITE)CALL FOPEN('Marginal likelihood simulations',2,.TRUE.)
      IF(NBURN.GT.NITER)GO TO 1030
      NRANGE=NITER
      CALL REPAR0
      NUSE=NITER-NBURN
      WRITE(6,5)NITER,NBURN,NUSE
    5 FORMAT(/,' Marginal likelihood analysis',/,
     1  '     Total iterations:',I10,/,   
     2  '   Burn-in iterations:',I10,/,
     3  '  Analysis iterations:',I10)
C  Section 1: Compute mean and variance of Theta, and constants.
      WTSUM=0.0D0
      CALL DSET(NPAR,0.0D0,TMU,1)
      CALL UM0SET(NPAR,NPAR,TVAR,LDK)
      BASE=0.0D0
C  Section 1a: Pass through file to cumulate sums.
      DO 15 ITER=1,NITER
      READ(1,*,END=1010,ERR=1010)JTER,WTLOG,PRI,PDAT
      READ(1,*,END=1010,ERR=1010)(PAR(IPAR),IPAR=1,NPAR)
      IF(ITER.LE.NBURN)GO TO 15
      CALL REPAR(THETA,NK)
      BASE=BASE+PRI+PDAT+0.5D0*NK
      WT=DEXP(WTLOG)
      WTSUM=WTSUM+WT
      DO 10 IK=1,NK
      AA=WT*THETA(IK)
      TMU(IK)=TMU(IK)+AA
      DO 10 JK=IK,NK
   10 TVAR(IK,JK)=TVAR(IK,JK)+AA*THETA(JK)
   15 CONTINUE
C  Section 1b:  Convert sums to moments.
      BASE=BASE/NUSE
      CALL DSCAL(NK,1.0D0/WTSUM,TMU,1)
      DO 20 IK=1,NK
      DO 20 JK=IK,NK
   20 TVAR(IK,JK)=TVAR(IK,JK)/WTSUM-TMU(IK)*TMU(JK)
      CALL DFULL(NK,TVAR,LDK)
      CALL DPDCK(TVAR,LDK,NK,INFO)
      IF(INFO.NE.0)GO TO 1020
      CALL DLINDS(NK,TVAR,LDK,TPRE,LDK)
C  Section 2: Prepare normal density constants for this problem.
  100 SET1=LSET+1
      DK=NK
C  Section 2a: Set p = .9, .8, ... , .1.
      DO 110 ISET=1,LSET
      P=(LSET+1-ISET)/SET1
      ADJ(ISET)=1.0D0/P
  110 CUT(ISET)=.5D0*DCHIIN(P,DK)
      LDUMMY=LRANGE(NK,THETA,LALL)
      IF(LALL)GO TO 200
C  Section 2b:  Integrate multivariate nomral over parameter space (big Theta).
      CALL GAU(NK,TMU,TPRE,LDK,THETA1,THETA2)
      CALL DSET(LSET,0.0D0,DEN,1)
      CALL DSET(LSET,0.0D0,TOP,1)
      DO 120 IRANGE=1,NRANGE
      CALL GAUA(NK,TMU,TPRE,LDK,THETA1,THETA2) 
      SIZE=-GAUK(NK,TMU,TPRE,LDK,THETA1)
      ADD=0.0D0
      IF(LRANGE(NK,THETA1,LALL))ADD=0.5D0
      IF(LRANGE(NK,THETA2,LALL))ADD=ADD+0.5D0
      DO 120 ISET=1,LSET
      IF(SIZE.GT.CUT(ISET))GO TO 120
      DEN(ISET)=DEN(ISET)+1.0D0
      TOP(ISET)=TOP(ISET)+ADD
  120 CONTINUE
      DO 130 ISET=1,LSET
      IF(TOP(ISET).LE.0.0D0)TOP(ISET)=-1.0D0
  130 ADJ(ISET)=ADJ(ISET)*DEN(ISET)/TOP(ISET)
  140 WRITE(6,150)(ADJ(ISET),ISET=1,LSET)
  150 FORMAT(/,' Adjustment factors for coverage',/,(5(1PE16.7)))
C  Section 3: Construct Gelfand-Dey posterior moment and print.
  200 REWIND(1)
      READ(1,*)NITER,NPAR
      CALL DSET(LSET,0.0D0,TOP,1)
      CALL DSET(LSET,0.0D0,DEN,1)
      IF(LWRITE)WRITE(2,205)NUSE,LSET
  205 FORMAT(2I16)
      WTSUM=0.0D0
C  Section 3a: Pass through file to cumulate sums.
      DO 230 ITER=1,NITER
      CALL RDSIM(1,JTER,WPP,PAR)
      WTLOG=WPP(1)
      PRI=WPP(2)
      PDAT=WPP(3)
      IF(ITER.LE.NBURN)GO TO 230
      CALL REPAR(THETA,NK)
      WT=DEXP(WTLOG)
      WTSUM=WTSUM+WT
      SIZE=-GAUK(NK,TMU,TPRE,LDK,THETA)
      RATIO=DEXP(BASE-SIZE-PRI-PDAT)
C  Loop for all LSET alternative values of p (to 210).
      DO 210 ISET=1,LSET
      G(ISET)=0.0D0
      IF(SIZE.GT.CUT(ISET))GO TO 210
      G(ISET)=ADJ(ISET)*RATIO
      TOP(ISET)=TOP(ISET)+WT*G(ISET)
      DEN(ISET)=DEN(ISET)+WT
  210 CONTINUE
      IF(LWRITE)CALL WTSIM(2,ITER,WPP,G)
  220 FORMAT(I16,1PE16.7,/,(5(1PE16.7)))
  230 CONTINUE
C  Section 3b:  Print marginal likelihood.
      AA=BASE-GAUN(NK,TPRE,LDK)
      WRITE(6,240)(ISET,(LSET+1-ISET)/SET1,DEN(ISET)/WTSUM,
     1  DLOG(WTSUM/TOP(ISET))+AA,ISET=1,LSET)
  240 FORMAT(/,10X,'Marginal Likelihood (ML) Computations',//,
     1  ' Version    Prob',8X,'Weighted fraction included',6X,
     2  'Log ML',/,(I5,0PF12.6,0PF23.6,1PE24.7))
      IF(LWRITE)WRITE(6,250)-AA
  250 FORMAT(/,' Inverse ML entries written to file',
     1  ' must be scaled by',/,' exp(',1PE16.7,
     2  ') to recover correct inverse ML.')
      TIME=CPSEC()
      WRITE(6,260)TIME
  260 FORMAT(/,'   CPU time:',F10.2)
      CALL NOWPR(6)
      STOP
 1010 WRITE(6,1015) JTER
 1015 FORMAT(' *** ERROR IN READING INPUT FILE ITERATION',I8)
      CALL TERM
 1020 WRITE(6,1025)
 1025 FORMAT(/,' SINGULAR PARAMETER VECTOR VARIANCE')
      CALL TERM
 1030 WRITE(6,1035)NBURN,NITER
 1035 FORMAT(/,' ***',I10,
     1  ' burn in requested but only',I10,' available.')
      CALL TERM
      END