!     Last change:  DOS   2 Aug 2000    9:26 pm
!                                       *** copyright 2000 ***
! *** filename quad4.f95                *** John F. Monahan **
!                                       **********************
      program quad4
!  demonstration of density estimation using importance sampling demo
!
!  compute posterior mean & variance for Example 10.1
!
      implicit none
      real t,s0,s1,s2,tm,se,z,wgt,mxw,sumw,sumww,sw,sz,sww,szz,szw
      real ran,gnroul,pstar
      integer i,np
!
  20  format(/'  Number of evaluations',i5/23x,'Normalization',5x,    &
     &       'Posterior',7x,'Posterior'/26x,'constant',9x,'mean',     &
     &       11x,'variance')
  21  format(' Importance Sampling  ',e12.4,2f16.8)
  22  format(f12.6,f14.9)
  25  format(/' Two Useful statistics, S1',f12.6,6x,'S2',f12.6)
  26  format(/8x,' mean vector',10x,'covariance matrix'/4x,'z',f14.4,    &
     &   6x,f14.6/4x,'w',f14.4,6x,2f14.6)
!
!                   put output here
      open( unit=6, file='quad4.out' )
!                   put theta's and weights  here
      open( unit=8, file='quad4.wgt' )
!
!                   mode of posterior for normal approx
      tm = .527167
!                   standard error from second deriv at mode of post
      se = 1./sqrt(42.036)
!
!                   set seed for uniform pseudorandom generator
      s0 = ran(1917515)
!
!                   do a single trial of 400 evaluations 
      np = 100
      write(6,20) np
!                   do the straight importance sampling       
      s0 = 0.
      s1 = 0.
      s2 = 0.
      mxw = 0.
      sumww = 0.
      sw = 0.
      sz = 0.
      sww = 0.
      szw = 0.
      szz = 0.
      do i = 1,np
!                   generate standard normal random variable
      z = gnroul(i)
!                   relocate and rescale by mode and 2nd deriv      
      t = tm + se*z
!                   divide by normal density
      wgt = pstar(t)*2.506628*exp(z*z/2.)
      write(8,22) t,wgt
!                   typical posterior calculations
      s0 = s0 + wgt
      s1 = s1 + wgt*t        
      s2 = s2 + wgt*t*t
!                   some statistics on weighting
      mxw = max( mxw, wgt )
      sumww = sumww + wgt*wgt
!                   means, variances, etc on z and w
      wgt = 10000.*wgt        ! rescale by 10000 for easier reading
      z = t*wgt
      sz = sz + z
      sw = sw + wgt
      if( i .gt. 1 ) szz = szz + (sz-i*z)*(sz-i*z)/(i*(i-1))
      if( i .gt. 1 ) sww = sww + (sw-i*wgt)*(sw-i*wgt)/(i*(i-1))
      if( i .gt. 1 ) szw = szw + (sz-i*z)*(sw-i*wgt)/(i*(i-1))
      end do  ! loop on i
      sumw = s0
!                  now convert to get posterior means and variances
      s1 = s1 / s0
      s2 = s2/s0 - s1*s1
!                  rescale s0 by se
      s0 = s0 * se / real(np)
      write(6,21) s0,s1,s2
!
!                  two useful statistics
      s1 = mxw / sumw
      s2 = sumww / (sumw*sumw)
      write(6,25) s1,s2
!                  get mean vector, covariance of (Z,W)
      sz = sz / real(np)
      sw = sw / real(np)
      szz = szz / real(np-1)
      sww = sww / real(np-1)
      szw = szw / real(np-1)
      write(6,26) sz,szz,sw,szw,sww
      stop
      end program quad4
      real function pstar(t)
!  posterior from logarithmic series example 10.1      
      real t
      pstar = 0.
!                   first take care of values out of range
      if( (t .le. 0.) .or. (t .ge. 1.) ) return
