!     Last change:  DOS   1 Aug 2000    8:09 pm
!                                       *** copyright 2000 ***
! *** filename ggchi.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PGGCHI
!  TEST OF RATIO OF UNIFORMS ALGORITHM FOR CHI DISDTRIBUTION
!                   let's check out 2**14 = 16384
      IMPLICIT NONE
      INTEGER, PARAMETER :: N = 16384
      INTEGER I,K,NDF
      REAL, DIMENSION(N) :: X
      REAL XS,X2S,FN,XI,XNMIP1,PHIXI,PHICXI,DKSP,DKSM,DKS,AD,ALPHA,DF
      REAL RAN,GCHIRV,PGAMA,QGAMA
!
  21  FORMAT(//' Sample Mean and Variance for',I12,' Chi Deviates'    &
     &  /9x,' with',I4,' degrees of freedom'                          &
     &  /'  Mean is',F16.8,4x,'Variance is',F16.8)
  22  FORMAT(' Goodness of Fit Statistics: '/4x,'Anderson-Darling is',&
     &F12.6/4x,'Compare to .05, .01 critical values of 2.492, 3.857'//&
     &  4x,'Kolmogorov-Smirnov   D+',F9.6,4x,'  D-',F9.6/             &
     &  4x,'SQRT(N)*MAX(D+,D-) is',F9.4/4x,'Compare to .05, .01',     &
     &  ' critical values of 1.36, 1.63')
!
!                   open output file
      OPEN( UNIT=6, FILE='ggchi.out' )
!                   let's check out 2**14 = 16384
      FN = REAL(N)
!                   initialize uniform generator
      XS = RAN(5151917)
!                   do a variety of shape parameters
      DO K = 1,5
      NDF = (1 + K*K )/2
      DF = REAL(NDF)
      ALPHA = DF/2.
!                   initialize moment sums
      X2S = 0.
      XS = 0.
      DO I = 1,N
      X(I) = GCHIRV(DF)
      XS = XS + X(I)
      IF( I .GT. 1 ) X2S = X2S + (XS - I*X(I))*(XS - I*X(I))/(I*(I-1))
      END DO  ! LOOP ON I
!                   get sample mean and variance
      XS = XS/FN
      X2S = X2S/(FN-1.)
      WRITE(6,21) N,NDF,XS,X2S
!                   sort the data -- get order statistics
      CALL HSORT(X,N)
!                   now compute K-S and A-D statistics
      DKSP = 0.
      DKSM = 0.
      AD = 0.
      DO I = 1,N
      XI = X(I)*X(I)/2.       ! transform to get gamma variates
      XNMIP1 = X(N+1-I)*X(N+1-I)/2.
      PHIXI = PGAMA(ALPHA,XI)
      PHICXI = QGAMA(ALPHA,XNMIP1)
      AD = AD - REAL(2*I-1)*( LOG(PHIXI) + LOG(PHICXI) )/FN - 1.
      DKSP = MAX( DKSP, REAL(I)/FN - PHIXI )
      DKSM = MAX( DKSM, PHIXI - REAL(I-1)/FN )
      END DO  ! LOOP ON I
!                   usual standardization of K-S
      DKS = SQRT(FN)*MAX( DKSP, DKSM )
      WRITE(6,22) AD,DKSP,DKSM,DKS
      END DO  ! LOOP ON K / ALPHA
      STOP
      END PROGRAM PGGCHI
      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
      real function pgama(a,x)
!  computes incomplete gamma function
!
!  integral from 0 to x of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then pgama = 0
!  if a <= 0 then pgama = -1
!  if convergence problems then pgama = -2
!
!  follows methods for regions I, III and tail as described in
!  Walter Gautschi (1979) "A Computational Procedure for Incomplete
!   Gamma Functions," ACM TOMS, Volume 5, pp. 466-481.
!
! J F Monahan (June 1997) Dept of Stat, NC State U, Raleigh, NC USA
!   recoded February, April 2000 for Fortran 95
      implicit none
      real, intent(in) :: a,x
      real, parameter :: eps = 1.e-6
      real t,r,st,ak,fk
      real algama
      integer k
!    first flag some out of range values
      pgama = -1.
      if( a .le. 0. ) return
      pgama = 0.
      if( x .le. 0. ) return
!
!    triage -- three methods depending on a and x
!
      if( a .gt. x +.25 ) then
