Derivative-Free Regression

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: DUD                                                 */
 /*   TITLE: Derivative-Free Regression                          */
 /* PRODUCT: IML                                                 */
 /*    KEYS: MATRIX REGR                                         */
 /*   PROCS: IML NLIN                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: LWB                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC: CONVERTED FROM MATRIX TO IML USING MATIML           */
 /*                                                              */
 /****************************************************************/
*-----------NONLINEAR LEAST SQUARES USING THE DUD METHOD------------*
| SEE MARY RALSTON AND ROBERT JENNRICH, TECHNOMETRICS 20(1) FEB 78  |
| PP 7-14, "DUD, A DERIVATIVE-FREE ALGORITHM FOR NONLINEAR LEAST    |
| SQUARES."   THE DATA IN THIS EXAMPLE IS THE POPULATION OF THE US  |
| 1790 TO 1970.  THE MODEL IS A SIMPLE EXPONENTIAL GROWTH CURVE.    |
*-----------------------------------------------------------JPSALL--*;

TITLE1 'SAMPLE PROGRAM:  DUD';
DATA USPOP;
   INPUT POP@@;
   RETAIN YEAR 1780;
   YEAR=YEAR+10;
   CARDS;
     3.929 5.308 7.239 9.638 12.866 17.069 23.191 31.443 39.818
     50.155 62.947 75.994 91.972 105.710 122.775 131.669 151.325
     179.323 203.211
;
RUN;
PROC PRINT;
     TITLE2 'GROWTH OF POPULATION IN THE US';
RUN;
PROC NLIN DATA=USPOP;
   PARMS C0=3.9 C1=.022;
   Z=EXP(C1*(YEAR-1790));
   MODEL POP=C0*Z;
   DER.C0=Z;
   DER.C1=(YEAR-1790)*C0*Z;
   TITLE2 'NONLINEAR LEAST SQUARES USING GAUSS-NEWTON';
RUN;

PROC IML;
  TITLE2 'NONLINEAR LEAST SQUARES USING THE DUD METHOD';

  *---ROUTINE TO COMPUTE RESPONSE FUNCTION AND SS---;
  START FUN;
     F= B[1,]# EXP( B[2,]#YEAR);
     SS= SSQ(Y-F);
  FINISH;

  START DUD;
     *---ROUTINE TO GET LEAST SQUARES ESTIMATES USING THE DUD METHOD---;
     *---SPECIFIED BY THE USER:
          ---B0   INITIAL PARAMETER ESTIMATES
          ---Y    RESPONSE VARIABLE
          ---FUN  ROUTINE TO EVALUATE F (RESPONSE) AND SS ;

     *---DIMENSIONS-------------------------;
     P= NROW(B0);                       *---NUMBER OF PARAMETERS;
     P2=2:P;                          *---INDICES FROM 2 TO P;
     N= NROW(Y);                        *---NUMBER OBSERVATIONS;

     *---GET STARTING VALUES----------------;
     DO I=1 TO P;                     *---P STARTING TRIALS;
        B=B0; B[I,]= B0[I,]#1.1;      *---PERTURB THE ITH PARM;
        BB=BB||B;                       *---STORE PERTURBED PARM;
        RUN FUN;                        *---EVALUATE RESPONSE;
        SSS=SSS||SS;                    *---STORE RESULTS;
        FF=FF||F;
     END;

     *---ORDER STARTING VALUES--------------;
     R= RANK(-SSS);                     *---IN DESCENDING ORDER;
     Z=SSS; SSS[,R]=Z;                  *---REORDER IN ORDER BY     ;
     Z=FF;                              *---RESIDUAL SS             ;
     FF[,R]=Z;
     Z=BB;
     BB[,R]=Z;

     *---PREPARE FOR ITERATIONS-------------;
     B=B0;
     RUN FUN;
     CRIT=1;
     PRINT "STARTING VALUES AND SSE";
     PRINT B0, SS;

     *---ITERATE UNTIL CONVERGED------------;
     DO IT=1 TO 20 WHILE(CRIT>1E-5);  *---ITERATE;
        OLDB=B;                             *---SAVE;
        OLDSS=SS;
        OLDF=F;
        J=FF-F* J(1,P);
        ALPHA= SOLVE(J`*J,J`*(Y-F));
        DP=BB-B* J(1,P);                  *---CHANGE IN PARMS;
        DELTA=DP*ALPHA;                     *---CORRECTION VECTOR;
        B=B+DELTA;                          *---APPLY THE CORRECTION;
        RUN FUN;                       *---EVALUATE RESPONSE AND SSE;

        *---STEP SHORTENING IF SSE DOES NOT IMPROVE-----;
        DO IS=1 TO 8 WHILE(SS>OLDSS);
           DELTA=DELTA#.5;           *---HALF OF THE PREVIOUS STEP;
           B=OLDB+DELTA;
           RUN FUN;
        END;

        *---TO THE NEXT ITERATION---;
        BB= BB[,P2]||OLDB;
        SSS= SSS[,P2]||OLDSS;
        FF= FF[,P2]||OLDF;
        CRIT= ABS(SS-OLDSS)/OLDSS;     *---CRITERION FOR CONVERGENCE;
        SHOW=SHOW//(IT||IS||CRIT||SS||(B`));
    END;

     *---PRINT OUT RESULTS------------------;
     SNAME={'ITER' 'HALVINGS' 'CRIT' 'SS' 'B1' 'B2' 'B3' 'B4' 'B5'};
     PRINT SHOW[COLNAME=SNAME L="ITERATIONS"];
     PRINT B[L="PARAMETER ESTIMATES"];
     S=(SS/(N-P))# DP* INV(J`*J)*DP`;
     PRINT S[L="COVARIANCE OF PARAMETER ESTIMATES"];
     RESID=Y-F;
     PRINT RESID[L="RESIDUALS"];
  FINISH DUD;

  *---PREPARE THE DATA, INITIAL VALUES, CALL DUD---;
  USE USPOP ;
  READ ALL INTO X ;
  Y= X[,1];
  YEAR= X[,2]-1790;
  B0={3.9 .022}`;

  RUN DUD;
quit;
title;title2;