! G22YDF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module g22ydfe_mod
! G22YDF Example Program Module:
! Parameters and User-defined Routines
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: fit_lm, read_line
! .. Parameters ..
Integer, Parameter, Public :: nin = 5, nout = 6
Contains
Subroutine read_line(v1)
! Read in a line from NIN and remove any comments
! .. Scalar Arguments ..
Character (*), Intent (Out) :: v1
! .. Local Scalars ..
Integer :: pend
! .. Intrinsic Procedures ..
Intrinsic :: adjustl, index
! .. Executable Statements ..
Read (nin,'(A200)') v1
pend = index(v1,'::')
If (pend/=0) Then
v1 = v1(1:pend-1)
End If
v1 = adjustl(v1)
Return
End Subroutine read_line
Subroutine fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)
! Perform a multiple linear regression using G02DAF
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: g02daf, g22znf, nag_wp
! .. Scalar Arguments ..
Type (c_ptr), Intent (In) :: hform
Integer, Intent (In) :: ip, mx, nobs
Character (*), Intent (In) :: intcpt
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: x(:,:), y(:)
Integer, Intent (In) :: isx(:)
Character (*), Intent (In) :: plab(:)
! .. Local Scalars ..
Real (Kind=nag_wp) :: rss, rvalue, tol
Integer :: i, idf, ifail, irank, ivalue, ldq, &
ldx, lwt, optype
Logical :: svd
Character (200) :: cvalue
Character (1) :: weight
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), q(:,:), &
res(:), se(:), wk(:), wt(:)
! .. Intrinsic Procedures ..
Intrinsic :: repeat, size, trim
! .. Executable Statements ..
ldx = size(x,1)
! We are assuming un-weighted data
weight = 'U'
lwt = 0
ldq = nobs
Allocate (b(ip),cov((ip*ip+ip)/2),h(nobs),p(ip*(ip+ &
2)),q(ldq,ip+1),res(nobs),se(ip),wk(ip*ip+5*(ip-1)),wt(lwt))
! Use suggested value for tolerance
tol = 0.000001E0_nag_wp
! Fit a regression model
ifail = 0
Call g02daf(intcpt,weight,nobs,x,ldx,mx,isx,ip,y,wt,rss,idf,b,se,cov, &
res,h,q,ldq,svd,irank,p,tol,wk,ifail)
! Get the formula for the model being fit
ifail = 0
Call g22znf(hform,'Formula',ivalue,rvalue,cvalue,optype,ifail)
! Display the results
Write (nout,*) 'Model: ', trim(cvalue)
Write (nout,*) ' Parameter Standard'
Write (nout,*) 'Coefficients Estimate Error'
Write (nout,*) repeat('-',51)
Do i = 1, ip
Write (nout,99997) plab(i), b(i), se(i)
End Do
Write (nout,*) repeat('-',51)
Write (nout,99998) 'Residual sum of squares = ', rss
Write (nout,99999) 'Degrees of freedom = ', idf
Return
99999 Format (1X,A,I9)
99998 Format (1X,A,F9.4)
99997 Format (1X,A30,1X,F7.3,5X,F7.3)
End Subroutine fit_lm
End Module g22ydfe_mod
Program g22ydfe
! .. Use Statements ..
Use g22ydfe_mod, Only: fit_lm, nin, nout, read_line
Use, Intrinsic :: iso_c_binding, Only: c_null_ptr, &
c_ptr
Use nag_library, Only: g22yaf, g22ybf, g22ycf, g22ydf, g22zaf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Type (c_ptr) :: hddesc, hform, hxdesc
Integer :: i, ifail, ip, lddat, ldx, lisx, &
lplab, lvinfo, lvnames, mx, nobs, &
nvar, sddat, sdx
Character (200) :: formula
Character (1) :: intcpt, what
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: dat(:,:), x(:,:), y(:)
Real (Kind=nag_wp) :: tx(0,0)
Integer, Allocatable :: isx(:), levels(:), vinfo(:)
Integer :: tisx(0), tvinfo(3)
Character (50), Allocatable :: plab(:), vnames(:)
Character (1) :: tplab(0)
! .. Executable Statements ..
Write (nout,*) 'G22YDF Example Program Results'
Write (nout,*)
hform = c_null_ptr
hddesc = c_null_ptr
hxdesc = c_null_ptr
! Skip heading in data file
Read (nin,*)
! Read in size of the data matrix and number of variable labels supplied
Read (nin,*) nobs, nvar, lvnames
! Read in number of levels and names for the variables
Allocate (levels(nvar),vnames(lvnames))
Read (nin,*) levels(1:nvar)
If (lvnames>0) Then
Read (nin,*) vnames(1:lvnames)
End If
! Create a description of the data matrix
ifail = 0
Call g22ybf(hddesc,nobs,nvar,levels,lvnames,vnames,ifail)
! Read in the data matrix and response variable
lddat = nobs
sddat = nvar
Allocate (dat(lddat,sddat),y(nobs))
Read (nin,*)(dat(i,1:nvar),y(i),i=1,nobs)
! Read in the formula for the full model, remove comments and parse it
Call read_line(formula)
ifail = 0
Call g22yaf(hform,formula,ifail)
! Calculate the size of the design matrix
ldx = 0
sdx = 0
ifail = 1
Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
If (ifail/=81 .And. ifail/=82 .And. ifail/=91 .And. ifail/=92) Then
! Redisplay any error messages not related to X being too small
ifail = 0
Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
End If
! Generate the design matrix
ldx = nobs
sdx = mx
Allocate (x(ldx,sdx))
ifail = 0
Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,x,ldx,sdx,mx,ifail)
what = 'X'
! Get size of output arrays
lvinfo = 3
lisx = 0
lplab = 0
ifail = 1
Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo, &
tvinfo,ifail)
If (ifail/=102) Then
! Redisplay the error message
ifail = 0
Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo, &
tvinfo,ifail)
End If
! Allocate output arrays
lisx = tvinfo(1)
lplab = tvinfo(2)
! We don't need VINFO as we are using labels in PLAB
lvinfo = 0
Allocate (isx(lisx),plab(lplab),vinfo(lvinfo))
! Get the ISX flag and parameter labels
ifail = 0
Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo, &
vinfo,ifail)
! Fit the full model and print the results
Call fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)
! Read in the formula for the sub-model, remove comments and parse it
Call read_line(formula)
ifail = 0
Call g22yaf(hform,formula,ifail)
! Get the ISX flag and parameter labels, as the new model has to be
! a sub-model of the original one, the output arrays, ISX, PLAB and VINFO
! can be reused as they will be of sufficient size
ifail = 0
what = 'S'
Call g22ydf(what,hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo, &
vinfo,ifail)
Write (nout,*)
Write (nout,*)
! Fit the sub-model and print the results
Call fit_lm(hform,intcpt,nobs,mx,x,isx,ip,y,plab)
! Clean-up the G22 handles
ifail = 0
Call g22zaf(hform,ifail)
Call g22zaf(hddesc,ifail)
Call g22zaf(hxdesc,ifail)
Deallocate (dat,x,y)
Deallocate (isx,levels,vinfo)
Deallocate (plab,vnames)
End Program g22ydfe