Program g08chfe

!     G08CHF Example Program Text
!     Mark 27.0 Release. NAG Copyright 2019.

!     .. Use Statements ..
Use nag_library, Only: g05kff, g05sff, g08chf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: genid = 1, lseed = 1, mstate = 17,   &
nin = 5, nout = 6, subid = -1
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: a2, aa2, beta, nupper, p, sa2, sbeta
Integer                          :: i, ifail, j, k, lstate, n, nsim,     &
n_pseudo
Logical                          :: issort
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: x(:), xsim(:), y(:)
Integer                          :: seed(lseed), state(17)
!     .. Intrinsic Procedures ..
Intrinsic                        :: exp, real, sum
!     .. Executable Statements ..
Write (nout,*) 'G08CHF Example Program Results'
Write (nout,*)

!     Skip heading in data file

!     Memory allocation
Allocate (x(n),y(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
lstate = mstate
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