Parsing a Model at Run Time

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: NLREG2                                              */
/*   TITLE: Parsing a Model at Run Time                         */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATRIX  SUGI6                                       */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: Rick Wicklin                UPDATE: SEP2013         */
/*     REF: Originally a full-screen editor by WMS              */
/*                                                              */
/****************************************************************/


proc iml;

/* _FUN and _DER are text strings that define model and deriv */
/* _parm contains parm names                                  */
/* _beta contains initial values for parameters               */
/* _k is the number of parameters                             */
start nlinit;
   _dep = "uspop";                      /* dependent variable */
   _fun = "a0*exp(a1*time)";    /* nonlinear regression model */
   /* deriv w.r.t. parameters */
   _der  = {"exp(a1*time)", "time*a0*exp(a1*time)"};
   _parm = {"a0", "a1"};               /* names of parameters */
   _beta = {3.9, 0};          /* initial guess for parameters */
   _k= nrow(_parm);                   /* number of parameters */
finish nlinit;

/* Generate the following modules at run time:                */
/* NLFIT: evaluate the model. After RUN NLFIT:                */
/*    _y contains response,                                   */
/*    _p contains predictor after call                        */
/*    _r contains residuals                                   */
/*    _sse contains sse                                       */
/* NLDERIV: evaluate derivs w.r.t params. After RUN NLDERIV:  */
/*    _x contains jacobian                                    */
start nlgen;
   call change(_fun, '*', '#', 0);  /* substitute '#' for '*' */
   call change(_der, '*', '#', 0);

   /* Write the NLFIT module at run time */
   call queue('start nlfit;');
   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;" );

   /* Write the NLDERIV function at run time */
   call queue('start nlderiv; free _NULL_; _x = ');
   do i=1 to _k;
      call queue("(", _der[i], ")||");
   end;
   call queue("_NULL_; finish;");

   call queue("resume;");    /* Pause to compile the functions */
   pause *;
finish nlgen;

/* Gauss-Newton nonlinear regression with Hartley step-halving */
start nlest;
   run nlfit;                /* f, r, and sse for initial beta */

   /* Gauss-Newton iterations to estimate parameters */
   do _iter=1 to 30 until(_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 nlfit;                         /* compute residual */
      _eps = abs((_lastsse-_sse)) / (_sse+1e-6);
      /* Hartley subiterations */
      do _subit=1 to 10 while(_sse>_lastsse);
         _delta = _delta/2;   /* halve the correction vector */
         _beta = _old+_delta; /* apply the halved correction */
         run nlfit;                        /* find sse et al */
      end;
      /* if no improvement after 10 halvings, exit iter loop */
      if _subit>10 then _eps=0;
   end;

   /* display table of results  */
   if _iter < 30 then do; /* convergence */
      _dfe = nrow(_y) - _k;
      _mse = _sse/_dfe;
      _std = sqrt(vecdiag(_xpxi)#_mse);
      _t = _beta/_std;
      _prob = 1 - cdf("F", _t#_t, 1, _dfe);
      print _beta[label="Estimate"] _std[label="Std Error"]
            _t[label="t Ratio"]     _prob[format=pvalue6.];
      print _iter[label="Iterations"] _lastsse[label="SSE"];
   end;
   else print "Convergence failed";
finish nlest;

/* main program: run nonlinear regression on data */
uspop = {3929, 5308, 7239, 9638, 12866, 17069, 23191, 31443, 39818,
         50155, 62947, 75994, 91972, 105710, 122775, 131669, 151325,
         179323, 203211}/1000;  /* US population, in thousands */
year = do(1790,1970,10)`;
time = year - 1790;

run nlinit; /* define strings that define the regression model */
run nlgen;  /* write modules that evaluate the model           */
run nlest;  /* compute param estimates, std errs, and p-values */