PROGRAM g02defe ! G02DEF Example Program Text ! Mark 23 Release. NAG Copyright 2011. ! .. Use Statements .. USE nag_library, ONLY : g02ddf, g02def, nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: nin = 5, nout = 6 ! .. Local Scalars .. REAL (KIND=nag_wp) :: rss, rsst, tol INTEGER :: i, idf, ifail, ip, irank, ldq, lwt, & m, n LOGICAL :: svd CHARACTER (1) :: weight ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), p(:), q(:,:), se(:), & wk(:), wt(:), x(:) ! .. Executable Statements .. WRITE (nout,*) 'G02DEF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) ! Read in the problem size READ (nin,*) n, m, weight IF (weight=='W' .OR. weight=='w') THEN lwt = n ELSE lwt = 0 END IF ldq = n ALLOCATE (b(m),cov(m*(m+1)/2),p(m*(m+2)),q(ldq,m+1),se(m),wk(m*m+5*m), & wt(n),x(n)) ! Read in the dependent variable, Y, and store in first column of Q READ (nin,*) q(1:n,1) ! Read in weights IF (lwt>0) THEN READ (nin,*) wt(1:n) END IF ! Use suggested value for tolerance tol = 0.000001E0_nag_wp ! Loop over each of the supplied variables ip = 0 U_LP: DO READ (nin,*,IOSTAT=ifail) x(1:n) IF (ifail/=0) THEN EXIT U_LP END IF ! Add the new variable to the model ifail = -1 CALL g02def(weight,n,ip,q,ldq,p,wt,x,rss,tol,ifail) IF (ifail/=0) THEN IF (ifail==3) THEN WRITE (nout,99999) ' * Variable ', ip, & ' is linear combination of previous columns' WRITE (nout,*) ' so it has not been added' WRITE (nout,*) CYCLE U_LP ELSE GO TO 20 END IF END IF ip = ip + 1 WRITE (nout,99999) 'Variable ', ip, ' added' ! Get G02DDF to recalculate RSS rsst = 0.0E0_nag_wp ! Calculate the parameter estimates ifail = 0 CALL g02ddf(n,ip,q,ldq,rsst,idf,b,se,cov,svd,irank,p,tol,wk,ifail) IF (svd) THEN WRITE (nout,*) 'Model not of full rank' WRITE (nout,*) END IF WRITE (nout,99998) 'Residual sum of squares = ', rsst WRITE (nout,99999) '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) WRITE (nout,*) END DO U_LP 20 CONTINUE 99999 FORMAT (1X,A,I0,A) 99998 FORMAT (1X,A,E13.4) 99997 FORMAT (1X,I6,2E20.4) END PROGRAM g02defe