!
!   monotone series for little g
!
      t = 1./a
      st = t
      do k = 1,50             ! iteration max is 50
      t = t*x/(a+real(k))
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 50 ) then
      pgama = st*exp(a*log(x)-x-algama(a))
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.
      t = 1.
      st = 1.
      do k = 1,50             ! iteration max is 50
      fk = real(k)
      ak = fk*(a-fk)/( (x+fk+fk-1.-a) * (x+fk+fk+1.-a) )
      r = -ak*(1.+r)/( 1.+ak*(1.+r) )
      t = t*r
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      pgama = 1. - st*exp(a*log(x)-x-algama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.
      r = 1.
      do k = 1,50             ! iteration max is 50
!                *** notice k+1 here
      fk = real(k+1)
      r = -x*r/fk
      t = (a+1.) * r / (a+fk)
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      pgama = (1.+1./a-st*x)*exp( a*log(x) - algama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
      end function pgama
      real function qgama(a,x)
!  computes incomplete gamma function
!
!  integral from x to infinity of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then qgama = 1
!  if a <= 0 then qgama = -1
!  if convergence problems then qgama = -2
!
!  follows methods for regions I, III and tail as described in
!  Walter Gautschi (1979) "A Computational Procedure for Incomplete
!   Gamma Functions," ACM TOMS, Volume 5, pp. 466-481.
!
! J F Monahan (June 1997) Dept of Stat, NC State U, Raleigh, NC USA
!   recoded February, April 2000 for Fortran 95
      implicit none
      real, intent(in) :: a,x
      real, parameter :: eps = 1.e-6
      real t,r,st,ak,fk
      real algama
      integer k
!    first flag some out of range values
      qgama = -1.
      if( a .le. 0. ) return
      qgama = 1.
      if( x .le. 0. ) return
!
!    triage -- three methods depending on a and x
!
      if( a .gt. x +.25 ) then
!
!   monotone series for little g
!
      t = 1./a
      st = t
      do k = 1,50             ! iteration max is 50
      t = t*x/(a+real(k))
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 50 ) then
      qgama = 1. - st*exp(a*log(x)-x-algama(a))
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.
      t = 1.
      st = 1.
      do k = 1,50             ! iteration max is 50
      fk = real(k)
      ak = fk*(a-fk)/( (x+fk+fk-1.-a) * (x+fk+fk+1.-a) )
      r = -ak*(1.+r)/( 1.+ak*(1.+r) )
      t = t*r
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      qgama = st*exp(a*log(x)-x-algama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.
      r = 1.
      do k = 1,50             ! iteration max is 50
!                *** notice k+1 here
      fk = real(k+1)
      r = -x*r/fk
      t = (a+1.) * r / (a+fk)
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      qgama = 1. - (1.+1./a-st*x)*exp( a*log(x) - algama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
      end function qgama
      REAL FUNCTION ALGAMA(X)
!  COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS
!   X .LE. 0.  RETURNS -1.0
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(5) :: P,Q         ! COEFS FOR RATIONAL FUNCTION
      REAL, DIMENSION(3) :: R           ! COEFS FOR STIRLING CORRECTION
      REAL XN,A,V
      REAL, PARAMETER :: LSQR2PI = .918938533
      DATA P/ -.510046056E1, -.368309363E2, -.664664448E2,  &
     &        -.195692802E3, -.185072903E1 /
      DATA Q/ -.114321842E2,  .368848969E2,  .16269403E2,   &
     &        -.195692802E3, 1.0 /
      DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 /      
!   USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS'
!   ED. BY JF HART,  SIAM SERIES,  WILEY(1968)
!   27 MAR 74    JFM W/RWS,RM
!                         ** MODIFICATIONS   JUNE 1988, 1993 **
!    RECODED OCTOBER 1999 FOR FORTRAN 95
      IF( X .LE. 0. ) THEN
      ALGAMA = -1.0            ! NOT DESIGNED FOR NEGATIVE ARGUMENTS
      RETURN
                        ENDIF
!                          FOR LARGE X (>8) USE STIRLING ROUTE
      IF( X .GT. 8. ) THEN
!                           LARGE X -- USE STIRLING WITH CORRECTION
      XN = X * X
!                           CORRECTION FOR STIRLING
      V = R(2) + ( R(3)/XN + R(1) )/XN                                  
      ALGAMA = (X-0.5)*LOG(X) - X + V/X + LSQR2PI
                      ELSE
!                          FOR SMALL X, USE REPRODUCTION FORMULA
!                           TO GET NUMBER BETWEEN 2 AND 3
      V = 1.0                                                             
      XN = X                                                              
      DO WHILE( XN .LT. 2. )
      V = V/XN
      XN = XN + 1.
      END DO  ! WHILE( XN .LT. 2. )
      DO WHILE( XN .GT. 3. )
      XN = XN - 1.
      V = V * XN
      END DO  ! WHILE( XN .GT. 3. )
      A = XN-2.0
!                           RATIONAL FUNCTION APPROXIMATION
      V = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) /  &
     &   (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5)))))
      ALGAMA = LOG(V)
                      ENDIF
      RETURN
      END FUNCTION ALGAMA
! *** end of filename ggchi.f95         *********************