# model must be of the form Y = F(t[k],P), where P are parameters
# and t is the independent variable. The sumindex subscript [k] MUST
# be attached to t in the model equation. For example,
# leastsqrs( y=a+b*x[k], x, k, [a,b] );
# Also valid is
# leastsqrs( y[k]=a+b*x[k], x, k, [a,b] );
#---------------------------------------------------------------------
leastsqrs := proc( model::`=`, indepvar::name,
sumindex::name, params::list )
local chi2, y, i, t, partials, k, eqs, solset, T, Y,
replace_denominators, sums, YTsums, n, p, nmax, N,
Tcount, Scount, ispoly, Delta, delta, tmp, goodsols;
global _subslist, time0;
if not has(model,sumindex) then
ERROR("model must be of the form Y=F(t[i],P)");
fi;
time0 := time();
i := sumindex;
t := indepvar;
if not type(t,indexed) then
t := t[i];
fi;
y := lhs(model);
if not type(y,indexed) then
y := y[i];
fi;
chi2 := Sum((y-rhs(model))^2,i=1..N);
if type(rhs(model),polynom) then
ispoly := true;
else
ispoly := false;
fi;
debug_print(procname,"i,t,y,ispoly,chi2",4,i,t,y,ispoly,chi2);
#-----------------------
# calculate the partials
#-----------------------
debug_print(procname,"Forming the normal equations...",3);
partials := [];
for k from 1 to nops(params) do
partials := [ op(partials), diff(chi2,params[k])=0 ];
od;
if printlevel >= 4 then
debug_print(procname,"partials",4);
for p in partials do
print(p);
od;
fi;
#-------------------------------------------
# form the summation substitutions variables
#-------------------------------------------
sums := convert( select( (x,y)->op(0,x)=y,
Page 5
