$ontext Non-Linear Least Squares Regression example Erwin Kalvelagen, may 2008 Reference: http://www.statsci.org/data/general/brunhild.html Blood Sulfate in a Baboon Named Brunhilda -------------------------------------------------------------------------------- Description The observed responses are Geiger counter counts (times 10-4) used to measure the amount of radioactively tagged sulfate drug in the blood of a baboon named Brunhilda after an injection of the drug. -------------------------------------------------------------------------------- Variable Description -------------------------------------------------------------------------------- Hours Time in hours since injection Sulfate Geiger counter counts × 10^-4 -------------------------------------------------------------------------------- Source Jennrich, R. I., and Bright, P. B. (1976). Fitting systems of linear differential equations using computer generated exact derivatives. Technometrics 18, 385-392. Jennrich, R. I. (1995). An Introduction to Computational Statistics. Prentice-Hall, Englewood Cliffs, New Jersey. Section 8.2.1. Analysis The sulfate decay curve is not well fitted by simple exponential decay, Sulfate = a exp(-b Hours). The model Suflate = a Hours^b is not bad, as the plot shows. Successive observations will be positively correlated. In this case the counts arrive according to a Poisson process. One could obtain a goodness of fit test for the fitted curve by comparing the errors to Poisson variation. $offtext set i /i1*i21/; table data(i,*) Hours Sulfate i1 2 15.11 i2 4 11.36 i3 6 9.77 i4 8 9.09 i5 10 8.48 i6 15 7.69 i7 20 7.33 i8 25 7.06 i9 30 6.7 i10 40 6.43 i11 50 6.16 i12 60 5.99 i13 70 5.77 i14 80 5.64 i15 90 5.39 i16 110 5.09 i17 130 4.87 i18 150 4.6 i19 160 4.5 i20 170 4.36 i21 180 4.27 ; parameters Hours(i) 'Time in hours since injection' Sulfate(i) 'Geiger counter counts × 10^-4' ; Hours(i) = data(i,'Hours'); Sulfate(i) = data(i,'Sulfate'); variables a 'parameter to be estimated' b 'parameter to be estimated' sse 'sum of squared errors' ; equation fit1(i) 'equation to fit sulfate=a*exp(-b*hours)' fit2(i) 'equation to fit sulfate=a*Hours^b' fit3(i) 'linearized equation to fit log(sulfate)=a+b*log(hours)' sumsq 'dummy objective' ; sumsq.. sse =n= 0; fit1(i).. Sulfate(i) =e= a * exp(-b * Hours(i)); fit2(i).. Sulfate(i) =e= a * Hours(i) ** b ; fit3(i).. log(Sulfate(i)) =e= a + b * log(Hours(i)) ; options nlp=nls, lp=ls; models m1 /sumsq,fit1/ m2 /sumsq,fit2/ m3 /sumsq,fit3/ ; solve m1 minimizing sse using nlp; solve m2 minimizing sse using nlp; display a.l,b.l; solve m3 minimizing sse using lp; scalar a2; a2 = exp(a.l); display a2,b.l;