 |
|
| General Statistics Examples |
Example 9.15 Full-Screen Nonlinear Regression
This example shows how to build a menu system that enables you to perform nonlinear regression from a menu. Six modules are stored on an IML storage disk. After you have stored them, use this example to try out the system. First, invoke IML and set up some sample data in memory, in this case the population of the U.S. from 1790 to 1970. Then invoke the module NLIN, as follows:
reset storage='nlin';
load module=_all_;
uspop = {3929, 5308, 7239, 9638, 12866, 17069, 23191, 31443,
39818, 50155, 62947, 75994, 91972, 105710, 122775, 131669,
151325, 179323, 203211}/1000;
year=do(1790,1970,10)`;
time=year-1790;
print year time uspop;
run nlin;
A menu similar to the following menu appears. The entry fields are shown by underscores here, but the underscores become blanks in the real session.
Nonlinear Regression Response function: ______________________________________________ Predictor function: _____________________________________________ Parameter Value Derivative : ________ ___________ __________________________________________ : ________ ___________ __________________________________________ : ________ ___________ __________________________________________ : ________ ___________ __________________________________________ : ________ ___________ __________________________________________ : ________ ___________ __________________________________________
Enter an exponential model and fill in the response and predictor expression fields. For each parameter, enter the name, initial value, and derivative of the predictor with respect to the parameter. Here are the populated fields:
Nonlinear Regression Response function: uspop_________________________________________ Predictor function: a0*exp(a1*time)______________________________ Parameter Value Derivative : a0______ ________3.9 exp(a1*time)_______________________________ : a1______ __________0 time*a0*exp(a1*time)_______________________ : ________ ___________ ___________________________________________ : ________ ___________ ___________________________________________ : ________ ___________ ___________________________________________ : ________ ___________ ___________________________________________
Now press the SUBMIT key. The model compiles, the iterations start blinking on the screen, and when the model has converged, the estimates are displayed along with their standard errors, 
test, and significance probability.
To modify and rerun the model, submit the following command:
run nlrun;
Here is the program that defines and stores the modules of the system.
/* Full-Screen Nonlinear Regression */ /* Six modules are defined, which constitute a system for */ /* nonlinear regression. The interesting feature of this */ /* system is that the problem is entered in a menu, and both */ /* iterations and final results are displayed on the same */ /* menu. */ /* */ /* Run this source to get the modules stored. Examples */ /* of use are separate. */ /* */ /* Caution: this is a demonstration system only. It does not */ /* have all the necessary safeguards in it yet to */ /* recover from user errors or rough models. */ /* Algorithm: */ /* Gauss-Newton nonlinear regression with step-halving. */ /* Notes: program variables all start with nd or _ to */ /* minimize the problems that would occur if user variables */ /* interfered with the program variables. */
/* Gauss-Newton nonlinear regression with Hartley step-halving */
/*---Routine to set up display values for new problem---*/
start nlinit;
window nlin rows=15 columns=80 color='green'
msgline=_msg cmndline=_cmnd
group=title +30 'Nonlinear Regression' color='white'
group=model / @5 'Response function:' color='white'
+1 nddep $55. color='blue'
/ @5 'Predictor function:' color='white'
+1 ndfun $55. color='blue'
group=parm0 // @5 'Parameter' color='white' @15 'Value'
@30 'Derivative'
group=parm1 // @5 'Parameter' color='white' @15 'Value'
group=parm2 // @5 'Parameter' color='white' @19 'Estimate'
@33 'Std Error'
@48 'T Ratio'
@62 'Prob>|T|'
group=parminit /@3 ':' color='white'
@5 ndparm $8. color='blue'
@15 ndbeta best12. @30 ndder $45.
group=parmiter / @5 _parm color='white'
@15 _beta best12. color='blue'
group=parmest / @5 _parm color='white'
@15 _beta best12. color='blue'
@30 _std best12.
@45 _t 10.4
@60 _prob 10.4
group=sse // @5 'Iteration =' color='white' _iter 5. color='blue'
' Stephalvings = ' color='white' _subit 3. color='blue'
/ @5 'Sum of Squares Error =' color='white' _sse best12.
color='blue';
nddep=cshape(' ',1,1,55,' ');
ndfun=nddep;
nd0=6;
ndparm=repeat(' ',nd0,1);
ndbeta=repeat(0,nd0,1);
ndder=cshape(' ',nd0,1,55,' ');
_msg='Enter New Nonlinear Problem';
finish nlinit; /* Finish module NLINIT */
/* Main routine */
start nlin;
run nlinit; /* initialization routine */
run nlrun; /* run routine */
finish nlin;
/* Routine to show each iteration */
start nliter;
display nlin.title noinput,
nlin.model noinput,
nlin.parm1 noinput,
nlin.parmiter repeat noinput,
nlin.sse noinput;
finish nliter;
/* Routine for one run */
start nlrun;
run nlgen; /* generate the model */
run nlest; /* estimate the model */
finish nlrun;
/* Routine to generate the model */
start nlgen;
/* Model definition menu */
display nlin.title, nlin.model, nlin.parm0, nlin.parminit repeat;
/* Get number of parameters */
t=loc(ndparm=' ');
if nrow(t)=0 then
do;
print 'no parameters';
stop;
end;
_k=t[1] -1;
/* Trim extra rows, and edit '*' to '#' */
_dep=nddep; call change(_dep,'*','#',0);
_fun=ndfun; call change(_fun,'*','#',0);
_parm=ndparm[1:_k,];
_beta=ndbeta[1:_k,];
_der=ndder [1:_k,];
call change(_der,'*','#',0);
/* Construct nlresid module to split up parameters and */
/* compute model */
call queue('start nlresid;');
do i=1 to _k;
call queue(_parm[i] ,"=_beta[",char(i,2),"] ;");
end;
call queue("_y = ",_dep,";",
"_p = ",_fun,";",
"_r = _y-_p;",
"_sse = ssq(_r);",
"finish;" );
/* Construct nlderiv function */
call queue('start nlderiv; _x = ');
do i=1 to _k;
call queue("(",_der[i] ,")#repeat(1,nobs,1)||");
end;
call queue(" nlnothin; finish;");
/* Pause to compile the functions */
call queue("resume;");
pause *;
finish nlgen; /* Finish module NLGEN */
/* Routine to do estimation */
start nlest;
/* Modified Gauss-Newton Nonlinear Regression */
/* _parm has parm names */
/* _beta has initial values for parameters */
/* _k is the number of parameters */
/* after nlresid: */
/* _y has response, */
/* _p has predictor after call */
/* _r has residuals */
/* _sse has sse */
/* after nlderiv */
/* _x has jacobian */
/* */
eps=1;
_iter = 0;
_subit = 0;
_error = 0;
run nlresid; /* f, r, and sse for initial beta */
run nliter; /* print iteration zero */
nobs = nrow(_y);
_msg = 'Iterating';
/* Gauss-Newton iterations */
do _iter=1 to 30 while(eps>1e-8);
run nlderiv; /* subroutine for derivatives */
_lastsse=_sse;
_xpxi=sweep(_x`*_x);
_delta=_xpxi*_x`*_r; /* correction vector */
_old = _beta; /* save previous parameters */
_beta=_beta+_delta; /* apply the correction */
run nlresid; /* compute residual */
run nliter; /* print iteration in window */
eps=abs((_lastsse-_sse))/(_sse+1e-6);
/* convergence criterion */
/* Hartley subiterations */
do _subit=1 to 10 while(_sse>_lastsse);
_delta=_delta*.5; /* halve the correction vector */
_beta=_old+_delta; /* apply the halved correction */
run nlresid; /* find sse et al */
run nliter; /* print subiteration in window */
end;
if _subit>10 then
do;
_msg = "did not improve after 10 halvings";
eps=0; /* make it fall through iter loop */
end;
end;
/* print out results */
_msg = ' ';
if _iter>30 then
do;
_error=1;
_msg = 'convergence failed';
end;
_iter=_iter-1;
_dfe = nobs-_k;
_mse = _sse/_dfe;
_std = sqrt(vecdiag(_xpxi)#_mse);
_t = _beta/_std;
_prob= 1-probf(_t#_t,1,_dfe);
display nlin.title noinput,
nlin.model noinput,
nlin.parm2 noinput,
nlin.parmest repeat noinput,
nlin.sse noinput;
finish nlest; /* Finish module NLEST */
/* Store the modules to run later */
reset storage='nlin';
store module=_all_;
|
|
Copyright © SAS Institute, Inc. All Rights Reserved.