PROGRAM g08chfe

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g05kff, g05sff, g08chf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: genid = 1, lseed = 1, nin = 5,       &
                                          nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: a2, aa2, beta, nupper, p, sa2, sbeta
       INTEGER                         :: i, ifail, j, k, lstate = 17, n,      &
                                          nsim, n_pseudo, subid = -1
       LOGICAL                         :: issort
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: x(:), xsim(:), y(:)
       INTEGER                         :: seed(lseed), state(17)
!      .. Intrinsic Functions ..
       INTRINSIC                          exp, real, sum
!      .. Executable Statements ..
       WRITE (nout,*) 'G08CHF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read number of observations
       READ (nin,*) n

!      Memory allocation
       ALLOCATE (x(n),y(n))

!      Read observations
       READ (nin,*) (x(i),i=1,n)

!      Maximum likelihood estimate of mean
       beta = sum(x(1:n))/real(n,kind=nag_wp)

!      PIT, using exponential CDF with mean beta
       DO i = 1, n
          y(i) = 1.0E0_nag_wp - exp(-x(i)/beta)
       END DO

!      Let g08chf sort the (approximately) uniform variates
       issort = .FALSE.

!      Calculate A-squared
       ifail = 0
       a2 = g08chf(n,issort,y,ifail)
       aa2 = (1.0E0_nag_wp+0.6E0_nag_wp/real(n,kind=nag_wp))*a2

!      Number of simulations
       nsim = 888

!      Generate exponential variates using a repeatable seed
       ALLOCATE (xsim(n*nsim))
       seed(1) = 206033
       ifail = 0
       CALL g05kff(genid,subid,seed,lseed,state,lstate,ifail)
       n_pseudo = n*nsim
       ifail = 0
       CALL g05sff(n_pseudo,beta,state,xsim,ifail)

!      Simulations loop
       nupper = 0.0E0_nag_wp
       DO j = 1, nsim
          k = (j-1)*n
!         Maximum likelihood estimate of mean
          sbeta = sum(xsim(k+1:k+n))/real(n,kind=nag_wp)
!         PIT
          DO i = 1, n
             y(i) = 1.0E0_nag_wp - exp(-xsim(k+i)/sbeta)
          END DO
!         Calculate A-squared
          ifail = 0
          sa2 = g08chf(n,issort,y,ifail)
          IF (sa2>aa2) THEN
             nupper = nupper + 1.0E0_nag_wp
          END IF
       END DO

!      Simulated upper tail probability value
       p = nupper/real(nsim+1,kind=nag_wp)

!      Results
       WRITE (nout,'(1X,A,E11.4)') &
          'H0: data from exponential distribution with mean', beta
       WRITE (nout,'(1X,A,1X,F8.4)') 'Test statistic, A-squared: ', a2
       WRITE (nout,'(1X,A,1X,F8.4)') 'Upper tail probability:    ', p

    END PROGRAM g08chfe