! Last change: DOS 31 Jul 2000 8:08 pm ! *** copyright 2000 *** ! *** filename quad3.f95 *** John F. Monahan ** ! ********************** program quad3 ! demonstration of importance sampling and antithetic variates ! ! compute posterior mean & variance for Example 10.1 ! implicit none integer, parameter :: np = 100 integer, parameter :: ntrial = 10 real, dimension(ntrial) :: pmnimp,pmnant real t,s0,s1,s2,tm,se,z,wgt real pstar,ran,gnroul integer i,ktrial ! 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(' with Antithetic Var ',e12.4,2f16.8) 23 format(/' Analysis of computed posterior means'/26x,'mean', & & 13x,'std error') 24 format(/' Trial',i3,' of',i3) ! ! put output here open( unit=6, file='quad3.out' ) ! ! set seed for uniform pseudorandom generator s0 = ran(5151917) ! mode of posterior for normal approx tm = .527167 ! standard error from second deriv at mode of post se = 1./sqrt(42.036) ! ! ! do 10 trials of 100 evaluations each write(6,20) np ! do ktrial = 1,ntrial write(6,24) ktrial,ntrial ! do the straight importance sampling first s0 = 0. s1 = 0. s2 = 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.) s0 = s0 + wgt s1 = s1 + wgt*t s2 = s2 + wgt*t*t end do ! loop on i ! 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 pmnimp(ktrial) = s1 ! save posterior mean for later analysis ! now include antithetic variates s0 = 0. s1 = 0. s2 = 0. do i = 1,np,2 ! 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.) s0 = s0 + wgt s1 = s1 + wgt*t s2 = s2 + wgt*t*t ! redo now for -z ! 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.) s0 = s0 + wgt s1 = s1 + wgt*t s2 = s2 + wgt*t*t end do ! loop on i ! 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,22) s0,s1,s2 pmnant(ktrial) = s1 ! save posterior mean for later analysis ! end do ! loop on ktrial ! analyze computed posterior mean ! get mean and std error ! first for importance sampling write(6,23) s1 = 0. s2 = 0. do i = 1,ntrial t = (pmnimp(i)-s1)/real(i) s2 = s2 + real(i-1)*(pmnimp(i)-s1)*t s1 = s1 + t end do ! loop on i s2 = sqrt( s2 / real(ntrial*(ntrial-1)) ) write(6,21) s1,s2 ! second for antithetic variates s1 = 0. s2 = 0. do i = 1,ntrial t = (pmnant(i)-s1)/real(i) s2 = s2 + real(i-1)*(pmnant(i)-s1)*t s1 = s1 + t end do ! loop on i s2 = sqrt( s2 / real(ntrial*(ntrial-1)) ) write(6,22) s1,s2 stop end program quad3 real function pstar(t) ! posterior from logarithmic series example 10.1 implicit none real, intent(in) :: 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 quad3.f95 *********************