PROGRAM g02gpfe

!      G02GPF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g02gaf, g02gpf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: a, eps, rss, s, tol
       INTEGER                         :: i, idf, ifail, ip, iprint, irank,    &
                                          ldv, ldx, loff, lt, lwk, lwt, m,     &
                                          maxit, n, tldx
       LOGICAL                         :: vfobs
       CHARACTER (1)                   :: error, link, mean, offset, weight
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), eta(:), off(:),        &
                                          pred(:), se(:), seeta(:), sepred(:), &
                                          t(:), twt(:), tx(:,:), ty(:),        &
                                          v(:,:), wk(:), wt(:), x(:,:), y(:)
       INTEGER, ALLOCATABLE            :: isx(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          count
!      .. Executable Statements ..
       WRITE (nout,*) 'G02GPF Example Program Results'
       WRITE (nout,*)

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

!      Read in training data for model that will be used for prediction
       READ (nin,*) link, mean, offset, weight, n, m, s

       IF (weight=='W' .OR. weight=='w') THEN
          lwt = n
       ELSE
          lwt = 0
       END IF
       tldx = n
       ALLOCATE (tx(tldx,m),ty(n),twt(lwt),isx(m))

!      Read in data
       IF (lwt>0) THEN
          READ (nin,*) (tx(i,1:m),ty(i),twt(i),i=1,n)
       ELSE
          READ (nin,*) (tx(i,1:m),ty(i),i=1,n)
       END IF

!      Read in variable inclusion flags
       READ (nin,*) isx(1:m)

!      Calculate IP
       ip = count(isx(1:m)>0)
       IF (mean=='M' .OR. mean=='m') THEN
          ip = ip + 1
       END IF

!      Read in power for exponential link
       IF (link=='E' .OR. link=='e') THEN
          READ (nin,*) a
       END IF

       ldv = n
       lwk = (ip*ip+3*ip+22)/2
       ALLOCATE (b(ip),se(ip),cov(ip*(ip+1)/2),v(ldv,ip+7),wk(lwk))

!      Read in the offset
       IF (offset=='Y' .OR. offset=='y') THEN
          READ (nin,*) v(1:n,7)
       END IF

!      Read in control parameters
       READ (nin,*) iprint, eps, tol, maxit

!      Call routine to fit generalized linear model, with Normal errors, 
!      to training data
       ifail = -1
       CALL g02gaf(link,mean,offset,weight,n,tx,tldx,m,isx,ip,ty,twt,s,a,rss, &
          idf,b,irank,se,cov,v,ldv,tol,maxit,iprint,eps,wk,ifail)
       IF (ifail/=0) THEN
          IF (ifail<6) THEN
             GO TO 20
          END IF
       END IF

!      Display parameter estimates for training data
       WRITE (nout,99999) 'Residual sum of squares = ', rss
       WRITE (nout,99998) 'Degrees of freedom      = ', idf
       WRITE (nout,*)
       WRITE (nout,*) '      Estimate     Standard error'
       WRITE (nout,*)
       WRITE (nout,99997) (b(i),se(i),i=1,ip)

!      Skip second lot of headings in data file
       READ (nin,*)

!      Read in size of prediction data
       READ (nin,*) n, vfobs, offset, weight

!      Used G02GAF to fit training model, so error structure is normal is
!      do not require T array
       error = 'N'
       lt = 0

       ldx = n
       IF (weight=='W' .OR. weight=='w') THEN
          lwt = n
       ELSE
          lwt = 0
       END IF
       IF (offset=='Y' .OR. offset=='y') THEN
          loff = n
       ELSE
          loff = 0
       END IF
       ALLOCATE (x(ldx,m),y(n),wt(lwt),off(loff),eta(n),seeta(n),pred(n), &
          sepred(n),t(lt))

!      Read in predication data
       IF (lwt>0) THEN
          READ (nin,*) (x(i,1:m),wt(i),i=1,n)
       ELSE
          READ (nin,*) (x(i,1:m),i=1,n)
       END IF

!      Read in offset for the prediction data
       IF (offset=='Y' .OR. offset=='y') THEN
          READ (nin,*) off(1:n)
       END IF

!      Call prediction routine
       ifail = -1
       CALL g02gpf(error,link,mean,offset,weight,n,x,ldx,m,isx,ip,t,off,wt,s, &
          a,b,cov,vfobs,eta,seeta,pred,sepred,ifail)
       IF (ifail/=0) THEN
          WRITE (nout,99995) ifail
          IF (ifail/=22) THEN
             GO TO 20
          END IF
       END IF

!      Display predicted values
       WRITE (nout,*)
       WRITE (nout,*) '  I      ETA         SE(ETA)      Predicted  ', &
          ' SE(Predicted)'
       WRITE (nout,*)
       WRITE (nout,99996) (i,eta(i),seeta(i),pred(i),sepred(i),i=1,n)

20     CONTINUE

99999  FORMAT (1X,A,E12.4)
99998  FORMAT (1X,A,I0)
99997  FORMAT (1X,2F14.4)
99996  FORMAT (1X,I3,')',F10.5,3F14.5)
99995  FORMAT (1X,A,I5)
    END PROGRAM g02gpfe