! Last change: DOS 31 Jul 2000 7:54 pm ! *** copyright 2000 *** ! *** filename chex102t.f95 *** John F. Monahan ** ! ********************** program chex102t ! ext102 -- integration of posterior from Chapter 9's Extended Example ! using triangular rules ! implicit none real, dimension(3) :: vx,vy real, dimension(6) :: ss,sm real, dimension(1000) :: xa,ya,wgt real t1,t2,rll,pstmx,f real rlglik integer kint,i,np ! 27 format(/' Midpoint rule',i4,' points -- Normalization constant',& & f12.8/10x,'Posterior mean',10x,'Covariance Matrix') 28 format(/'Simpson''s rule',i4,' points -- Normalization constant',& & f12.8/10x,'Posterior mean',10x,'Covariance Matrix') 29 format(' theta(1)',f12.6,6x,f12.6/ & & ' theta(2)',f12.6,6x,2f12.6) ! ! vertices of triangular parameter space data vx/ 0., 1., 0./ data vy/ 0., 0., 1./ ! ! output unit open( unit=6, file='chex102t.out' ) ! ! max of log posterior pstmx = -28.0271 ! first the triangular integration do kint = 10,25,5 ! initialize to zero sm = 0. ss = 0. ! get triangle abscissas and weights call triqds(kint,vx,vy,xa,ya,wgt) ! how many points for Simpson np = (3*kint*(kint+1)+2)/2 do i = 1,np t1 = xa(i) t2 = ya(i) rll = rlglik(t1,t2) - pstmx if( rll .gt. -30. ) then f = exp(rll)*wgt(i) else f = 0. end if ! ( rll .ggt. -30. ) ! Simpson rule -- use all points ss(1) = ss(1) + f ss(2) = ss(2) + f*t1 ss(3) = ss(3) + f*t2 ss(4) = ss(4) + f*t1*t1 ss(5) = ss(5) + f*t1*t2 ss(6) = ss(6) + f*t2*t2 if( i .le. kint*kint ) then ! Midpoint rule -- use just first kint*kint points f = f*4./3. sm(1) = sm(1) + f sm(2) = sm(2) + f*t1 sm(3) = sm(3) + f*t2 sm(4) = sm(4) + f*t1*t1 sm(5) = sm(5) + f*t1*t2 sm(6) = sm(6) + f*t2*t2 end if ! ( i .le. kint*kint ) end do ! loop on i ! posterior moments -- midpoint first sm(2) = sm(2) / sm(1) sm(3) = sm(3) / sm(1) sm(4) = sm(4) / sm(1) - sm(2)*sm(2) sm(5) = sm(5) / sm(1) - sm(2)*sm(3) sm(6) = sm(6) / sm(1) - sm(3)*sm(3) write(6,27) kint*kint,sm(1) write(6,29) sm(2),sm(4),sm(3),sm(5),sm(6) ! posterior moments -- Simpson ss(2) = ss(2) / ss(1) ss(3) = ss(3) / ss(1) ss(4) = ss(4) / ss(1) - ss(2)*ss(2) ss(5) = ss(5) / ss(1) - ss(2)*ss(3) ss(6) = ss(6) / ss(1) - ss(3)*ss(3) write(6,28) np,ss(1) write(6,29) ss(2),ss(4),ss(3),ss(5),ss(6) end do ! loop on kint stop end program chex102t 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 SUBROUTINE TRIQDS(N,VX,VY,XA,YA,W) ! COMPUTES THE ABSCISSAS AND WEIGHTS FOR INTEGRATING OVER A TRIANGLE ! WITH VERTICES (VX(1),VY(1)), (VX(2),VY(2)), AND (VX(3),VY(3)) ! USING SIMPSON'S ANALOGUE WITH N SUBDIVISIONS ON A SIDE ! ! N INTEGER NUMBER OF SUBDIVISIONS ON A SIDE ! VX REAL VECTOR X-COORDINATES OF VERTICES OF TRIANGLE ! VY REAL VECTOR Y-COORDINATES OF VERTICES OF TRIANGLE ! XA REAL VECTOR X-COORDINATES OF ABSCISSAS ! YA REAL VECTOR Y-COORDINATES OF ABSCISSAS ! W REAL VECTOR WEIGHTS CORRESPONDING TO ABSCISSAS ! ! NOTE 1) IF MIDPOINT ANALOGUE IS DESIRED, THEN TAKE THE FIRST N*N ! ABSCISSAS, AND CHANGE THE WEIGHTS BY MULTIPLYING BY (4/3) ! 2) TOTAL NUMBER OF POINTS IS (3*N*(N+1)+2)/2 ! ! J F MONAHAN (OCT,1989) DEPT OF STATISTICS, NCSU, RALEIGH, NC 27695 ! Recoded November 1999, April 2000 for Fortran 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(3), INTENT(IN) :: VX,VY REAL, DIMENSION( (3*N*(N+1)+2)/2 ), INTENT(OUT) :: XA,YA,W REAL A,B,C,D,E,F,DET,FN,R13N,R23N,WC,X,Y INTEGER NP1,I,J,KP,KW ! ! FIRST GET THE TRANSFORMATION FOR CHANGING FROM THE ! FUNDAMENTAL TRIANGLE WITH VERTICES (0,0),(1,0), (1,1) ! TO THE NEW TRIANGLE E = VX(1) F = VY(1) A = VX(2) - VX(1) C = VY(2) - VY(1) B = VX(3) - VX(2) D = VY(3) - VY(2) ! TRANSFORMATION NEW X = A*X + B*Y + E ! NEW Y = C*X + D*Y + F ! DET = ABS( A*D - B*C ) FN = REAL(N) R13N = 1./(3.*FN) R23N = 2./(3.*FN) ! GET THE WEIGHT FOR THE CENTROIDS WC = (3.*DET)/(8*N*N) ! FIRST GET THE CENTROIDS KP = 0 DO I = 1,N DO J = 1,N KP = KP + 1 IF( I .GE. J ) THEN X = (I-1)/FN + R23N Y = (J-1)/FN + R13N ELSE X = (J-1)/FN + R13N Y = (I-1)/FN + R23N END IF ! ( I .GE. J ) ! TRANSFORM FROM FUNDAMENTAL TRIANGLE XA(KP) = A*X + B*Y + E YA(KP) = C*X + D*Y + F W(KP) = WC END DO ! LOOP ON J END DO ! LOOP ON I ! NEXT DO THE VERTICES NP1 = N + 1 WC = WC/9. DO I = 1,NP1 DO J = 1,I KP = KP + 1 KW = 6 ! BE CAREFUL TO COUNT TRIANGLES 6 OR 3 OR 1 IF( (I.EQ.1) .OR. (I.EQ.NP1) ) KW = KW / 2 IF( (J.EQ.1) .OR. (J.EQ.I) ) KW = KW / 2 ! YES, WE WANT INTEGER DIVISION HERE (3/2 = 1) X = (I-1)/FN Y = (J-1)/FN XA(KP) = A*X + B*Y + E YA(KP) = C*X + D*Y + F W(KP) = KW*WC END DO ! LOOP ON J END DO ! LOOP ON I RETURN END SUBROUTINE TRIQDS ! *** end of filename chex102t.f95 *********************