! Last change: DOS 2 Aug 2000 9:20 pm ! *** copyright 2000 *** ! *** filename mixed1.f95 *** John F. Monahan ** ! ********************** program mixed1 ! mixed1 -- integration of posterior from Chapter 9's Extended Example ! (modified) using mixed spherical-radial quadrature ! ! centered at posterior mode ! and scaled posterior normal approximation ! implicit none integer, parameter :: npr = 256 ! number of radial points integer, parameter :: q = 4 ! number of rotations integer, parameter :: km = 8 ! number of pts on sphere/circle integer, parameter :: nf = 6 ! number of functions real, parameter :: twopi = 6.2831853 real, dimension(nf) :: s,ss,sij,si,ssi real t1,t2,tm1,tm2,h11,h12,h22,rll,pstmx,r,u,f,cs,sn,shift real ran,rlglik integer i,j,k,l ! 22 format(4x,2f16.8) 24 format(f9.4,2f16.8) 27 format(/'Relative standard error in normalization constant',f8.4) 28 format(/'Mixed Spherical-Radial',i6,' points' & & /' norm constant',f12.8,'*exp(',f10.4,' )' & & /10x,'Posterior mean',10x,'Covariance Matrix') 29 format(' theta(1)',f12.6,6x,f12.6/ & & ' theta(2)',f12.6,6x,2f12.6) ! ! mode of posterior data tm1,tm2/ .265782, .110980 / ! hessian matrix from normal approximation data h11,h12,h22/ 215.1117, 59.4624, 472.0310/ ! max of log posterior data pstmx/ -28.0271 / ! ! output unit for most results open( unit=6, file='mixed1.out' ) ! output unit for diagnostics open( unit=8, file='mixed1.dgn' ) ! ! initialize random number generator r = ran(5151917) ! initialize sums s = 0. ! find Cholesky factor of hessian h11 = sqrt(h11) h12 = h12 / h11 h22 = sqrt( h22 - h12*h12 ) ! write(6,22) h11,h12,h22 ! r = 0. f = 1. rll = 0. write(8,24) r,f,rll ! ! initialize overall sums ! first Simpson weight for mode u = 8./real(6*npr) s(1) = s(1) + u s(2) = s(2) + u*tm1 s(3) = s(3) + u*tm2 s(4) = s(4) + u*tm1*tm1 s(5) = s(5) + u*tm1*tm2 s(6) = s(6) + u*tm2*tm2 ! ss = 0. ! ! npr = number of radial points do i = 1,npr ! get radius r = 8.*real(i)/real(npr) ! Simpson's rule weights u = 4.*8./real(6*npr) if( 2*(i/2) .eq. i ) u = 2.*8./real(6*npr) if( i .eq. npr ) u = 8./real(6*npr) ! ! initialize sums for this r si = 0. ssi = 0. ! ! q = number of rotations do j = 1,q ! here just a shift on circle shift = ran(j) ! ! initialize sum for this rotation sij = 0. ! km = number of points on sphere (circle) do k = 1,km cs = r*cos( twopi*(shift + real(k))/real(km) ) sn = r*sin( twopi*(shift + real(k))/real(km) ) ! first L-inv-transpose * z t2 = sn / h22 t1 = (cs - h12*t2) / h11 ! now center at mode t1 = t1 + tm1 t2 = t2 + tm2 ! evaluate log posterior - max rll = rlglik(t1,t2) - pstmx f = 0. if( rll .gt. -30. ) f = exp(rll) ! ! add up contribution from this point sij(1) = sij(1) + f sij(2) = sij(2) + f*t1 sij(3) = sij(3) + f*t2 sij(4) = sij(4) + f*t1*t1 sij(5) = sij(5) + f*t1*t2 sij(6) = sij(6) + f*t2*t2 ! end do ! loop on k (around circle) ! ! get the means and variances over rotations j do l = 1,nf sij(l) = sij(l) / real(km) si(l) = si(l) + sij(l) if( j .ne. 1 ) ssi(l) = ssi(l) + & & (real(j)*sij(l)-si(l))*(real(j)*sij(l)-si(l)) / real(j*(j-1)) end do ! loop on nf ! done with rotations end do ! loop on j ! add up contribution from this radius r(i) do l = 1,nf si(l) = si(l) / real(q) ssi(l) = ssi(l) / real( q*(q-1) ) s(l) = s(l) + si(l)*u*r ss(l) = ss(l) + ssi(l)*u*u*r*r end do ! loop on l f = sqrt( ssi(1) ) write(8,24) r,si(1),f ! done with this r end do ! loop on i ! ! posterior moments s(2) = s(2) / s(1) s(3) = s(3) / s(1) s(4) = s(4) / s(1) - s(2)*s(2) s(5) = s(5) / s(1) - s(2)*s(3) s(6) = s(6) / s(1) - s(3)*s(3) ! ! get relative error f = sqrt( ss(1) ) / s(1) write(6,27) f ! fix normalization constant s(1) = twopi * s(1) / (h11 * h22) write(6,28) npr*q*km,s(1),pstmx write(6,29) s(2),s(4),s(3),s(5),s(6) ! done stop end program mixed1 REAL FUNCTION RLGLIK(T1,T2) ! COMPUTE LOG-LIKELIHOOD FOR EXTENDED EXAMPLE ! *** modified from Example 9.4 by dropping sample size *** IMPLICIT NONE REAL, INTENT(IN) :: T1,T2 REAL P1,P2,P3,P4 INTEGER N1,N2,N3,N4 DATA N1,N2,N3,N4/ 1, 10, 4, 9/ ! **************** numbers are changed from Chapter 9 P1 = 2.*T1*T2 P2 = T1*(2. - T1 - 2.*T2) P3 = T2*(2. - T2 - 2.*T1) P4 = (1.-T1-T2)*(1.-T1-T2) ! CHECK IF OUT OF RANGE ! ASSIGN TINY LIKELIHOOD IF OUTSIDE TRIANGLE RLGLIK = -1.E30 IF( P1 .LE. 0. ) RETURN IF( P2 .LE. 0. ) RETURN IF( P3 .LE. 0. ) RETURN IF( P4 .LE. 0. ) RETURN ! RLGLIK = N1*LOG(P1) + N2*LOG(P2) + N3*LOG(P3) + N4*LOG(P4) ! RETURN END FUNCTION RLGLIK 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 mixed1.f95 *********************