PROGRAM ROBUST
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (LD=50000,LDPAR=300,LVAL=100)
      DIMENSION G(LD),WT(LD)
      DIMENSION PAR(LDPAR),WPP(5),BETA(LVAL),DPAR(LVAL),
     1  QDOWN(LVAL),QUP(LVAL)
      COMMON/ROBBLK/AK
      EXTERNAL ROBEQ
      CALL INIT
C  Read in and check user information
      WRITE(6,10)
   10 FORMAT(10X,'Prior density ratio class robustness',/)
      READ(5,*)ND,NBURN,IFUNC
      CALL PINTI('Number of d values',ND,1,LVAL)
      CALL PLISTR(5,ND,DPAR)
      WRITE(6,20)NBURN,IFUNC
   20 FORMAT('  Number of burn-in iterations:',I10,/,
     1       '          Function of interest:',I10)
      CALL RDSIM0(1,LDPAR,NITER,NPARS)
      CALL PINTI('Burn-in iterations',NBURN,0,NITER)
      CALL PINTI('Function of interest',IFUNC,1,NPARS)
      NUSE=NITER-NBURN
      NSKIP=(NUSE-1)/LD+1
      IF(NBURN.EQ.0)GO TO 40
      DO 30 IBURN=1,NBURN
   30 CALL RDSIM(1,JTER,WPP,PAR)
C  Pass through file to get weights and values
   40 M=0
      DO 50 ITER=1,NUSE,NSKIP
      M=M+1
      CALL RDSIM(1,JTER,WPP,PAR)
      WTLOG=WPP(1)
      IF(ITER.EQ.1)WTMAX=WTLOG
      IF(WTLOG.GT.WTMAX)WTMAX=WTLOG
      WT(M)=WTLOG
   50 G(M)=PAR(IFUNC)
C  Convert from log weights to weights and find moments
      MTOT=M
      DO 60 M=1,MTOT
   60 WT(M)=DEXP(WT(M)-WTMAX)
      CALL ROB1(MTOT,WT,G,ND,DPAR,QUP,AMU,ASD)
      CALL DSCAL(MTOT,-1.0D0,G,1)
      CALL ROB1(MTOT,WT,G,ND,DPAR,QDOWN,AA,ASD)
      CALL DSCAL(ND,-1.0D0,QDOWN,1)
      DO 70 ID=1,ND
      AK=DEXP(DPAR(ID))
C      CALL DZREAL(ROBEQ,1.0D-6,1.0D-6,1.0D-6,1.0D-6,1,1000,0.0D0,F,INFO)
      INFO=100
      A1=0D0
      F=3D0
      CALL DZBREN(ROBEQ,1.0D-6,1.0D-6,A1,F,INFO)
   70 BETA(ID)=(AK-1.0D0)*PHI(F)/(AK+(1.0D0-AK)*DNORDF(F))
   75 format(6(1pe14.5))
      WRITE(6,80)IFUNC,AMU,ASD,(DPAR(ID),QDOWN(ID),QUP(ID),
     1             AMU-ASD*BETA(ID),AMU+ASD*BETA(ID),ID=1,ND)
   80 FORMAT(15X,'Robustness analysis of function',I3,//,
     1  ' Posterior mean',1PE15.6,'   Stan dev',1PE15.6,//,24X,'Exact',
     2  '            DeRobertis & Hartigan asymptotic' ,/,
     3  9X,'d    Lower bound    Upper bound',
     4  '    Lower bound    Upper bound',/,(0PF10.6,4(1PE15.6))) 
      END
C
      REAL*8 FUNCTION PHI(X)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/CONST/DLOG2,DLG2PI,HLG2PI
      PHI=DEXP(-HLG2PI-.5D0*X**2)
      RETURN
      END
C
      REAL*8 FUNCTION ROBEQ(F)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/ROBBLK/AK
      ROBEQ=(AK-1.0D0)*(PHI(F)+F*DNORDF(F))-AK*F
      RETURN
      END
C    
      SUBROUTINE ROB1(MTOT,WT,G,ND,DPAR,Q,AMU,ASD)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=50000)
      DIMENSION WT(MTOT),G(MTOT),DPAR(ND),Q(ND)
      DIMENSION D(LD),T(LD),NV1(LD)
      DO 5 M=1,MTOT
    5 NV1(M)=M
      CALL DSVRGP(MTOT,G,G,NV1)
      W1=0.0D0
      WG1=0.0D0
      WG2=0.0D0 
      DO 10 M=1,MTOT
      W=WT(NV1(M))
      WG=W*G(M)
      W1=W1+W
      WG1=WG1+WG
      WG2=WG2+WG*G(M)
      D(M)=W1
   10 T(M)=WG1
      AMU=WG1/W1
      ASD=DSQRT(WG2/W1-AMU**2)
      DO 50 ID=1,ND
      AK=DEXP(DPAR(ID))
      L=MTOT/2
      LSTEP=MTOT/2
   20 LSTEP=(LSTEP-1)/2+1
      QL=(T(L)+AK*(WG1-T(L)))/(D(L)+AK*(W1-D(L)))
      IF(QL.GT.G(L+1))GO TO 30
      IF(QL.LT.G(L))GO TO 40
      Q(ID)=QL
      GO TO 50
   30 L=L+LSTEP
      GO TO 20
   40 L=L-LSTEP
      GO TO 20
   50 CONTINUE
      RETURN
      END