--  General Non-Linear Least Squares Curve Fitting Program
--  by John C. Pezzullo, Ph.D.;  NIH Perinatology Research Branch, Georgetown Univ.
--  JohnP71@aol.com
-- 3/26/01 modified to use solve() routine in SVD Xfun

--  Define the function to be fitted
f(x,p) = p[3] + ( p[1] - p[3] ) / ( 2 ^ ( x / p[2] ) )  --  Exponential Decay

Par = {110,4.,7}  --  Provide guesses for parameters here

--  Provide observed values and precisions here
x  = {  0,   1,   2,  3,  4,  5,  6,  7,  8,  9}
yo = {112, 108, 102, 98, 94, 89, 87, 84, 83, 82}
yErr := 1 + 0.0 * yo:
--  (for unweighted LSQ fit, use  yErr := 1 + 0 * yo)

--  Computer figures these out
nPar:3.00000000;  nPts:10.00000000;  DegFrdm:7.00000000
w := 1 / yErr^2:--  Weight factors are inversely proportional to square of std errs of y's

--  Iterate to find Least Squares Parameters
dPar:={0.1,0.1,0.1}, Cycles:=0,
     (Cycles:=Cycles+1, yc:=ycf, dy:=yo-yc, D:=Df, A:=Af, R:=Rf, solve(A,R),,
      dPar:=_X, Par:=Par+dPar) while sum(abs(dPar[i]),i,1,nPar)>1e-12:

Cycles:15.00000000--  # of iterations it took to converge

--  Fitted Parameters
table Par
 112.94809928
 4.65076188
 69.70573693

--  Agreement between observed and calculated values
Nmin=1;  Nmax=nPts;  Nsteps=nPts
--table  x[N],     yo[N],     yc[N],     dy[N]

plot {x,yo};  plot f(X,Par) --  Graph the results

--  Function definitions
nPts = count(x)                   --  # of data points
nPar = count(Par)                 --  # of adjustable parameters
DegFrdm= nPts-nPar                --  degrees of freedom
ycf = f(x,Par)                    --  Calculated dependent variable

pm(p,k)[i] = 1.01*p[i] when i=k, p[i] otherwise   --  increment the kth param by 1%
d(x,p)[k] = 100*(f(x,pm(p,k))-f(x,p))/p[k]        --  Deriv w/resp to kth param
Df[i,k] = d(x[i],Par)[k]  dim[nPts,nPar]          --  Derivs of each pt w/resp to each param
Af[j,k] = sum(w[i]*D[i,j]*D[i,k],i,1,nPts) dim[nPar,nPar]  --  Normal Equations Matrix
Rf[k] = sum(w[i]*D[i,k]*dy[i],i,1,nPts)  dim[nPar]         --  Norm Equations Right Hand Side