! Last change: DOS 2 Aug 2000 9:22 pm ! *** copyright 2000 *** ! *** filename mixed2.f95 *** John F. Monahan ** ! ********************** program mixed2 ! mixed2 -- integration of posterior from Example 10.3 ! a variance components problem -- see also chex122a, chex122b ! using mixed spherical-radial quadrature ! ! centered at posterior mode ! and scaled posterior normal approximation ! implicit none integer, parameter :: d = 3 ! number of parameters integer, parameter :: npr = 256 ! number of radial points integer, parameter :: q = 16 ! number of rotations integer, parameter :: km = 2*d ! number of pts on sphere integer, parameter :: nf = ((d+2)*(d+1))/2 ! number of functions real(KIND=8), parameter :: twopi = 6.2831853071795864d0 real(KIND=8), dimension(nf) :: s,sij,si,ss,ssi real(KIND=8), dimension(d) :: t,tm real(KIND=8), dimension( (d*(d+1))/2 ) :: h real(KIND=8), dimension(d,2*d) :: a real(KIND=8) rll,pstmx,f,r,u,det,fj real(KIND=8) rlgpst real ran integer i,j,k,l,idet ! 22 format(4x,2f16.8) 23 format(f9.4,4f12.8) 24 format(f9.4,2f16.8) 25 format(i12,f12.8) 27 format(/'Relative standard error in normalization constant',f8.4) 28 format(/'Mixed Spherical-Radial',i6,' points -- norm constant', & & e12.6/10x,'Posterior mean',10x,'Covariance Matrix') 29 format(' theta(1)',f12.6,6x,f12.6/ & & ' theta(2)',f12.6,6x,2f12.6/ & & ' theta(3)',f12.6,6x,3f12.6) ! ! INTERFACE BLOCK INTERFACE subroutine spkblh(p,a) integer, intent(in) :: p real(KIND=8), dimension(:,:), intent(out) :: a end subroutine spkblh END INTERFACE ! ! mode of posterior data tm/1.200431d0, 0.905818d0, 7.327349d0 / ! hessian matrix from normal approximation data h/ 4.733423d0,1.007230d0, 6.689120d0,.494502d0, -.701696d0,& & 3.559180d0 / ! max of log posterior data pstmx/ 18.466364d0 / ! ! ! output unit for most results open( unit=6, file='mixed2.out' ) ! output unit for diagnostics open( unit=8, file='mixed2.dgn' ) ! ! initialize random number generator r = ran(5151917) ! find Cholesky factor of hessian call chlzky(h,d,det,idet) ! write(6,22) h write(6,25) idet,det ! r = 0. f = 1. rll = 0. write(8,24) r,f,rll ! ! initialize overall sums ! first Simpson weight for mode u = 16.d0/real(3*npr,8) !KIND s(1) = s(1) + u s(2) = s(2) + u*tm(1) s(3) = s(3) + u*tm(2) s(4) = s(4) + u*tm(3) s(5) = s(5) + u*tm(1)*tm(1) s(6) = s(6) + u*tm(1)*tm(2) s(7) = s(7) + u*tm(2)*tm(2) s(8) = s(8) + u*tm(1)*tm(3) s(9) = s(9) + u*tm(2)*tm(3) s(10) = s(10) + u*tm(3)*tm(3) ! ss = 0. ! ! npr = number of radial points do i = 1,npr ! get radius r = 16.d0*real(i,8)/real(npr,8) !KIND ! Simpson's rule weights u = 4.d0*16.d0/real(3*npr,8) !KIND if( 2*(i/2) .eq. i ) u = 2.d0*16.d0/real(3*npr,8) !KIND if( i .eq. npr ) u = 16.d0/real(3*npr,8) !KIND ! ! initialize sums for this r si = 0. ssi = 0. ! ! q = number of rotations do j = 1,q fj = real(j,8) !KIND ! get a set of points call spkblh(d,a) ! ! initialize sum for this rotation sij = 0. ! km = number of points on sphere do k = 1,km ! first L-inv-transpose * z t(1) = a(1,k) t(2) = a(2,k) t(3) = a(3,k) call chlzih(h,3,t) ! now center at mode t(1) = r*t(1) + tm(1) t(2) = r*t(2) + tm(2) t(3) = r*t(3) + tm(3) ! evaluate log posterior - max rll = rlgpst(t) - pstmx f = 0. if( rll .lt. 30.d0 ) f = exp(-rll) ! ! add up contribution from this point sij(1) = sij(1) + f sij(2) = sij(2) + f*t(1) sij(3) = sij(3) + f*t(2) sij(4) = sij(4) + f*t(3) sij(5) = sij(5) + f*t(1)*t(1) sij(6) = sij(6) + f*t(1)*t(2) sij(7) = sij(7) + f*t(2)*t(2) sij(8) = sij(8) + f*t(1)*t(3) sij(9) = sij(9) + f*t(2)*t(3) sij(10) = sij(10) + f*t(3)*t(3) ! end do ! loop on k ! ! get the means and variances over rotations j do l = 1,nf sij(l) = sij(l) / real(km,8) !KIND si(l) = si(l) + sij(l) if( j .ne. 1 ) ssi(l) = ssi(l) + & & (fj*sij(l)-si(l))*(fj*sij(l)-si(l))/(fj*(fj-1.d0)) end do ! loop on l ! 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,8) !KIND ssi(l) = ssi(l) / real( q*(q-1),8 ) !KIND 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 ! ! get relative error f = sqrt( ss(1) ) / s(1) write(6,27) f ! posterior moments s(2) = s(2) / s(1) s(3) = s(3) / s(1) s(4) = s(4) / s(1) s(5) = s(5) / s(1) - s(2)*s(2) s(6) = s(6) / s(1) - s(2)*s(3) s(7) = s(7) / s(1) - s(3)*s(3) s(8) = s(8) / s(1) - s(2)*s(4) s(9) = s(9) / s(1) - s(3)*s(4) s(10) = s(10) / s(1) - s(4)*s(4) ! ! fix normalization constant s(1) = 2.* twopi * s(1) * exp(-pstmx) / (det*2.**idet) write(6,28) npr*q*km,s(1) write(6,29) s(2),s(5),s(3),s(6),s(7),s(4),(s(l),l=8,10) ! done stop end program mixed2 REAL(KIND=8) function rlgpst(t) ! - log( posterior ) for exchange variance components problem ! data are 100*(rate-1.6) for USD/GBP rate in May 1997 ! ! t1 error/within variance component ! t2 trader/tmt variance component ! t3 mean parameter ! implicit none real(KIND=8), dimension(3), intent(in) :: t integer, parameter :: p = 5 ! number of groups integer, dimension(p) :: ni real(KIND=8), dimension(p) :: yb real(KIND=8) w,a1,b1,a2,b2,b0,phi0,sf,sg ! integer i,nbig ! data for problem ! p number of groups ! ni number of obs in group i ! nbig sum of ni ! w within/error sum of squares ! a1,b1 prior for error component is inversegamma( a1, b1 ) ! a2,b2 prior for treatment component is inversegamma( a2, b2 ) ! b0,phi0 prior for mean in Normal( b0, phi0 ) ! yb mean for group i data ni/ 2, 1, 7, 3, 2 / data yb/ 6.8d0, 10.3d0, 5.9d0, 6.6d0, 8.5d0/ ! w = 13.10d0 ! within ss nbig = 15 a1 = 1. b1 = 1. a2 = 4. b2 = 4. b0 = 8. phi0 = 16. ! assign outlanding posterior (tiny) if out of range rlgpst = 1.d6 if( t(1) .le. 0.d0 ) return if( t(2) .le. 0.d0 ) return ! first part of posterior rlgpst = (a2 + p/2. + 1)*log(t(2)) + (a1+nbig/2.+1.)*log(t(1)) sg = 0. sf = 0. do i = 1,p sg = sg + log( ni(i)/t(1) + 1./t(2) ) sf = sf + (t(3)-yb(i))*(t(3)-yb(i))/(t(2)+t(1)/ni(i)) end do ! loop on i rlgpst = rlgpst + b2/t(2) + (b1+w/2.)/t(1) & & + (t(3)-b0)*(t(3)-b0)/(2.*phi0) & & + sg/2. + sf/2. ! ! using - log( posterior) *** reminder rlgpst = rlgpst return end function rlgpst subroutine spkblh(p,a) ! compute random rotation of antipodal points on surface of a sphere ! in p dimensions -- gives 2*p points stored as columns in a ! ! p integer dimension of sphere ! a real points on surface: point k is in a( ,k) ! ! J F Monahan (June 1994) Dept of Statistics, N C State U, Raleigh 27695 ! recoded June 1998 ! Recoded February, April 2000 for Fortran 95 implicit none integer, intent(in) :: p real(KIND=8), dimension(:,:), intent(out) :: a ! dummy dimensions real(KIND=8), dimension(p) :: x integer i,ii,j,k,p2 real(KIND=8) s,d real gnroul ! compute some constants p2 = p + p ! initialize to zero first a = 0. ! now put 1 and -1 for poles do j = 1,p a(j,j) = 1. a(j,j+p) = -1. end do ! loop on j ! matrix a is now I(p), -I(p) do ii = 2,p i = p - ii + 1 ! generate random rotation by Householder transforms ! defined by iid normal elements s = 0. do j = i,p x(j) = gnroul(i+j) s = s + x(j)*x(j) end do ! loop on j d = s s = sqrt(s) d = d - s*x(i) ! use negative to insure R x(i) = x(i) - s ! has positive diagonals ! x now holds vector u, d holds constant ! ! apply rotation to axial points do k = 1,p2 ! Householder requires one inner product s = 0. do j = i,p s = s + x(j)*a(j,k) end do ! loop on j s = s / d do j = i,p a(j,k) = a(j,k) - s*x(j) end do ! loop on j ! done with this rotation end do ! loop on j end do ! loop on ii return end subroutine spkblh 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 SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL(KIND=8), INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL(KIND=8), INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL(KIND=8) T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A( (I*(I+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + J ) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A( (I*(I+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY SUBROUTINE CHLZIH(A,N,Y) ! *** NOTICE TRANSPOSE *** ! OVERWRITES Y WITH INVERSE( L' ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX ! AND A VECTOR Y, THAT IS INVERSE( A ) * Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y) ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,II,IP1,J ! ! SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP ! ! DO THE LAST ONE SEPARATELY Y(N) = Y(N) / A( (N*(N+1))/2 ) ! IF N = 1 THEN QUIT ELSE START LOOPING IF( N .EQ. 1 ) RETURN DO II = 2,N ! I BELOW NOW GOES FROM N-1 DOWN TO 1 I = N + 1 - II IP1 = I + 1 DO J = IP1,N Y(I) = Y(I) - A( (J*(J-1))/2 + I )*Y(J) ENDDO ! LOOP ON J Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON II,I RETURN END SUBROUTINE CHLZIH ! *** end of filename mixed2.f95 *********************