! Last change: DOS 3 Aug 2000 5:19 pm ! *** copyright 2000 *** ! *** filename chex133.f95 *** John F. Monahan ** ! ********************** program chex133 ! Example 13.3 -- Metropolis within Gibbs for pump example ! implicit none integer, parameter :: nobs = 1000 ! number of observations integer, parameter :: p = 10 ! number of pumps (J) integer, parameter :: pp2 = p+2 ! total number of parameters integer, dimension(p) :: y ! failures real, DIMENSION(p) :: x,f real, DIMENSION(pp2) :: th,thmean real, DIMENSION(pp2,pp2) :: thcov real a,b real c,ft,anew,alpha,sumf,slf,v,lpstrold,lpstrnew real ran,gchirv,algama ! integer i,j,t,kaccept ! 21 format(i6,12f6.3) 28 format(/'Metropolis within Gibbs -- Example 13.3'/ & & 10x,'with',i8,' observations and',i8,' acceptances',& & //10x,'Posterior mean',10x,'Covariance Matrix') 29 format(' alpha',f12.6,6x,f12.6/ & & ' beta ',f12.6,6x,2f12.6) ! ! data for problem data x/ 94.3, 15.7, 62.9, 126., 5.24, 31.4, 1.05, 1.05, 2.1, 10.5/ data y/ 5, 1, 5, 14, 3, 19, 1, 1, 4, 22 / ! ! output unit for most results open( unit=6, file='chex133.out' ) ! output unit for diagnostics open( unit=8, file='chex133.dgn' ) ! ! prior parameter c = .1 ! initialize v = ran(5151917) thmean = 0. thcov = 0. kaccept = 0 ! starting values a = .7 b = .9 ! ! main iteration loop do t = 1,nobs ft = real(t) ! ! first the distribution of phi's given others slf = 0. sumf = 0. do j = 1,p v = gchirv( 2.*(real(y(j))+a) ) ! v is chi( 2(y(j)+a) ) f(j) = (v*v/2.) / (x(j)+b) ! phi(j) is gamma(y(j)+a,x(j)+b) slf = slf + log(f(j)) sumf = sumf + f(j) th(j+2) = f(j) ! copy into param vector end do ! loop on j ! generate new beta ! ! *** gchirv requires argument greater than 1 *** not guaranteed *** ! should crash if fails due to sqrt(negative) v = gchirv( 2.*(real(p)*a+c) ) ! v is chi( 2(Ja+c) ) b = (v*v/2.) / (1.+sumf) ! beta is gamma(Ja+c,1+sumf) th(2) = b ! copy into param vector ! ! last the distribution of alpha given others ! ! pstarold changes each time lpstrold = a*p*log(b) + a*slf - a - p*algama(a) anew = a + (ran(t) -.5) ! random walk to get candidate lpstrnew = -1.e10 if( anew .gt. 0. ) lpstrnew = & & anew*p*log(b) + anew*slf - anew - p*algama(anew) alpha = 1. if( lpstrnew .lt. lpstrold ) alpha = exp( lpstrnew - lpstrold ) ! got alpha; ready for trial if( ran(t) .lt. alpha ) then ! acceptance kaccept = kaccept + 1 a = anew end if ! ( ran(t) .lt. alpha ) th(1) = a ! copy into param vector ! ! write out results to file for later analysis write(8,21) t,th ! ! compute statistics do i = 1,pp2 thmean(i) = thmean(i) + th(i) if( t .gt. 1 ) then do j = 1,i thcov(i,j) = thcov(i,j) + & & (ft*th(i)-thmean(i))*(ft*th(j)-thmean(j))/(ft*(ft-1.)) end do ! loop on j end if ! ( t .gt. 1 ) end do ! loop on i end do ! loop on t -- iterations ft = float(nobs) thmean = thmean / ft thcov = thcov / (ft-1.) write(6,21) nobs,thmean do i = 1,pp2 write(6,21) i,(thcov(i,j),j=1,i) end do ! loop on i write(6,28) nobs,kaccept write(6,29) (thmean(i),(thcov(i,j),j=1,i),i=1,2) stop end program chex133 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 GCHIRV(A) ! ! ALGORITHM TO GENERATE RANDOM VARIABLES WITH THE CHI DISTRIBUTION ! WITH A DEGREES OF FREEDOM, FOR A GREATER THAN OR EQUAL TO ONE ! ! J F MONAHAN (MAY, 1986) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA ! RECODED (FEB, 2000) FOR FORTRAN 95 ! IMPLICIT NONE REAL, INTENT(IN) :: A REAL U,V,Z,ZZ,RNUM,W,S,VMAX REAL, SAVE :: ALPHM1,BETA,VMIN,VDIF ! REAL, PARAMETER :: EMHLF = .6065307 ! EXP(-1/2) REAL, PARAMETER :: RSQRT2 = .7071068 ! 1/SQRT(2) REAL, PARAMETER :: EMHLF4 = .1516327 ! EXP(-1/2)/4 REAL, PARAMETER :: EQTRT2 = 2.568051 ! 2*EXP(1/4) REAL, PARAMETER :: C = 1.036961 ! 4*EXP(-1.35) ! REAL, SAVE :: ALPHA = 0. ! GIVE INITIAL VALUE & SAVE REAL RAN ! IS THIS ALPHA THE SAME AS THE LAST ONE? IF( A .NE. ALPHA ) THEN ! DO A LITTLE SETUP ALPHA = A ALPHM1 = ALPHA - 1. BETA = SQRT( ALPHM1 ) ! GET DIMENSIONS OF BOX VMAX = EMHLF * ( RSQRT2 + BETA )/( .5 + BETA ) VMIN = -BETA IF( BETA .GT. 0.483643 ) VMIN = EMHLF4/ALPHA - EMHLF VDIF = VMAX - VMIN END IF ! ( A .NE. ALPHA) ! START ( AND RESTART ) ALGORITHM HERE DO ! NOTE UNRESTRICTED DO U = RAN(1) V = VDIF*RAN(2) + VMIN Z = V / U GCHIRV = Z + BETA ! RETURN ON SUCCESS ! DO SOME QUICK REJECT CHECKS FIRST IF( Z .LE. - BETA ) CYCLE ! REJECT ZZ = Z * Z RNUM = 2.5 - ZZ IF( V .LT. 0. ) RNUM = RNUM + ZZ * Z / ( 3. * (Z + BETA ) ) ! DO QUICK INNER CHECK IF( U .LT. RNUM/EQTRT2 ) RETURN ! ACCEPT IF( ZZ .GT. C / U + 1.4 ) CYCLE ! REJECT ! ABOVE WAS KNUTH'S NORMAL OUTER CHECK W = 2. * LOG( U ) ! NOW THE REAL CHECK S = - ( ZZ / 2. + Z * BETA ) IF( BETA .GT. 0. ) S = ALPHM1*LOG(1.+Z/BETA) + S IF( W .LE. S ) RETURN ! ACCEPT END DO ! UNRESTRICTED DO RETURN END FUNCTION GCHIRV 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 REAL FUNCTION ALGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(5) :: P,Q ! COEFS FOR RATIONAL FUNCTION REAL, DIMENSION(3) :: R ! COEFS FOR STIRLING CORRECTION REAL XN,A,V REAL, PARAMETER :: LSQR2PI = .918938533 DATA P/ -.510046056E1, -.368309363E2, -.664664448E2, & & -.195692802E3, -.185072903E1 / DATA Q/ -.114321842E2, .368848969E2, .16269403E2, & & -.195692802E3, 1.0 / DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0. ) THEN ALGAMA = -1.0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8. ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = R(2) + ( R(3)/XN + R(1) )/XN ALGAMA = (X-0.5)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.0 XN = X DO WHILE( XN .LT. 2. ) V = V/XN XN = XN + 1. END DO ! WHILE( XN .LT. 2. ) DO WHILE( XN .GT. 3. ) XN = XN - 1. V = V * XN END DO ! WHILE( XN .GT. 3. ) A = XN-2.0 ! RATIONAL FUNCTION APPROXIMATION V = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) / & & (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5))))) ALGAMA = LOG(V) ENDIF RETURN END FUNCTION ALGAMA ! *** end of filename chex132.f95 *********************