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;