! Last change: DOS 2 Aug 2000 9:34 pm ! *** copyright 2000 *** ! *** filename rho3.f95 *** John F. Monahan ** ! ********************** program rho3 ! demonstration of spline density estimation ! ! get the distribution of the sample correlation coefficient ! for a value of rho not equal to zero ! ! implicit none integer, parameter :: n = 12 ! sample size for r integer, parameter :: nrep = 1000 ! replication size integer, parameter :: np = 100 ! density evaluations integer, parameter :: mknots = 33 ! max knots for spline real :: r ! no need to store real, dimension(mknots) :: bins,yk,h,m real, dimension(np,3) :: fhat real a,b,c,fnm1,fnm2,rho,onemrr,sr1mr2,rm,rv,y real ran,gchirv,gnroul,spled integer i,j,nkots,nbins ! 21 format(f8.4,3f16.6) ! ! put output here open( unit=6, file='rho3.out' ) nkots = 33 ! start with full number of knots nbins = nkots - 1 ! initialize bin count and fill yk yk(1) = -1. bins(nkots) = 0. do j = 1,nbins bins(j) = 0. yk(j+1) = real(2*j-nbins) / real(nbins) end do ! loop on j ! number of observations in sample fnm1 = real(n-1) fnm2 = real(n-2) ! initialize generator rm = ran(5151917) rm = 0. rv = 0. ! true value of correlation rho = .5 onemrr = 1. - rho*rho sr1mr2 = sqrt( onemrr ) ! generate observations do i = 1,nrep ! get sample correlation coefficient in three pieces ! a is chi with n-1 df ! b is normal( 0, 1) ! c is chi with n-2 df a = gchirv(fnm1) b = gnroul(i) c = gchirv(fnm2) r = ( a*rho + sr1mr2*b ) / & & sqrt( rho*rho*a*a + 2*rho*sr1mr2*a*b + onemrr*(b*b+c*c) ) rm = rm + r if( i .gt. 1 ) rv = rv + (rm - i*r)*(rm - i*r)/(i*(i-1)) ! place observation in correct bin j = int( ((1.+r)/2.)*nbins ) + 2 bins(j) = bins(j) + 1. end do ! loop on i ! have data, now do some density estimation ! ! first the set-up for spline ! do j = 2,nkots bins(j) = bins(j) + bins(j-1) end do ! loop on j do j = 1,nkots bins(j) = bins(j)/real(nrep) end do ! loop on j ! natural spline since expect zero derivs at endpoints call splstn(yk,bins,m,h,nkots) ! now evaluate density estimate at lots of points do i = 1,np y = real(2*i-1-np)/real(np) fhat(i,1) = spled(y,yk,bins,m,h,nkots) end do ! loop on i ! do it again with fewer bins nkots = 17 ! fewer knots do j = 2,nkots bins(j) = bins(2*j-1) yk(j) = yk(2*j-1) end do ! loop on j ! natural spline since expect zero derivs at endpoints call splstn(yk,bins,m,h,nkots) ! now evaluate density estimate at lots of points do i = 1,np y = real(2*i-1-np)/real(np) fhat(i,2) = spled(y,yk,bins,m,h,nkots) end do ! loop on i ! do it again with even fewer bins nkots = 9 do j = 2,nkots bins(j) = bins(2*j-1) yk(j) = yk(2*j-1) end do ! loop on j ! natural spline since expect zero derivs at endpoints call splstn(yk,bins,m,h,nkots) ! now evaluate density estimate at lots of points do i = 1,np y = real(2*i-1-np)/real(np) fhat(i,3) = spled(y,yk,bins,m,h,nkots) end do ! loop on i ! write out results do i = 1,np y = real(2*i-1-np)/real(np) write(6,21) y,fhat(i,1),fhat(i,2),fhat(i,3) end do ! loop on i stop end program rho3 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 SPLSTN(XK,Y,M,H,N) ! NATURAL CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( N ), INTENT(IN) :: XK,Y REAL, DIMENSION( N ), INTENT(OUT) :: M,H REAL, DIMENSION(3*N) :: W REAL S INTEGER J,NM1 ! NM1 = N-1 ! COMPUTE DIFFERENCES H(J) AND FIRST DIFS OF Y'S DO J = 2,N H(J) = XK(J)-XK(J-1) M(J) = (Y(J)-Y(J-1))/H(J) ! S(J) END DO ! LOOP ON J ! FILL MATRIX TRIDIAGONAL MATRIX A IN VECTOR W ! DIAGONAL W(J) = A(J,J) = 2 ! SUBDIAGONAL W(N+J) = A(J,J-1) = LAMBDA(J) ! SUPERDIAGONAL W(2*N+J) = A(J,J+1) = MU(J) ! ! M(J) NOW HOLDS FIRST DIFFERENCES ! ! FIX TWO ENDS DO J = 2,NM1 S = H(J)+H(J+1) W(2*N+J) = H(J+1)/S ! SUPERDIAGONAL W(N+J) = H(J)/S ! SUBDIAGONAL W(J) = 2. ! DIAGONAL M(J) = 6.*(M(J+1)-M(J))/S ! SECOND DIFFERENCES END DO ! LOOP ON J ! FINISH ENDS W(1) = 2. W(N) = 2. W(2*N+1) = 0. W(2*N) = 0. M(1) = 0. M(N) = 0. ! SOLVE TRIDIAGONAL SYSTEM CALL STRIDS(W,M,N) RETURN END SUBROUTINE SPLSTN SUBROUTINE STRIDS(W,X,N) ! SOLVES TRIDIAGONAL SYSTEM: A*(X-OUT)=(X-IN) ! WHERE DIAGONAL ELEMENTS IN W(I) = A(I,I) ! SUBDIAGONAL ELEMENTS IN W(N+I) = A(I,I-1) ! SUPERDIAGONAL ELEMENTS IN W(2*N+I) = A(I,I+1) ! ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( 3*N ), INTENT(IN OUT) :: W REAL, DIMENSION( N ), INTENT(IN OUT) :: X INTEGER I,K ! COMPUTE LU FACTORIZATION (L IS UNIT LOWER) DO K = 2,N W(N+K) = W(N+K)/W(K-1) W(K) = W(K) - W(N+K)*W(2*N+K-1) END DO ! LOOP ON K ! SOLVE LOWER TRIANGULAR (BIDIAGONAL) SYSTEM DO I = 2,N X(I) = X(I) - W(N+I)*X(I-1) END DO ! LOOP ON I ! SOLVE UPPER TRIANGULAR SYSTEM (START AT BOTTOM) X(N) = X(N) / W(N) DO K = 2,N I = N+1-K X(I) = ( X(I) - W(2*N+I)*X(I+1) ) / W(I) END DO ! LOOP ON K RETURN END SUBROUTINE STRIDS REAL FUNCTION SPLEV(X,XK,Y,M,H,N) ! EVALUATES SPLINE AT X XK(JC-1) LT X LE XK(JC) ! CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,M,H REAL DA,DB INTEGER, SAVE :: J = 0 INTEGER JFIND ! INITIAL VALUE -- RETURN THIS IF OUT OF RANGE SPLEV = 0. ! IF X OUT OF RANGE, SPLEV=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! IF LAST CALL IN SAME INTERVAL, NO NEED TO SEARCH IF( ( J .LE. 1 ) .OR. ( J .GT. N ) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ELSE IF( (X .LE. XK(J-1)) .OR. (X .GT. XK(J)) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ENDIF ENDIF ! COMPUTE DIFFERENCES FROM NEARBY KNOTS DB = X-XK(J-1) DA = XK(J)-X ! EVALUATE CUBIC SPLINE SPLEV = ( (M(J-1)*DA**3 + M(J)*DB**3)/6. & & + (Y(J-1)-M(J-1)*H(J)*H(J)/6.)*DA & & + (Y(J)-M(J)*H(J)*H(J)/6.)*DB ) / H(J) RETURN END FUNCTION SPLEV REAL FUNCTION SPLED(X,XK,Y,M,H,N) ! EVALUATES DERIVATIVE OF SPLINE AT X XK(JC-1) LT X LE XK(JC) ! CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED FEBRUARY, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,M,H REAL DA,DB INTEGER, SAVE :: J = 0 INTEGER JFIND ! INITIAL VALUE -- RETURN THIS IF OUT OF RANGE SPLED = 0. ! IF X OUT OF RANGE, SPLEV=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! IF LAST CALL IN SAME INTERVAL, NO NEED TO SEARCH IF( ( J .LE. 1 ) .OR. ( J .GT. N ) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ELSE IF( (X .LE. XK(J-1)) .OR. (X .GT. XK(J)) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ENDIF ENDIF ! COMPUTE DIFFERENCES FROM NEARBY KNOTS DB = X-XK(J-1) DA = XK(J)-X ! EVALUATE DERIVATIVE OF CUBIC SPLINE SPLED = 0.5*( M(J)*DB*DB - M(J-1)*DA*DA )/H(J) & & + (Y(J)-Y(J-1))/H(J) - H(J)*(M(J)-M(J-1))/6. RETURN END FUNCTION SPLED INTEGER FUNCTION JFIND(X,XK,N) ! FIND J SUCH THAT XK(J-1) LT X LE XK(J) ! JFIND=2 IF X EQ XK(1) JFIND=0 IF X IS OUT OF RANGE ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK INTEGER JU,JL,I JFIND = 0 ! IF X OUT OF RANGE, JFIND=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! CHECK SPECIAL CASE IF( X .EQ. XK(1) ) THEN JFIND = 2 RETURN END IF ! INITIALIZE POINTERS JU = N JL = 1 DO WHILE( JU-JL .GT. 1 ) ! TEST THIS ONE I = ( JU + JL ) / 2 IF( X .LE. XK(I) ) THEN JU = I ELSE JL = I END IF END DO ! WHILE( JU-JL .GT. 1 ) JFIND = JU RETURN END FUNCTION JFIND ! *** end of filename rho3.f95 *********************