! H05ABF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
Module h05abfe_mod
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: f, prepare_user_arrays
! .. Parameters ..
Integer, Parameter, Public :: nin = 5, nout = 6
Contains
Subroutine f(m,drop,lz,z,la,a,score,iuser,ruser,info)
! Score calculating function required by H05ABF
! M,DROP,LZ,Z,LA,A,SCORE IUSER and RUSER are all as described in the
! documentation of H05ABF.
! This particular example finds the set, of a given size, of explanatory
! variables that best fit a response variable when a linear regression
! model is used. Therefore the feature set is the set of all the
! explanatory variables and the best set of features is defined as set of
! explanatory variables that gives the smallest residual sums of squares.
! See the documentation for G02DAF for details on linear regression
! models.
! .. Use Statements ..
Use nag_library, Only: g02daf
! .. Scalar Arguments ..
Integer, Intent (In) :: drop, la, lz, m
Integer, Intent (Inout) :: info
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (Out) :: score(max(la,1))
Integer, Intent (In) :: a(la), z(lz)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: rss, tol
Integer :: ex, ey, i, idf, ifail, &
inv_drop, ip, irank, ldq, ldx, &
n, sx, sy
Logical :: svd
Character (1) :: mean, weight
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), &
q(:,:), res(:), se(:), wk(:), &
wt(:)
Integer, Allocatable :: isx(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, count, max
! .. Executable Statements ..
Continue
info = 0
n = iuser(1)
ldq = n
ldx = n
! No intercept term and no weights
mean = 'Z'
weight = 'U'
! Allocate various arrays required by G02DAF
Allocate (b(m),cov((m*m+m)/2),h(n),p(m*(m+2)),q(ldq,m+1),res(n),se(m), &
wk(m*m+5*(m-1)),wt(1),isx(m))
! Set up the initial feature set. If DROP = 0, this is the Null set (i.e.
! no features). If DROP = 1 then this is the full set (i.e. all features)
isx(1:m) = drop
! Add (if DROP = 0) or remove (if DROP = 1) the all the features specified
! in Z
inv_drop = abs(drop-1)
Do i = 1, lz
isx(z(i)) = inv_drop
End Do
! Get the start and end of X and Y in RUSER
sx = iuser(4)
ex = sx + m*n - 1
sy = iuser(3)
ey = sy + n - 1
! Extract some parameters from RUSER
tol = ruser(iuser(2))
Do i = 1, max(la,1)
If (la>0) Then
If (i>1) Then
! Reset the feature altered at the last iteration
isx(a(i-1)) = drop
End If
! Add or drop the I'th feature in A
isx(a(i)) = inv_drop
End If
ip = count(isx(1:m)==1)
! Fit the regression model
ifail = 0
Call g02daf(mean,weight,n,ruser(sx:ex),ldx,m,isx,ip,ruser(sy:ey),wt, &
rss,idf,b,se,cov,res,h,q,ldq,svd,irank,p,tol,wk,ifail)
! Return the score (the residual sums of squares)
score(i) = rss
End Do
! Keep track of the number of subsets evaluated
iuser(5) = iuser(5) + max(1,la)
End Subroutine f
Subroutine prepare_user_arrays(m,iuser,ruser)
! Populate the user arrays
! M is the maximum number of features (as per H05ABF)
! In this example RUSER holds the data required for the linear regression
! (the matrix of explanatory variables, X, the vector of response values,
! Y, and the tolerance, TOL) and IUSER holds number of observations and
! the location of the various elements in RUSER
! .. Scalar Arguments ..
Integer, Intent (In) :: m
! .. Array Arguments ..
Real (Kind=nag_wp), Allocatable, Intent (Out) :: ruser(:)
Integer, Allocatable, Intent (Out) :: iuser(:)
! .. Local Scalars ..
Integer :: i, ierr, j, liuser, lruser, n, &
p1, p2
Character (200) :: line
! .. Executable Statements ..
Continue
! Read in the number of observations for the data used in the linear
! regression skipping any headings or blank lines
Do
Read (nin,*,Iostat=ierr) line
If (ierr/=0) Exit
Read (line,*,Iostat=ierr) n
If (ierr==0) Exit
End Do
! Allocate space for the user arrays
liuser = 5
lruser = n + n*m + 1
Allocate (ruser(lruser),iuser(liuser))
! Number of observations
iuser(1) = n
! Location of TOL in RUSER
iuser(2) = 1
! Start of Y in RUSER
iuser(3) = 2
! Start of X in RUSER
iuser(4) = iuser(3) + n
! Keep track of the number of subsets evaluated
iuser(5) = 0
! Read in the tolerance for the regression
Read (nin,*) ruser(iuser(2))
! Read in the data
p1 = iuser(3)
p2 = iuser(4)
Do i = 0, n - 1
Read (nin,*)(ruser(p2+i+j*n),j=0,m-1), ruser(p1+i)
End Do
End Subroutine prepare_user_arrays
End Module h05abfe_mod
Program h05abfe
! .. Use Statements ..
Use nag_library, Only: h05abf, nag_wp
Use h05abfe_mod, Only: f, nin, nout, prepare_user_arrays
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: gamma
Integer :: i, ifail, ip, j, la, m, mincnt, &
mincr, mip, nbest
! .. Local Arrays ..
Real (Kind=nag_wp) :: acc(2)
Real (Kind=nag_wp), Allocatable :: bscore(:), ruser(:)
Integer, Allocatable :: a(:), bz(:,:), ibz(:), iuser(:), &
z(:)
Logical, Allocatable :: mask(:)
! .. Intrinsic Procedures ..
Intrinsic :: max, pack
! .. Executable Statements ..
Write (nout,*) 'H05ABF Example Program Results'
Write (nout,*)
! Skip headings in data file
Read (nin,*)
Read (nin,*)
! Read in the problem size
Read (nin,*) m, ip, nbest
! Read in the control parameters for the subset selection
Read (nin,*) mincr, mincnt, gamma, acc(1:2)
! Allocate memory required by the subset selection routine
mip = m - ip
Allocate (z(mip),a(max(nbest,m)),bz(mip,nbest),bscore(max(nbest,m)))
! Prepare the user workspace arrays IUSER and RUSER
Call prepare_user_arrays(m,iuser,ruser)
! Call the forward communication best subset routine
ifail = -1
Call h05abf(mincr,m,ip,nbest,la,bscore,bz,f,mincnt,gamma,acc,iuser, &
ruser,ifail)
If (ifail/=0 .And. ifail/=42) Then
! An error occurred
Go To 100
End If
! Titles
Write (nout,99999) ' Score Feature Subset'
Write (nout,99999) ' ----- --------------'
! Display the best subsets and corresponding scores. H05AAF returns a list
! of features excluded from the best subsets, so this is inverted to give
! the set of features included in each subset
Allocate (ibz(m),mask(m))
ibz(1:m) = (/(i,i=1,m)/)
Do i = 1, la
mask(1:m) = .True.
Do j = 1, mip
mask(bz(j,i)) = .False.
End Do
Write (nout,99998) bscore(i), pack(ibz,mask)
End Do
Write (nout,*)
If (ifail==42) Then
Write (nout,99997) nbest, &
' subsets of the requested size do not exist, only ', la, &
' are displayed.'
End If
Write (nout,99996) iuser(5), ' subsets evaluated in total'
100 Continue
99999 Format (1X,A)
99998 Format (1X,E12.5,100(1X,I5))
99997 Format (1X,I0,A,I0,A)
99996 Format (1X,I0,A)
End Program h05abfe