Program g02dcfe ! G02DCF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g02daf, g02dcf, g02ddf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: rss, tol, wt, y Integer :: i, idf, ifail, ip, irank, ix, ldq, & ldxm, lwt, m, n Logical :: svd Character (1) :: mean, update, weight ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: b(:), cov(:), h(:), p(:), q(:,:), & res(:), se(:), wk(:), wtm(:), x(:), & xm(:,:), ym(:) Integer, Allocatable :: isx(:) ! .. Intrinsic Procedures .. Intrinsic :: count ! .. Executable Statements .. Write (nout,*) 'G02DCF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) n, m, weight, mean If (weight=='W' .Or. weight=='w') Then lwt = n Else lwt = 0 End If ldxm = n Allocate (xm(ldxm,m),ym(n),wtm(lwt),isx(m),x(m)) ! Read in data If (lwt>0) Then Read (nin,*)(xm(i,1:m),ym(i),wtm(i),i=1,n) Else Read (nin,*)(xm(i,1:m),ym(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 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 initial model using G02DAF ifail = 0 Call g02daf(mean,weight,n,xm,ldxm,m,isx,ip,ym,wtm,rss,idf,b,se,cov,res, & h,q,ldq,svd,irank,p,tol,wk,ifail) ! Display results from initial model fit Write (nout,*) 'Results from initial model fit using G02DAF' If (svd) Then Write (nout,*) Write (nout,*) 'Model not of full rank' End If 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) ! Updating data is held in X consecutively ix = 1 ! Add or delete observations supplied in the data file u_lp: Do Read (nin,*,Iostat=ifail) update If (ifail/=0) Then Exit u_lp End If If (lwt>0) Then Read (nin,*) x(1:m), y, wt Else Read (nin,*) x(1:m), y End If ! Update the regression ifail = 0 Call g02dcf(update,mean,weight,m,isx,q,ldq,ip,x,ix,y,wt,rss,wk,ifail) ! Display titles and update observation count Write (nout,*) Select Case (update) Case ('a','A') Write (nout,*) 'Results from adding an observation using G02DCF' n = n + 1 Case ('d','D') Write (nout,*) 'Results from dropping an observation using G02DCF' n = n - 1 Case Default Write (nout,*) 'Unknown update flag read in from data file' Go To 100 End Select ! Recalculate the parameter estimates etc ifail = 0 Call g02ddf(n,ip,q,ldq,rss,idf,b,se,cov,svd,irank,p,tol,wk,ifail) ! Display updated results 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) End Do u_lp 100 Continue 99999 Format (1X,A,E12.4) 99998 Format (1X,A,I4) 99997 Format (1X,I6,2E20.4) End Program g02dcfe