Program g02dafe ! G02DAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g02buf, g02daf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: aic, arsq, en, mult, rsq, rss, sw, tol Integer :: i, idf, ifail, ip, irank, ldq, ldx, & lwt, m, n Logical :: svd Character (1) :: mean, weight ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), q(:,:), & res(:), se(:), wk(:), wt(:), x(:,:), & y(:) Real (Kind=nag_wp) :: c(1), wmean(1) Integer, Allocatable :: isx(:) ! .. Intrinsic Procedures .. Intrinsic :: count, log, real ! .. Executable Statements .. Write (nout,*) 'G02DAF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) ! Read in the problem size Read (nin,*) n, m, weight, mean If (weight=='W' .Or. weight=='w') Then lwt = n Else lwt = 0 End If ldx = n Allocate (x(ldx,m),y(n),wt(lwt),isx(m)) ! Read in data 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(1:m)>0) If (mean=='M' .Or. mean=='m') Then ip = ip + 1 End If ldq = n Allocate (b(ip),cov((ip*ip+ip)/2),h(n),p(ip*(ip+ & 2)),q(ldq,ip+1),res(n),se(ip),wk(ip*ip+5*(ip-1))) ! Use suggested value for tolerance tol = 0.000001E0_nag_wp ! Fit general linear regression model ifail = -1 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) If (ifail/=0) Then If (ifail/=5) Then Go To 100 End If End If ! Calculate (weighted) total sums of squares, adjusted for mean if required ! If in G02DAF, an intercept is added to the regression by including a ! column of 1's in X, rather than by using the MEAN argument then ! MEAN = 'M' should be used in this call to G02BUF. ifail = 0 Call g02buf(mean,weight,n,1,y,n,wt,sw,wmean,c,ifail) ! Get effective number of observations (=N if there are no zero weights) en = real(idf+irank,kind=nag_wp) ! Calculate R-squared, corrected R-squared and AIC rsq = 1.0_nag_wp - rss/c(1) If (mean=='M' .Or. mean=='m') Then mult = (en-1.0E0_nag_wp)/(en-real(irank,kind=nag_wp)) Else mult = en/(en-real(irank,kind=nag_wp)) End If arsq = 1.0_nag_wp - mult*(1.0_nag_wp-rsq) aic = en*log(rss/en) + 2.0_nag_wp*real(irank,kind=nag_wp) ! Display results If (svd) Then Write (nout,99999) 'Model not of full rank, rank = ', irank Write (nout,*) End If Write (nout,99998) 'Residual sum of squares = ', rss Write (nout,99999) 'Degrees of freedom = ', idf Write (nout,99998) 'R-squared = ', rsq Write (nout,99998) 'Adjusted R-squared = ', arsq Write (nout,99998) 'AIC = ', aic Write (nout,*) Write (nout,*) 'Variable Parameter estimate ', 'Standard error' Write (nout,*) If (ifail==0) Then Write (nout,99997)(i,b(i),se(i),i=1,ip) Else Write (nout,99996)(i,b(i),i=1,ip) End If Write (nout,*) Write (nout,*) ' Obs Residuals H' Write (nout,*) Write (nout,99997)(i,res(i),h(i),i=1,n) 100 Continue 99999 Format (1X,A,I4) 99998 Format (1X,A,E12.4) 99997 Format (1X,I6,2E20.4) 99996 Format (1X,I6,E20.4) End Program g02dafe