! G22YBF Example Program Text
! Mark 26.1 Release. NAG Copyright 2017.
Module g22ybfe_mod
! G22YBF Example Program Module:
! Parameters and User-defined Routines
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: construct_labels, 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 ..
Continue
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 construct_labels(ip,plab,vnames,vinfo)
! Construct a set of parameter labels using the information
! stored in VINFO
! .. Scalar Arguments ..
Integer, Intent (In) :: ip
! .. Array Arguments ..
Integer, Intent (In) :: vinfo(:)
Character (*), Allocatable, Intent (Out) :: plab(:)
Character (*), Intent (In) :: vnames(:)
! .. Local Scalars ..
Integer :: b, i, j, k, li, vi
Character (200) :: line, term
! .. Intrinsic Procedures ..
Intrinsic :: adjustl, trim
! .. Executable Statements ..
Continue
Allocate (plab(ip))
k = 1
Do j = 1, ip
b = vinfo(k)
k = k + 1
If (b==0) Then
line = 'Intercept'
Else
line = ''
Do i = 1, b
vi = vinfo(k)
li = vinfo(k+1)
If (li/=0) Then
Write (term,99999) trim(adjustl(vnames(vi))), li
Else
Write (term,99998) trim(adjustl(vnames(vi)))
End If
If (i==1) Then
line = trim(line) // trim(term)
Else
line = trim(line) // ' . ' // trim(term)
End If
! We are ignoring the contrast identifier when constructing
! these labels
k = k + 3
End Do
End If
plab(j) = line
End Do
99999 Format (A,' (Level ',I0,')')
99998 Format (A)
Return
End Subroutine construct_labels
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 ..
Continue
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 g22ybfe_mod
Program g22ybfe
! .. Use Statements ..
Use g22ybfe_mod, Only: construct_labels, 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
! .. 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,*) 'G22YBF 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 model formula, remove comments and parse it
Call read_line(formula)
ifail = 0
Call g22yaf(hform,formula,ifail)
! Start of constructing the design matrix ...
! Calculate the size of design matrix
ldx = 0
sdx = 0
ifail = 1
Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
If (ifail/=91) Then
! redisplay the error message
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)
! ... End of constructing the design matrix
! Start of getting the ISX vector and information on parameter labels ...
! Get size of output arrays used by G22YDF
lvinfo = 3
lisx = 0
lplab = 0
ifail = 1
Call g22ydf(hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo,tvinfo, &
ifail)
If (ifail/=92) Then
! redisplay the error message
ifail = 0
Call g22ydf(hform,hxdesc,intcpt,ip,lisx,tisx,lplab,tplab,lvinfo, &
tvinfo,ifail)
End If
! Allocate output arrays (we already know that LISX = MX, but G22YDF
! returns it just in case)
lisx = tvinfo(1)
lvinfo = tvinfo(3)
! We don't need PLAB as we are constructing our own labels from VINFO
Allocate (isx(lisx),vinfo(lvinfo))
! Get the ISX flag and parameter labels
ifail = 0
Call g22ydf(hform,hxdesc,intcpt,ip,lisx,isx,lplab,tplab,lvinfo,vinfo, &
ifail)
! Construct some verbose labels for the parameters
Call construct_labels(ip,plab,vnames,vinfo)
! ... End of getting the ISX vector and information on parameter labels
! Fit a regression 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 g22ybfe