! Last change: DOS 2 Aug 2000 9:32 pm ! *** copyright 2000 *** ! *** filename rho2.f95 *** John F. Monahan ** ! ********************** program rho2 ! demonstration of orthogonal series 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 ! evaluations of density est ! number of terms in orthogonal series integer, parameter :: max = 15 real, dimension(0:max) :: v,cc real, dimension(1000) :: r real a,b,c,fnm1,fnm2,rho,onemrr,sr1mr2,y,rm,rv real ran,gchirv,gnroul,ortpe integer i,j ! 21 format(f8.4,3f16.6) ! ! put output here open( unit=6, file='rho2.out' ) ! number of observations in sample fnm1 = real(n-1) fnm2 = real(n-2) ! initialize generator r(1) = 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, get some simple statistics rm = rm / real(nrep) rv = sqrt( rv / real(nrep-1) ) ! ! *** no need to sort the data *** ! call hsort(r,nrep) ! have data, now do some density estimation ! ! first the set-up for orthogonal series ! ! compute Fourier coefficients cc = 0. ! initialize to zero do i = 1,nrep y = r(i) call ortpv(y,1,max,v) do j = 0,max cc(j) = cc(j) + v(j) end do ! loop on j end do ! loop on i ! put on normalization constant (actually twice) call ortpn(1,max,v) ! coefficients have an extra sqrt(norm) do j = 0,max cc(j) = cc(j) / ( float(nrep)*v(j) ) end do ! loop on jntinue ! now evaluate density estimate at lots of points do i = 1,np y = float(2*i-1-np)/float(np) a = ortpe(y,1,4,cc) b = ortpe(y,1,8,cc) c = ortpe(y,1,12,cc) write(6,21) y,a,b,c end do ! loop on i stop end program rho2 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 REAL FUNCTION ORTPE(X,KIND,NMAX,C) ! USES CLENSHAW'S METHOD TO EVALUATE A FUNCTION GIVEN AS A FINITE ! EXPANSION IN ORTHOGONAL POLYNOMIALS WITH COEFFICIENTS STORED ! IN C(0),...,C(NMAX) ! ! KIND=1 LEGENDRE W(X)=1 ! KIND=2 CHEBYSHEV, FIRST W(X)=1/SQRT(1-X**2) ! KIND=3 CHEBYSHEV, SECOND W(X)=SQRT(1-X**2) ! KIND=4 STANDARD FOURIER SERIES W(X)=1. ! KIND=5 LAGUERRE (STANDARD, ALPHA=0) W(X)=EXP(-X) ! KIND=6 HERMITE W(X)=EXP(-X**2) ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, INTENT(IN) :: X INTEGER, INTENT(IN) :: KIND,NMAX REAL, DIMENSION(0:NMAX), INTENT(IN) :: C REAL, PARAMETER :: PI = 3.141593 ! 3.14159265358979 REAL YNEW,YOLD,YOLDR INTEGER L,LL,K ! YOLD = 0. YOLDR = 0. ORTPE = 0. SELECT CASE (KIND) ! LEGENDRE CASE (1) DO L = 0,NMAX K = NMAX - L YNEW = REAL(2*K+1)*X*YOLD/REAL(K+1) - & & REAL(K+1)*YOLDR/REAL(K+2) + C(K) YOLDR = YOLD YOLD = YNEW END DO ! LOOP ON L ORTPE = YNEW ! CHEBYSHEV, FIRST CASE (2) DO L = 1,NMAX K = NMAX + 1 - L YNEW = 2.*X*YOLD - YOLDR + C(K) YOLDR = YOLD YOLD = YNEW END DO ! LOOP ON L ORTPE = X*YOLD - YOLDR + C(0) ! CHEBYSHEV, SECOND CASE (3) DO L = 0,NMAX K = NMAX - L YNEW = 2.*X*YOLD - YOLDR + C(K) YOLDR = YOLD YOLD = YNEW END DO ! LOOP ON L ORTPE = YNEW ! STANDARD FOURIER SERIES CASE (4) ORTPE = 0. DO LL = 1,NMAX L = NMAX + 1 - LL K = (L+1)/2 IF( L .EQ. 2*K ) THEN ! EVEN ORTPE = ORTPE + C(L)*SIN(K*PI*X) ELSE ORTPE = ORTPE + C(L)*COS(K*PI*X) END IF ! ( L .EQ. 2*K ) END DO ! LOOP ON LL ORTPE = ORTPE + C(0) ! LAGUERRE CASE (5) DO L = 0,NMAX K = NMAX - L YNEW = (REAL(2*K+1)-X)*YOLD/REAL(K+1) - & & REAL(K+1)*YOLDR/REAL(K+2) + C(K) YOLDR = YOLD YOLD = YNEW END DO ! LOOP ON L ORTPE = YNEW ! HERMITE CASE (6) DO L = 0,NMAX K = NMAX - L YNEW = 2.*X*YOLD - 2.*REAL(K+1)*YOLDR + C(K) YOLDR = YOLD YOLD = YNEW END DO ! LOOP ON L ORTPE = YNEW CASE DEFAULT END SELECT RETURN END FUNCTION ORTPE SUBROUTINE ORTPV(X,KIND,NMAX,V) ! COMPUTES VALUES OF ORTHOGONAL POLYS FROM 0 THRU NMAX ! IN V(0),,.V(NMAX) ! ! KIND=1 LEGENDRE W(X)=1 ! KIND=2 CHEBYSHEV, FIRST W(X)=1/SQRT(1-X**2) ! KIND=3 CHEBYSHEV, SECOND W(X)=SQRT(1-X**2) ! KIND=4 STANDARD FOURIER SERIES W(X)=1. ! KIND=5 LAGUERRE (STANDARD, ALPHA=0) W(X)=EXP(-X) ! KIND=6 HERMITE W(X)=EXP(-X**2) ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, INTENT(IN) :: X INTEGER, INTENT(IN) :: KIND, NMAX REAL, DIMENSION(0:NMAX), INTENT(OUT) :: V REAL, PARAMETER :: PI = 3.141593 ! 3.14159265358979 INTEGER N,K SELECT CASE (KIND) ! LEGENDRE CASE(1) V(0) = 1. V(1) = X DO N = 2,NMAX V(N) = ( REAL(2*N-1)*X*V(N-1) - REAL(N-1)*V(N-2) )/REAL(N) END DO ! LOOP ON N ! CHEBYSHEV, FIRST CASE (2) V(0) = 1. V(1) = X DO N = 2,NMAX V(N) = 2.*X*V(N-1) - V(N-2) END DO ! LOOP ON N ! CHEBYSHEV, SECOND CASE (3) V(0) = 1. V(1) = 2.*X DO N = 2,NMAX V(N) = 2.*X*V(N-1) - V(N-2) END DO ! LOOP ON N ! STANDARD FOURIER SERIES CASE (4) V(0) = 1. DO N = 1,NMAX K = (N+1)/2 IF( N .EQ. 2*K ) THEN ! EVEN V(N) = SIN(K*PI*X) ELSE V(N) = COS(K*PI*X) END IF ! ( N .EQ. 2*K ) END DO ! LOOP ON N ! LAGUERRE CASE (5) V(0) = 1. V(1) = 1.-X DO N = 2,NMAX V(N) = ( (REAL(2*N-1)-X)*V(N-1) - REAL(N-1)*V(N-2) ) / REAL(N) END DO ! LOOP ON N ! HERMITE CASE (6) V(0) = 1. V(1) = 2.*X DO N = 2,NMAX V(N) = 2.*( X*V(N-1) - REAL(N-1)*V(N-2) ) END DO ! LOOP ON N CASE DEFAULT END SELECT RETURN END SUBROUTINE ORTPV REAL FUNCTION WEIGHT(X,KIND) ! COMPUTES WEIGHT FUNCTION FOR ORTHOGONAL POLYS ! KIND=1 LEGENDRE W(X)=1 ! KIND=2 CHEBYSHEV, FIRST W(X)=1/SQRT(1-X**2) ! KIND=3 CHEBYSHEV, SECOND W(X)=SQRT(1-X**2) ! KIND=4 STANDARD FOURIER SERIES W(X)=1. ! KIND=5 LAGUERRE (STANDARD, ALPHA=0) W(X)=EXP(-X) ! KIND=6 HERMITE W(X)=EXP(-X**2) ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE REAL, INTENT(IN) :: X INTEGER, INTENT(IN) :: KIND SELECT CASE (KIND) ! LEGENDRE CASE (1) WEIGHT = 1. ! CHEBYSHEV, FIRST CASE (2) WEIGHT = 1./ SQRT(1.-X*X) ! CHEBYSHEV, SECOND CASE (3) WEIGHT = SQRT(1.-X*X) ! STANDARD FOURIER SERIES CASE (4) WEIGHT = 1. ! LAGUERRE CASE (5) IF( X .GT. 85. ) THEN WEIGHT = 0. ! UNDERFLOW ELSE WEIGHT = EXP(-X) ! OK END IF ! ( X .GT. 85. ) ! HERMITE CASE (6) IF( X .GT. 9. ) THEN WEIGHT = 0. ! UNDERFLOW ELSE WEIGHT = EXP(-X*X) ! OK END IF ! ( X .GT. 9. ) CASE DEFAULT END SELECT RETURN END FUNCTION WEIGHT SUBROUTINE ORTPN(KIND,NMAX,V) ! COMPUTES NORMALIZATION CONSTANTS FOR ORTHOGONAL POLYS FROM ! 0 THRU NMAX AND STORES THEM IN V(0),...,V(NMAX) ! KIND=1 LEGENDRE W(X)=1 ! KIND=2 CHEBYSHEV, FIRST W(X)=1/SQRT(1-X**2) ! KIND=3 CHEBYSHEV, SECOND W(X)=SQRT(1-X**2) ! KIND=4 STANDARD FOURIER SERIES W(X)=1. ! KIND=5 LAGUERRE (STANDARD, ALPHA=0) W(X)=EXP(-X) ! KIND=6 HERMITE W(X)=EXP(-X**2) ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: KIND INTEGER, INTENT(IN) :: NMAX REAL, DIMENSION(0:NMAX), INTENT(OUT) :: V REAL, PARAMETER :: PI = 3.141593 ! 3.14159265358979 REAL, PARAMETER :: PIO2 = 1.570796 ! 1.5707963267949 REAL, PARAMETER :: SQRTPI = 1.772454 ! 1.77245385090552 INTEGER N SELECT CASE (KIND) ! LEGENDRE CASE (1) V(0) = 2. DO N = 1,NMAX V(N) = 2./REAL(2*N+1) END DO ! LOOP ON N ! CHEBYSHEV, FIRST CASE (2) V(0) = PI DO N = 1,NMAX V(N) = PIO2 END DO ! LOOP ON N ! CHEBYSHEV, SECOND CASE (3) V(0) = PIO2 DO N = 1,NMAX V(N) = PIO2 END DO ! LOOP ON N ! FOURIER CASE (4) V(0) = 2. DO N = 1,NMAX V(N) = 1. END DO ! LOOP ON N ! LAGUERRE CASE (5) V(0) = 1. DO N = 1,NMAX V(N) = 1. END DO ! LOOP ON N ! HERMITE CASE (6) V(0) = SQRTPI DO N = 1,NMAX V(N) = 2.*N*V(N-1) END DO ! LOOP ON N CASE DEFAULT END SELECT RETURN END SUBROUTINE ORTPN ! *** end of filename rho2.f95 *********************