! Last change: DOS 3 Aug 2000 5:17 pm ! *** copyright 2000 *** ! *** filename chex132.f95 *** John F. Monahan ** ! ********************** program chex132 ! chex132.f95 -- Ramus height (jaw) -- Example 13.2 ! compute posterior distribution using Gibbs sampling ! implicit none integer, parameter :: n = 20 ! sample size integer, parameter :: p = 4 ! dimension integer, parameter :: nrep = 4000 ! MC replication size integer, parameter :: q = 2 ! dimension of beta real, dimension(q) :: beta ! parameters real, dimension((p*(p+1))/2) :: omega,e,estst,sigma real, dimension(p) :: ybar,ybmm,uc real, dimension(p,p) :: omhalf real, dimension(q) :: bw,z,betam real, dimension((q*(q+1))/2) :: xox,bcov real det,s,s1,s2,s3,fk real ran,gchirv,gnroul ! integer i,j,k,l,krep,idet ! 20 format(' Jaw Example -- Example 13.2 with',i6,' replications') 21 format(i4,2x,2f8.3,2x,5f8.3/10x,5f8.3,4x,3f8.3/12x,7f8.3) 22 format(' Mean beta(1)',f12.4,4x,'Cov Mx',f12.6/ & & 6x,'beta(2)',f12.4,10x,2f12.6) ! data ybar/ 48.66, 49.62, 50.57, 51.45/ data e/120.2695, 117.5925, 122.5375, 109.763, 116.915,& & 131.442, 105.415, 112.545, 131.98, 141.83 / ! starting values from mode of posterior data beta/ 33.818619, 1.855641 / ! output for summary open( unit=6, file='chex132.out' ) ! output for further analysis open( unit=8, file='chex132.dgn' ) ! write(6,20) nrep ! initializations: rng, mean, cov s1 = ran(5151917) betam = 0. bcov = 0. ! *** main loop *** do krep = 1,nrep fk = real(krep) ! first get differences from means do i = 1,p ybmm(i) = ybar(i) - beta(1) - beta(2)*(15. + real(i))/2. end do ! loop on i ! now get omega and generate Cholesky factor k = 0 do i = 1,p do j = 1,i k = k + 1 estst(k) = e(k) + real(n)*ybmm(i)*ybmm(j) if( i .eq. j ) estst(k) = estst(k) + 1. ! add prior matrix C end do ! loop on j end do ! loop on i call chlzky(estst,p,det,idet) ! generate Wishart factor do j = 1,p ! fill by column of cholesky factor of omega do i = 1,p if( i .lt. j ) uc(i) = 0. if( i .eq. j ) uc(i) = gchirv( real(n+p-i+1) ) if( i .gt. j ) uc(i) = gnroul(i) end do ! loop on i call chlzih(estst,p,uc) ! fill column do i = 1,p omhalf(i,j) = uc(i) end do ! loop on i end do ! loop on j ! get omega k = 0 do i = 1,p do j = 1,i k = k + 1 s = 0. do l = 1,p s = s + omhalf(i,l)*omhalf(j,l) end do ! loop on l omega(k) = s end do ! loop on j end do ! loop on i ! write out omega now that we've got it ! ! get x'omega x -- first x'omhalf in s1,s2 xox = 0. bw = 0. do j = 1,p s1 = 0. s2 = 0. s3 = 0. do k = 1,p s1 = s1 + omhalf(k,j) s2 = s2 + ((15.+real(k))/2.) * omhalf(k,j) s3 = s3 + ybar(k) * omhalf(k,j) end do ! loop on k xox(1) = xox(1) + s1*s1 xox(2) = xox(2) + s1*s2 xox(3) = xox(3) + s2*s2 bw(1) = bw(1) + s1*s3 bw(2) = bw(2) + s2*s3 end do ! loop on j ! factor x'omega x call chlzky(xox,q,det,idet) ! solve to get beta-wiggle call chlzhi(xox,q,bw) call chlzih(xox,q,bw) ! get normal( beta-wiggle, (n xox)inv ) xox = xox * sqrt(real(n)) ! rescale Cholesky factor z(1) = gnroul(1) z(2) = gnroul(2) call chlzih(xox,q,z) beta(1) = bw(1) + z(1) beta(2) = bw(2) + z(2) ! get sigma from omega sigma = omega call chlzky(sigma,p,det,idet) ! Cholesky factor call chlzoi(sigma,p) ! inverse of omega ! write out new beta, omega, sigma write(8,21) krep,beta,omega,sigma ! compute simple statistics betam = betam + beta if( krep .gt. 1 ) then bcov(1) = bcov(1) + & & (fk*beta(1)-betam(1))*(fk*beta(1)-betam(1))/(fk*(fk-1.)) bcov(2) = bcov(2) + & & (fk*beta(1)-betam(1))*(fk*beta(2)-betam(2))/(fk*(fk-1.)) bcov(3) = bcov(3) + & & (fk*beta(2)-betam(2))*(fk*beta(2)-betam(2))/(fk*(fk-1.)) end if ! ( krep .gt. 1 ) end do ! loop on krep fk = real(nrep) betam = betam/fk bcov = bcov / (fk-1.) write(6,22) betam(1),bcov(1),betam(2),bcov(2),bcov(3) stop end program chex132 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 SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL, 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, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL 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 CHLZHI(A,N,Y) ! ! 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 QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE ! MATRIX A, THAT IS Y' 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 THE SUM OF SQUARES OF THE ELEMENTS OF 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, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,J,IM1 ! ! SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN ! DO I = 1,N ! BRANCH AROUND ON THE FIRST ONE IF( I .GT. 1 ) THEN IM1 = I - 1 DO J = 1,IM1 Y(I) = Y(I) - A( (I*(I-1))/2 + J )*Y(J) ENDDO ENDIF Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZHI 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, DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL, 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 SUBROUTINE CHLZOI(A,N) ! ! OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER ! TRIANGULAR (N BY N) MATRIX STORED IN A ! ! TO COMPUTE THE INVERSE OF A POSITIVE DEFINITE MATRIX A ! FIRST: FACTOR A = L * L' WHERE L IS LOWER TRIANGULAR ! USING CHLSKY(A,N,D,IDET) ! SECOND: COMPUTE THE INVERSE(A) = INVERSE(L') * INVERSE(L) ! USING CHLSOI(A,N) ! ! J F MONAHAN (JULY,1984) DEPT OF STATISTICS, N C S U, RALEIGH, N C 2769 ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL S INTEGER I,J,K,IM1,KP1 ! FIRST INVERT L IN PLACE DO K = 1,N ! SOLVE L * X = E(K) A( (K*(K+1))/2 ) = 1./A( (K*(K+1))/2 ) ! BRANCH AROUND IF K EQUALS N IF( K .LT. N ) THEN KP1 = K + 1 DO I = KP1,N S = 0. IM1 = I - 1 DO J = K,IM1 S = S - A( (I*(I-1))/2 + J )*A( (J*(J-1))/2 + K ) ENDDO ! LOOP ON J A( (I*(I-1))/2 + K ) = S/A( (I*(I+1))/2 ) ENDDO ! LOOP ON I ENDIF ENDDO ! LOOP ON K ! DONE WITH INVERSION, NOW COMPUTE INNER PRODUCT DO J = 1,N DO I = J,N ! GET (I,J) ELEMENT OF MATRIX IS THE INNER PRODUCT ! OF COLUMN I AND COLUMN J S = 0. DO K = I,N S = S + A( (K*(K-1))/2 + I )*A( (K*(K-1))/2 + J ) ENDDO ! LOOP ON K ! STORE THE INNER PRODUCT A( (I*(I-1))/2 + J ) = S ENDDO ! LOOP ON I ENDDO ! LOOP ON J RETURN END SUBROUTINE CHLZOI ! *** end of filename chex132.f95 *********************