!                   evaluate function
      pstar = (1.-t)*( t**16 )/ ( (-log(1.-t))**10 )
      return
      end function pstar
      REAL FUNCTION GNROUL(IX)
!  RATIO OF UNIFORMS ALGORITHM FOR GENERATING NORMAL(0,1) RV'S
!   USING LEVA'S QUADRATIC INNER AND OUTER BOUNDS
!
!  A J KINDERMAN AND J F MONAHAN (1977) "COMPUTER GENERATION OF RANDOM
!   VARIABLES USING THE RATIO OF UNIFORM DEVIATES," ACM TRANSACTIONS ON
!   MATHEMATICAL SOFTWARE, VOLUME 3, PP.257-260
!
!  J L LEVA (1992) "A FAST NORMAL RANDOM NUMBER GENERATOR," ACM
!   TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 18, PP. 449-453.
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IX         ! DUMMY ARGUMENT
      REAL, PARAMETER :: R = 1.7155277  ! FIRST CONSTANT R = SQRT(2/E)
!                   CONSTANTS S,T CENTER OF ELLIPSES
      REAL, PARAMETER :: S = .449871
      REAL, PARAMETER :: T = -.386595
!                   SHAPE PARAMETERS OF ELLIPSES
      REAL, PARAMETER :: A = .19600
      REAL, PARAMETER :: B = .25472
!                   RADII OF ELLIPSES
      REAL, PARAMETER :: R1 = .27597
      REAL, PARAMETER :: R2 = .27846
      REAL RAN
      REAL U,V,X,Y,Q
!                   START / RESTART HERE
      DO            ! NOTE UNRESTRICTED DO
!                   GENERATE (U,V) UNIFORMLY IN BOX
      U = RAN(1)
      V = R*(RAN(2) - 0.5)
      X = U - S
      Y = ABS(V) - T
!                   COMPUTE QUADRATIC PIECE
      Q = X*X + Y*(A*Y - B*X)
!                   COMPUTE DEVIATE BEFORE TESTS
      GNROUL = V / U
!                   INNER BOUND -- QUICK ACCEPT CHECK
      IF( Q .LE. R1 ) RETURN
!                   OUTER BOUND -- QUICK REJECT CHECK
      IF( A .GT. R2 ) CYCLE
!                   FINAL RATIO OF UNIFORMS CHECK
      IF( GNROUL*GNROUL .LE. -4.*LOG(U) ) RETURN
!                   REJECT -- START OVER
      END DO        ! END OF UNRESTRICTED DO
      END FUNCTION GNROUL
      REAL FUNCTION RAN(IDUM)
!  PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR
!    LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR
!     X(N+1) = MOD( 16807*X(N), 2**31-1 )
!
!  P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION
!   SPRINGER-VERLAG, NEW YORK.  PP. 200-202.
!  UPDATED VERSION OF
!  LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' 
!   ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138.
!
!  ARGUMENT 
!   IDUM    INTEGER  FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      REAL, PARAMETER :: RP = 4.656612875E-10     ! 1/P
      INTEGER, PARAMETER :: A = 16807             ! MULTIPLIER
      INTEGER, PARAMETER :: B = 127773            ! B = P / A
      INTEGER, PARAMETER :: C = 2836              ! C = P MOD A
      INTEGER, PARAMETER :: P = 2147483647        ! MODULUS 2**31 - 1
      INTEGER, SAVE :: IX = 0
      INTEGER K1
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( IX .EQ. 0) IX = IDUM
!                   WRITE NUMBER AS ALPHA*2**16 + BETA
      K1 = IX / B
      IX = A*( IX - K1*B) - K1*C
!                   ABOVE DOES A*IX MOD B -K1*C
      IF( IX .LT. 0 ) IX = IX + P 
!                   RP BELOW IS RECIPROCAL OF P
      RAN = REAL(IX)*RP
      RETURN
      END FUNCTION RAN
! *** end of filename quad4.f95         *********************