! Last change: DOS 2 Aug 2000 9:30 pm ! *** copyright 2000 *** ! *** filename rho1.f95 *** John F. Monahan ** ! ********************** program rho1 ! demonstration of kernel 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 real, dimension(1000) :: r real a,b,c,fnm1,fnm2,rho,onemrr,sr1mr2,rm,rv,y,h1,h2,h3 real ran,gchirv,gnroul,eknlde ! passing a function to eknlde real, external :: ds3uns integer i ! 21 format(f8.4,3f16.6) ! ! put output here open( unit=6, file='rho1.out' ) ! 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(i) = ( a*rho + sr1mr2*b ) / & & sqrt( rho*rho*a*a + 2*rho*sr1mr2*a*b + onemrr*(b*b+c*c) ) rm = rm + r(i) if( i .gt. 1 ) rv = rv + (rm - i*r(i))*(rm - i*r(i))/(i*(i-1)) end do ! loop on i ! have data, now do some density estimation ! ! first sort the data call hsort(r,nrep) ! construct three values of h rv = sqrt( rv / float(nrep-1) ) rm = real(nrep) ** .2 ! scale factor 3 since kernel has variance 1/9 h1 = 3. * .6 * rv / rm h2 = 3. * .9 * rv / rm h3 = 3. * 1.2 * rv / rm ! now evaluate density estimate at lots of points do i = 1,np y = real(2*i-1-np)/real(np) a = eknlde(ds3uns,y,h1,r,nrep) b = eknlde(ds3uns,y,h2,r,nrep) c = eknlde(ds3uns,y,h3,r,nrep) write(6,21) y,a,b,c end do ! loop on i stop end program rho1 REAL FUNCTION EKNLDE(K,X,H,XK,N) ! EVALUATES KERNEL DENSITY ESTIMATE AT X USING KERNEL K WITH PARAM H ! *** KERNEL MUST HAVE FINITE SUPPORT *** LIMITED TO (-1,+1) *** ! *** VALUES IN XK MUST BE SORTED IN NONDECREASING ORDER ! J F MONAHAN JAN 1982 DEPT OF STAT, N C S U, RALEIGH, NC 27650 ! RECODED FEBRUARY 2000 FOR FORTRAN 95 IMPLICIT NONE REAL, INTENT(IN) :: X,H INTEGER, INTENT(IN) :: N REAL, DIMENSION(N) :: XK REAL K INTEGER JXL,JXR,I INTEGER IFIND ! FIRST FIND RANGE OF RELEVANT POINTS EKNLDE = 0. JXL = IFIND(X-H,XK,N) JXR = IFIND(X+H,XK,N) IF(JXR .EQ. JXL) RETURN ! IF NONE, RETURN ZERO JXL = JXL + 1 DO I = JXL,JXR EKNLDE=EKNLDE+K((X-XK(I))/H) END DO ! LOOP ON I EKNLDE = EKNLDE/(N*H) RETURN END FUNCTION EKNLDE REAL FUNCTION DS3UNS(Z) ! DENSITY FUNCTION OF THE SUM OF 3 UNIFORMS ! EACH INDEPENDENT ON ( -1/3, +1/3 ), SO SUPPORT ONLY ( -1, +1 ) IMPLICIT NONE REAL, INTENT(IN) :: Z REAL AZ AZ = ABS(Z) DS3UNS = 0. IF( AZ .GE. 1. ) RETURN ! SUPPORT ONLY ( -1, +1 ) IF( 3.*AZ .LT. 1. ) THEN DS3UNS = 1.125*(1.-3.*AZ*AZ) ELSE DS3UNS = 1.6875*(AZ-1.)*(AZ-1.) END IF ! ( 3.*AZ .LT. 1. ) RETURN END FUNCTION DS3UNS INTEGER FUNCTION IFIND(X,XK,N) ! FIND I SUCH THAT XK(I) LE X LT XK(I+1) ! IFIND=0 IF X LT XK(1) ! IFIND=N IF X GT XK(N) ! J F MONAHAN JAN 1982 DEPT OF STAT, N C S U, RALEIGH, NC 27650 ! RECODED FEBRUARY, 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION(N), INTENT(IN) :: XK INTEGER IL,IU ! SPECIAL CASES IFIND = 0 IF( X .LT. XK(1) ) RETURN IFIND = N IF( X .GE. XK(N) ) RETURN IL = 1 IU = N DO WHILE( IU-IL .GT. 1 ) IFIND = ( IU + IL ) / 2 IF( X .EQ. XK(IFIND) ) RETURN IF( X .LT. XK(IFIND) ) THEN IU = IFIND ELSE IL = IFIND END IF ! ( X .LT. XK(IFIND) ) END DO ! WHILE ( IU-IL .GT. 1 ) IFIND = IL RETURN END FUNCTION IFIND 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 SUBROUTINE HSORT(K,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K REAL KK INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HSORT 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 rho1.f95 *********************