Program g05pwfe
! G05PWF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
! .. Use Statements ..
Use nag_library, Only: g02gbf, g02gpf, g05kff, g05pwf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: lseed = 1, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: a, dev, eps, s, tol
Integer :: fn, fp, genid, i, idf, ifail, ip, &
iprint, irank, ldv, ldx, lstate, &
lwk, m, maxit, n, nn, np, nsamp, nt, &
nv, obs_val, pred_val, samp, sordx, &
subid, tn, tp, uset, usey
Logical :: vfobs
Character (1) :: errfn, link, mean, offset, weight
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: b(:), cov(:), eta(:), pred(:), &
se(:), seeta(:), sepred(:), t(:), &
v(:,:), wk(:), x(:,:), y(:)
Real (Kind=nag_wp) :: off(1), wt(1)
Integer, Allocatable :: isx(:), state(:)
Integer :: seed(lseed)
! .. Intrinsic Procedures ..
Intrinsic :: count, int, real
! .. Executable Statements ..
Write (nout,*) 'G05PWF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
! Set variables required by the regression (G02GBF) ...
! Read in the type of link function, whether a mean is required
! and the problem size
Read (nin,*) link, mean, n, m
! Set storage order for G05PWF (pick the one required by G02GBF and
! G02GPF)
sordx = 1
ldx = n
Allocate (x(ldx,m),y(n),t(n),isx(m))
! This example is not using an offset or weights
offset = 'N'
weight = 'U'
! Read in data
Read (nin,*)(x(i,1:m),y(i),t(i),i=1,n)
! Read in variable inclusion flags
Read (nin,*) isx(1:m)
! Read in control parameters for the regression
Read (nin,*) iprint, eps, tol, maxit
! Calculate IP
ip = count(isx(1:m)>0)
If (mean=='M' .Or. mean=='m') Then
ip = ip + 1
End If
! ... End of setting variables required by the regression
! Set variables required by data sampling routine (G05PWF) ...
! Read in the base generator information and seed
Read (nin,*) genid, subid, seed(1:lseed)
! Will always have a Y and T variable
usey = 1
uset = 1
! Query the required size of the STATE array
lstate = 0
Allocate (state(lstate))
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
! Reallocate STATE
Deallocate (state)
Allocate (state(lstate))
! Initialize the generator to a repeatable sequence
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
! Read in the size of the training set required
Read (nin,*) nt
! Read in the number of sub-samples we will use */
Read (nin,*) nsamp
! ... End of setting variables required by data sampling routine
! Set variables required by prediction routine (G02GPF) ...
! Regression is performed using G02GBF so error structure is binomial
errfn = 'B'
! This example does not use the predicted standard errors, so
! it doesn't matter what VFOBS is set to
vfobs = .False.
! ... End of setting variables required by prediction routine
! Calculate the size of the validation dataset
nv = n - nt
! Allocate arrays
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))
Allocate (eta(nv),seeta(nv),pred(nv),sepred(nv))
! Initialize counts
tp = 0
tn = 0
fp = 0
fn = 0
! Loop over each sample
Do samp = 1, nsamp
! Split the data into training and validation datasets
ifail = 0
Call g05pwf(nt,n,m,sordx,x,ldx,usey,y,uset,t,state,ifail)
! Call routine to fit generalized linear model, with Binomial errors
! to training data
ifail = -1
Call g02gbf(link,mean,offset,weight,nt,x,ldx,m,isx,ip,y,t,wt,dev,idf, &
b,irank,se,cov,v,ldv,tol,maxit,iprint,eps,wk,ifail)
If (ifail/=0) Then
If (ifail<6) Then
Go To 100
End If
End If
! Predict the response for the observations in the validation dataset
ifail = 0
Call g02gpf(errfn,link,mean,offset,weight,nv,x(nt+1,1),ldx,m,isx,ip, &
t(nt+1),off,wt,s,a,b,cov,vfobs,eta,seeta,pred,sepred,ifail)
! Count the true/false positives/negatives
Do i = 1, nv
obs_val = int(y(nt+i))
If (pred(i)>=0.5_nag_wp) Then
pred_val = 1
Else
pred_val = 0
End If
Select Case (obs_val)
Case (0)
! Negative
Select Case (pred_val)
Case (0)
! True negative
tn = tn + 1
Case (1)
! False positive
fp = fp + 1
End Select
Case (1)
! Positive
Select Case (pred_val)
Case (0)
! False negative
fn = fn + 1
Case (1)
! True positive
tp = tp + 1
End Select
End Select
End Do
End Do
! Display results
np = tp + fn
nn = fp + tn
Write (*,99998) ' Observed'
Write (*,99998) ' --------------------------'
Write (*,99998) 'Predicted | Negative Positive Total'
Write (*,99998) '--------------------------------------'
Write (*,99997) 'Negative |', tn, fn, tn + fn
Write (*,99997) 'Positive |', fp, tp, fp + tp
Write (*,99997) 'Total |', nn, np, nn + np
Write (*,*)
If (np/=0) Then
Write (nout,99999) 'True Positive Rate (Sensitivity):', &
real(tp,kind=nag_wp)/real(np,kind=nag_wp)
Else
Write (nout,99998) &
'True Positive Rate (Sensitivity): No positives in data'
End If
If (nn/=0) Then
Write (nout,99999) 'True Negative Rate (Specificity):', &
real(tn,kind=nag_wp)/real(nn,kind=nag_wp)
Else
Write (nout,99998) &
'True Negative Rate (Specificity): No negatives in data'
End If
100 Continue
99999 Format (1X,A,F5.2)
99998 Format (1X,A)
99997 Format (1X,A,1X,I5,5X,I5,5X,I5)
End Program g05pwfe