{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 287 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 288 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 291 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 292 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 293 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 1 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "File: depart\\math\\maple \\misc\\nlfit.mws Date: 2-dec-1999" }}}{EXCHG {PARA 18 "" 0 "" {TEXT -1 54 "Nonlinear Fitting Procedure for 1 Independent Variable" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 15 "This procedure " }{TEXT 292 5 "nlf it" }{TEXT -1 330 ", implemented by David Holmgren, is an implementati on of an algorithm by Jefferys to estimate non-linear parameters with \+ a method due to Levenberg and Marquardt. The user must specify extern ally a formula that constitutes an objective function in terms of a si ngle independent variable and parameters. Within a second procedure \+ " }{TEXT 293 2 "gr" }{TEXT -1 316 ", derivatives of the formula with r espect to the parameters are calculated algebraically, thus making use of Maple's symbolic capability within an essentially numerical applic ation of fitting data, and also these derivatives and the formula itse lf are evaluated with current values of parameters as a fit progresses ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "restart: with(linalg): Digits := 20:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4477 "nlfit : = proc(f, X, x::vector, y::vector, wt::vector, P::\{name,list(name)\}, p::vector, tol)\n local gradf,nobs,npar,k,yck,xk,r,a,drv,niter,j,n,co rr,alpha,beta,\n cvm,b,p1,errors,dp,res,nu,res1,sigma1,reso;\n # In itialize.\n if type(P,name) then RETURN(nlfit(f,X,x,y,wt,[seq(P[k],k=1 ..vectdim(p))],p,tol)) fi;\n gradf := linalg[grad](f,P);\n reso := 1.0 *10^9;\n nobs := vectdim(x);\n npar := vectdim(p);\n r := vector(nobs) ;\n a := matrix(nobs,npar,0.0);\n corr := matrix(npar,npar,0.0); \n dr v := vector(npar);\n p1 := vector(npar);\n errors := vector(npar);\n n iter := 50;\n nu := 1000;\n # Begin iterations.\nfor n from 1 to niter \n do\n print(` iteration No. `,n);\n# Compute initial fit and residu als (observed - calculated)...\n for k from 1 to nobs\n do\n xk := x[k];\n yck,drv := gr(f,gradf,X,P,npar,xk,p); r[k] := evalf(sqrt(wt [k])*(y[k] - yck));\n# Compute gradient (drv) and set up equations of \+ condition (a)...\n# Note use of weighting here.\n for j from 1 to np ar\n do\n a[k,j] := sqrt(wt[k])*drv[j];\n od:\n od:\n# sum of \+ squares of residuals...\n res1 := evalf(dotprod(r,r));\n print(` initi al SSR = `,res1);\n# Compute matrix of normal equations (alpha)...\n a lpha := evalf(evalm(transpose(a) &* a));\n# Compute RHS vector (beta) \+ for normal equations...\n beta := evalf(evalm(transpose(a) &* r));\n# \+ Note use of L-M factor nu in this next loop...\n# This increases as on e approaches solution.\n for j from 1 to npar\n do\n alpha[j,j] := ev alf(alpha[j,j]*(1. + 1./nu));\n od:\n# Solve the system of normal equa tions for parameter corrections...\n dp := evalf(linsolve(alpha,beta)) ;\n# Compute sum of squares of residuals at new point p1 = p + dp...\n p1 := evalf(evalm(p + dp));\n for k from 1 to nobs\n do\n xk := x[k] ;\n yck,drv := gr(f,gradf,X,P,npar,xk,p1); r[k] := evalf(sqrt(wt[k])* (y[k]-yck));\n od:\n res := evalf(dotprod(r,r));\n print(`new SSR = `, res);\n# L-M step: modify nu according to whether the sum of\n# square s increases or decreases after this correction...\n# wrong direction ( elements not corrected)...\n if res1 <= res then\n nu := nu/10.;\n el se\n# Right direction (correct elements)...\n nu := nu*10.;\n p := e valf(evalm(p + dp));\n fi:\n# test for convergence...\n if abs(res - r eso) <= tol then\n break\n fi;\n reso := res; \+ \nod:\n# covariance matr ix...\n cvm := evalf(inverse(alpha));\n# Standard error of a single ob servation of unit weight...\n sigma1 := sqrt(res/(nobs-npar));\n# erro rs of parameters ...\n for j from 1 to npar\n do\n errors[j] := evalf (sigma1*sqrt(cvm[j,j]));\n od:\nfor j from 1 to npar \+ \+ do \+ \+ for k from 1 to npar \+ \+ do \+ \+ corr[j,k] := cvm[j,k]/sqrt(cvm[j,j]*cvm[k,k]); \+ \+ od: \+ od: \+ \+ # print res ults...\n# vector of corrections...\n print(` vector of corrections `) ;\n print(dp);\n# final set of orbital elements...\n print(` final val ues of parameters `);\n print(p);\n print(` estimated standard errors \+ `);\n print(errors);\n# standard deviation of one observation of unit \+ weight...\n print(` standard deviation of fit = `,sigma1);\nprint(` ma trix of correlation coefficients `); \+ prin t(corr); \+ end: \+ \+ \ngr := proc(f,G,x,a,n,xval,p::vector) \n # procedure to evaluate derivatives and formula \nlocal k,y,i,drv,point;\npoint := \+ \{x=xval, seq(a[i]=p[i],i=1..n)\};\ndrv:=vector(n);\nfor k from 1 to n do \n drv[k] := evalf(eval(G[k],point)) \nod:\ny := evalf(eval(f,po int)); \nRETURN(y,drv);\nend:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 183 "To illustrate use of this procedure we prepare data in two sets. Her e is the first objective function or model, which is simply a cosine f unction of variable x with four parameters, " }{XPPEDIT 18 0 "a[1];" " 6#&%\"aG6#\"\"\"" }{TEXT -1 4 " -- " }{XPPEDIT 18 0 "a[4];" "6#&%\"aG6 #\"\"%" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "f := a[1] + a[2]*cos( a[3]*x + a[4] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 23 "We select 25 values of " }{TEXT 256 1 "x" }{TEXT -1 65 " \+ at which to calculate this function, which we store as a vector " } {TEXT 257 2 "xx" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "xx := vector(25): for j from 1 to 25 do xx[j] := 0.1+0.05*j od ; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "For purpose of test we sele ct values of parameters " }{XPPEDIT 18 0 "a[1];" "6#&%\"aG6#\"\"\"" } {TEXT -1 4 " -- " }{XPPEDIT 18 0 "a[4];" "6#&%\"aG6#\"\"%" }{TEXT -1 29 ", which we store as a vector " }{TEXT 258 1 "p" }{TEXT -1 66 ", co ntaining the number of parameters and their values as a list. " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "p := vector(4,[3,5,7,2]):" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 25 "We calculate as a vector " } {TEXT 259 1 "y" }{TEXT -1 58 " the corresponding 25 values of the func tion to be fitted." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 101 "y := vector(25): \nfor j from 1 to 25 do y[j] := evalf(eval(f, \{x=xx[j], \+ seq(a[k]=p[k], k=1..4)\})) od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 191 "We provide as a vector w the weights of these values; in general \+ such weights are evaluated from standard deviations of the correspondi ng values of y, but in this case we assume unit weights." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "w := vector(25): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "for j from 1 to 25 do w[j] := 1.0 o d:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "The same effect is achieved in this case with a single statement." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "w:=vector(25, 1):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "We provide as a vector initial estimates of parameters. " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "p2 := vector(4, [3.3,4.8,7.1 ,2.5]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 " In invoking proce dure " }{TEXT 294 5 "nlfit" }{TEXT -1 12 " hereunder, " }{TEXT 265 1 " f" }{TEXT -1 30 " is the name of the function, " }{TEXT 266 1 "X" } {TEXT -1 69 " is the name of the independent variable appearing in thi s function, " }{TEXT 267 1 "x" }{TEXT -1 36 " is the array of observed values of " }{TEXT 268 1 "X" }{TEXT -1 2 ", " }{TEXT 269 1 "y" } {TEXT -1 72 " is the array of observed values of the dependent variabl e or function, " }{TEXT 270 2 "wt" }{TEXT -1 102 " is an array of weig hts (usually inverse of standard deviations of each value of dependent variable), " }{TEXT 271 1 "P" }{TEXT -1 52 " is the name of parameter s as a vectorial quantity, " }{TEXT 272 1 "p" }{TEXT -1 63 " is a vect or of initial estimates of values of parameters, and " }{TEXT 273 3 "t ol" }{TEXT -1 110 " is a tolerance that serves as a criterion of conve rgence of the solution. For the given function fun above, " }{TEXT 280 1 "f" }{TEXT -1 3 " = " }{TEXT 281 3 "fun" }{TEXT -1 2 ", " } {TEXT 282 1 "X" }{TEXT -1 3 " = " }{TEXT 283 1 "x" }{TEXT -1 2 ", " } {TEXT 284 1 "y" }{TEXT -1 54 " = whatever name of array provided -- in this case y, " }{TEXT 285 2 "wt" }{TEXT -1 51 " = whatever name of ar ray provided -- in this case " }{TEXT 286 1 "w" }{TEXT -1 2 ", " } {TEXT 287 1 "P" }{TEXT -1 3 " = " }{TEXT 288 1 "a" }{TEXT -1 2 ", " } {TEXT 289 1 "p" }{TEXT -1 79 " = whatever name is provided for a vecto r of initial estimates -- in this case " }{TEXT 290 2 "p2" }{TEXT -1 87 ". For a rigorous first test we set a tolerance to correspond to 1 8 significant digits." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "to l := 10^(-(Digits - 2));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "With \+ these settings we initiate activity with procedure " }{TEXT 274 5 "nlf it" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "nlfit (f,x,xx,y,w,a,p2,tol);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "According to output above, aft er five iterations the sum of squared residuals decreases from 115 to \+ " }{XPPEDIT 18 0 "10^(-34);" "6#)\"#5,$\"#M!\"\"" }{TEXT -1 62 ". In \+ the final iteration the indicated corrections, of order " }{XPPEDIT 18 0 "10^(-12);" "6#)\"#5,$\"#7!\"\"" }{TEXT -1 130 ", in the display ed vector are applied to corresponding parameters that are altered fro m their initial values, provided in vector " }{TEXT 260 2 "p2" }{TEXT -1 110 ", to final values displayed above with their estimated standar d errors. The standard deviation of the fit, 3 " }{XPPEDIT 18 0 "10^( -18);" "6#)\"#5,$\"#=!\"\"" }{TEXT -1 80 ", reflects the magnitude of \+ the tolerance that is set before invoking procedure " }{TEXT 261 5 "nl fit" }{TEXT -1 54 ". The final values of parameters differ by at most 2 " }{XPPEDIT 18 0 "10^(-18);" "6#)\"#5,$\"#=!\"\"" }{TEXT -1 172 " f rom the values used to create the set of data fitted to the objective \+ function; because the data are essentially perfect, within the error o f rounding inherent in use of " }{TEXT 262 5 "float" }{TEXT -1 38 "s u sed to form the values of function " }{TEXT 263 1 "f" }{TEXT -1 11 " i n vector " }{TEXT 264 1 "y" }{TEXT -1 161 ", the procedure reproduces \+ essentially perfectly those values of parameters. The matrix of coeff icients of correlation between parameters shows that parameters " } {XPPEDIT 18 0 "a[1];" "6#&%\"aG6#\"\"\"" }{TEXT -1 4 " -- " }{XPPEDIT 18 0 "a[3];" "6#&%\"aG6#\"\"$" }{TEXT -1 100 " are essentially entirel y uncorrelated, as their coefficients are near zero; in contrast param eters " }{XPPEDIT 18 0 "a[3];" "6#&%\"aG6#\"\"$" }{TEXT -1 5 " and " } {XPPEDIT 18 0 "a[4];" "6#&%\"aG6#\"\"%" }{TEXT -1 91 " are moderately \+ closely anti-correlated, as their coefficient of correlation is about \+ -0.9." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 154 " In a normal course of events error of measurements inevitably arises in an experiment; w e simulate this phenomenon by adding to each generated value " }{TEXT 275 1 "y" }{TEXT -1 13 " of function " }{TEXT 276 1 "f" }{TEXT -1 169 " a random number uniformly distributed in a range [-0.25 .. 0.25]. A function to generate random integers uniformly distributed in a range [-500..500] is the following:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "ran := rand(-500..500):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "We can test this generator by printing some values." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "for j to 25 do ran()/2000.0 od;" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "With this generator of random num bers to add " }{TEXT 277 5 "noise" }{TEXT -1 107 " to our input data f or the dependent variable, we generate values of the same objective fu nction as before." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 176 "for j from 1 to 25 do \n y[j] :=evalf(eval(f,\{x=xx[j],seq(a[k]=p[k],k=1.. 4)\})) + ran()/2000.0 \+ od;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "With such imprecise data our former tolerance is inappropriate; we reset " }{TEXT 278 3 "tol" }{TEXT -1 61 " as follows and restate o ur initial estimates of parameters. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "tol := 10^(-6); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "p2 := vector(4, [3.2,4.7,7.3,1.9]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 82 "After these preparations we again invoke proced ure nlfit to repeat the fit of our " }{TEXT 279 5 "noisy" }{TEXT -1 6 " data." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "nlfit(f,x,xx,y,w ,a,p2,tol);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 304 "Examining these r esults we discover that, after only four iterations, final values of p arameters reproduce, within two estimated standard errors in most case s, the values used to generate input data. We plot both fitted points and curve of the fitted function to indicate visually the goodness of fit. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "with(plots):\npo i := pointplot(zip((xx,y) -> [xx,y],xx,y), symbol=circle, colour=blue) : " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "cur := plot(eval(f, \+ \{seq(a[j]=p2[j],j=1..4)\}), x=0..1.3):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "display(poi, cur);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 138 "To assess whether this fit suffers from systematic distortion \+ we plot residuals, which are differences between observed and fitted v alues." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 399 "df := vector(25) : \+ \nfor j to 25 do \n df[j] := y[j] - evalf(eval(f, \{x=xx[j], seq( a[k]=p2[k], k=1..4)\})) \n od: \+ \+ \+ " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "pointplot(zip((xx,df) -> [xx,df],xx,df), symbol=circl e, colour=orange);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "The residuals show no systematic t rend with increasing " }{TEXT 291 1 "x" }{TEXT -1 128 ", only a random scatter within a range [-0.25, 0.25] expected from the nature of appl ication of the generator of random numbers." }}}}{MARK "0 0 0" 56 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }