このアルゴリズムは非線形問題の近似の際、２次の項まで考慮しているため、比較的発散しにくいのが特徴です。

> restart;

> genfit:=module()
export mnlfit;
local gry,fitstats;
option package;

mnlfit := proc(f, X, x::array, y::vector, wt::vector,
P::{name,list(name)}, p::vector, tol)
beta,cvm,b,p1,errors,dp,res,nu,res1,sigma1,reso,nxvar,
yc,sw,sdf2w,sfo,sfo2,fstat;
# Initialize.
#print(f,X,x,y,wt,P,p,tol);
if type(P,name) then
#print(P);
#print([seq(P[k],k=1..linalg[vectdim](p))]);
RETURN(mnlfit(f,X,x,y,wt,[seq(P[k],k=1..linalg[vectdim](p))],p,tol));

fi;

reso := 1.0*10^9;
nobs := linalg[vectdim](y);
npar := linalg[vectdim](p);
nxvar := linalg[coldim](x);
xk := vector(nxvar);
r := vector(nobs);
a := matrix(nobs,npar,0.0);
corr := matrix(npar,npar,0.0);
drv := vector(npar);
p1 := vector(npar);
errors := vector(npar);
niter := 50;
nu := 1000;
# Begin iterations.
for n from 1 to niter do
#print(` iteration No. `,n);
# Compute initial fit and residuals ( = observed - calculated).
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
r[k] := evalf(sqrt(wt[k])*(y[k] - yck));
for j from 1 to npar do
a[k,j] := sqrt(wt[k])*drv[j];
od:
od:

# sum of squares of residuals
res1 := evalf(linalg[dotprod](r,r));
#print(` initial SSR = `, evalf(res1,5));
# Compute matrix alpha of normal equations.
alpha := evalf(evalm(linalg[transpose](a) &* a));
# Compute vector beta for right sides of normal equations.
beta := evalf(evalm(linalg[transpose](a) &* r));
# Note use of factor nu in this next loop; this increases as one approaches solution.
for j from 1 to npar do
alpha[j,j] := evalf(alpha[j,j]*(1. + 1./nu));
od:
# Solve a system of normal equations to correct parameters.
dp := evalf(linalg[linsolve](alpha, beta));
# Compute sum of squares of residuals at new point p1 = p + dp.
p1 := evalf(evalm(p + dp));
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
# xk := x[k];
r[k] := evalf(sqrt(wt[k])*(y[k]-yck));
od:
res := evalf(linalg[dotprod](r,r));

#print(`new SSR = `, evalf(res,5));
# Modify nu according to whether the sum of
# squares increases or decreases after this correction.
# For the wrong direction, elements are not corrected.
if res1 <= res then
nu := nu/10.;
else
# For the right direction, correct the elements.
nu := nu*10.;
p := evalf(evalm(p + dp));
fi:
# Test for convergence.
if abs(res - reso) <= tol then
break
fi;
reso := res;
od:
# covariance matrix
cvm := evalf(linalg[inverse](alpha));
# standard error of a single observation of unit weight
sigma1 := sqrt(res/(nobs - npar));
# errors of parameters
for j from 1 to npar do
errors[j] := evalf(sigma1*sqrt(cvm[j,j]));
od:
for 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:
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
od;
sdf2w := add((y[j] - yc[j])^2/wt[j], j=1..nobs):
fstat := (nobs-npar)/npar*(sfo2-sdf2w-sfo^2/sw)/sdf2w; # Print results.
# vector of corrections
#print(` vector of corrections `);
#print(seq(evalf(dp[j], 5), j=1..npar));
# final set of elements
#print(` final values of parameters `);
#print(` No. value error `);
#j := 'j': for j from 1 to npar do
# print(j, evalf(p[j], max(trunc(ln(abs(p[j]/errors[j]))/2.2+3), 4)),
# evalf(errors[j], 3))
# p(j):=evalf(p[j]);
#od;
RETURN(seq(p[j],j=1..npar));
end:
gry := proc(f,G,x,a,n,xval,p::vector)
# procedure to evaluate derivatives and formula
local k,y,i,drv,point;
point := {x=xval[1], seq(a[i]=p[i], i=1..n)};
drv := vector(n);
for k from 1 to n do
drv[k] := evalf(eval(G[k],point))
od:
y := evalf(eval(f,point));
RETURN(y,drv);
end:
fitstats:=proc(P)
local sigma1,fstat,npar,corr,nobs,xk,x,y,yc,p,wt,j,k,drv;
# standard deviation of one observation of unit weight
genfit:=module()
export mnlfit;
local gry,fitstats;
option package;

mnlfit := proc(f, X, x::array, y::vector, wt::vector,
P::{name,list(name)}, p::vector, tol)
beta,cvm,b,p1,errors,dp,res,nu,res1,sigma1,reso,nxvar,
yc,sw,sdf2w,sfo,sfo2,fstat;
# Initialize.
#print(f,X,x,y,wt,P,p,tol);
if type(P,name) then
#print(P);
#print([seq(P[k],k=1..linalg[vectdim](p))]);
RETURN(mnlfit(f,X,x,y,wt,[seq(P[k],k=1..linalg[vectdim](p))],p,tol));

