PROGRAM g02dkfe ! G02DKF Example Program Text ! Mark 23 Release. NAG Copyright 2011. ! .. Use Statements .. USE nag_library, ONLY : g02daf, g02dkf, nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: nin = 5, nout = 6 ! .. Local Scalars .. REAL (KIND=nag_wp) :: rss, tol INTEGER :: i, iconst, idf, ifail, ip, irank, & ldc, ldq, ldx, lwk, lwt, m, n, tlwk LOGICAL :: svd CHARACTER (1) :: mean, weight ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: b(:), c(:,:), cov(:), h(:), p(:), & q(:,:), res(:), se(:), wk(:), wt(:), & x(:,:), y(:) INTEGER, ALLOCATABLE :: isx(:) ! .. Intrinsic Functions .. INTRINSIC count ! .. Executable Statements .. WRITE (nout,*) 'G02DKF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) n, m, weight, mean ! Read in data IF (weight=='W' .OR. weight=='w') THEN lwt = n ELSE lwt = 0 END IF ldx = n ALLOCATE (x(ldx,m),isx(m),y(n),wt(lwt)) IF (lwt>0) THEN READ (nin,*) (x(i,1:m),y(i),wt(i),i=1,n) ELSE READ (nin,*) (x(i,1:m),y(i),i=1,n) END IF ! Read in variable inclusion flags READ (nin,*) isx(1:m) ! Calculate IP ip = count(isx>0) IF (mean=='M' .OR. mean=='m') THEN ip = ip + 1 END IF lwk = 5*(ip-1) + ip*ip ldq = n ldc = ip ALLOCATE (b(ip),se(ip),cov(ip*(ip+1)/2),res(n),h(n),q(ldq,ip+1),p(2*ip+ & ip*ip),wk(lwk),c(ldc,ip)) ! Use suggested value for tolerance tol = 0.000001E0_nag_wp ! Fit general linear regression model ifail = 0 CALL g02daf(mean,weight,n,x,ldx,m,isx,ip,y,wt,rss,idf,b,se,cov,res,h,q, & ldq,svd,irank,p,tol,wk,ifail) ! Display the initial estimates WRITE (nout,*) 'Estimates from G02DAF' WRITE (nout,*) WRITE (nout,99999) 'Residual sum of squares = ', rss WRITE (nout,99998) 'Degrees of freedom = ', idf WRITE (nout,*) WRITE (nout,*) 'Variable Parameter estimate Standard error' WRITE (nout,*) WRITE (nout,99997) (i,b(i),se(i),i=1,ip) ! Read in constraints iconst = ip - irank READ (nin,*) (c(i,1:iconst),i=1,ip) tlwk = 2*ip*ip + ip*iconst + 2*iconst*iconst + 4*iconst IF (tlwk>lwk) THEN ! Reallocate workspace DEALLOCATE (wk) ALLOCATE (wk(tlwk)) END IF ! Re-estimate the parameters given the constraints ifail = 0 CALL g02dkf(ip,iconst,p,c,ldc,b,rss,idf,se,cov,wk,ifail) ! Display the results WRITE (nout,*) WRITE (nout,*) 'Estimates from G02DKF using constraints' WRITE (nout,*) WRITE (nout,*) 'Variable Parameter estimate Standard error' WRITE (nout,*) WRITE (nout,99997) (i,b(i),se(i),i=1,ip) 99999 FORMAT (1X,A,E13.3) 99998 FORMAT (1X,A,I4) 99997 FORMAT (1X,I6,2E20.3) END PROGRAM g02dkfe