fi;

reso := 1.0*10^9;
nobs := linalg[vectdim](y);
npar := linalg[vectdim](p);
nxvar := linalg[coldim](x);
xk := vector(nxvar);
r := vector(nobs);
a := matrix(nobs,npar,0.0);
corr := matrix(npar,npar,0.0);
drv := vector(npar);
p1 := vector(npar);
errors := vector(npar);
niter := 50;
nu := 1000;
# Begin iterations.
for n from 1 to niter do
#print(` iteration No. `,n);
# Compute initial fit and residuals ( = observed - calculated).
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
r[k] := evalf(sqrt(wt[k])*(y[k] - yck));
for j from 1 to npar do
a[k,j] := sqrt(wt[k])*drv[j];
od:
od:

# sum of squares of residuals
res1 := evalf(linalg[dotprod](r,r));
#print(` initial SSR = `, evalf(res1,5));
# Compute matrix alpha of normal equations.
alpha := evalf(evalm(linalg[transpose](a) &* a));
# Compute vector beta for right sides of normal equations.
beta := evalf(evalm(linalg[transpose](a) &* r));
# Note use of factor nu in this next loop; this increases as one approaches solution.
for j from 1 to npar do
alpha[j,j] := evalf(alpha[j,j]*(1. + 1./nu));
od:
# Solve a system of normal equations to correct parameters.
dp := evalf(linalg[linsolve](alpha, beta));
# Compute sum of squares of residuals at new point p1 = p + dp.
p1 := evalf(evalm(p + dp));
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
# xk := x[k];
r[k] := evalf(sqrt(wt[k])*(y[k]-yck));
od:
res := evalf(linalg[dotprod](r,r));

#print(`new SSR = `, evalf(res,5));
# Modify nu according to whether the sum of
# squares increases or decreases after this correction.
# For the wrong direction, elements are not corrected.
if res1 <= res then
nu := nu/10.;
else
# For the right direction, correct the elements.
nu := nu*10.;
p := evalf(evalm(p + dp));
fi:
# Test for convergence.
if abs(res - reso) <= tol then
break
fi;
reso := res;
od:
# covariance matrix
cvm := evalf(linalg[inverse](alpha));
# standard error of a single observation of unit weight
sigma1 := sqrt(res/(nobs - npar));
# errors of parameters
for j from 1 to npar do
errors[j] := evalf(sigma1*sqrt(cvm[j,j]));
od:
for 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:
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
od;
sdf2w := add((y[j] - yc[j])^2/wt[j], j=1..nobs):
fstat := (nobs-npar)/npar*(sfo2-sdf2w-sfo^2/sw)/sdf2w; # Print results.
# vector of corrections
#print(` vector of corrections `);
#print(seq(evalf(dp[j], 5), j=1..npar));
# final set of elements
#print(` final values of parameters `);
#print(` No. value error `);
#j := 'j': for j from 1 to npar do
# print(j, evalf(p[j], max(trunc(ln(abs(p[j]/errors[j]))/2.2+3), 4)),
# evalf(errors[j], 3))
# p(j):=evalf(p[j]);
#od;
RETURN(seq(p[j],j=1..npar));
end:
gry := proc(f,G,x,a,n,xval,p::vector)
# procedure to evaluate derivatives and formula
local k,y,i,drv,point;
point := {x=xval[1], seq(a[i]=p[i], i=1..n)};
drv := vector(n);
for k from 1 to n do
drv[k] := evalf(eval(G[k],point))
od:
y := evalf(eval(f,point));
RETURN(y,drv);
end:
fitstats:=proc(P)
local sigma1,fstat,npar,corr,nobs,xk,x,y,yc,p,wt,j,k,drv;
# standard deviation of one observation of unit weight
print(` standard deviation of fit = `, evalf(sigma1, 4));
print(` F statistic = `, evalf(fstat, 4));
print(` matrix of correlation coefficients `);
j := 'j': for j from 1 to npar do
print( seq(evalf(corr[j,k],5), k=1..j) ) od;
print(` table of residuals `);
print(` x obs `,` y obs `,` y calc `,
` difference `,` ratio `);
j := 'j': for j from 1 to nobs do
for k from 1 to nxvar do
xk[k] := x[j,k] od;
print(xk,y[j],yc[j],evalf(y[j]-yc[j], 5), evalf((y[j]-yc[j])*sqrt(wt[j]), 5) );
od:
end:
end module: print(` standard deviation of fit = `, evalf(sigma1, 4));
print(` F statistic = `, evalf(fstat, 4));
print(` matrix of correlation coefficients `);
j := 'j': for j from 1 to npar do
print( seq(evalf(corr[j,k],5), k=1..j) ) od;
print(` table of residuals `);
print(` x obs `,` y obs `,` y calc `,
` difference `,` ratio `);
j := 'j': for j from 1 to nobs do
for k from 1 to nxvar do
xk[k] := x[j,k] od;
```Warning, `genfit` is implicitly declared local to procedure `fitstats`