非線形最小二乗法プロシージャ

Mapleには非線形最小二乗法のプログラムが直接組み込まれていないため、ここでは Maple開発元(Waterloo Mapl Inc.,社) より公開されている
http://www.maplesoft.com/apps/categories/data_analysis_stats/stats/html/genfit_6.html
より、Levenberg-Marquardt 法のアルゴリズムを使用しています。

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

> 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)
local gradf,nobs,npar,k,yck,xk,r,a,drv,niter,j,n,corr,alpha,
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;

gradf := linalg[grad](f,P);
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;
yck,drv := gry(f,gradf,X,P,npar,xk,p);
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];
yck,drv := gry(f,gradf,X,P,npar,xk,p1);
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;
yc[k],drv := gry(f,gradf,X,P,npar,xk,p)
od;
sdf2w := add((y[j] - yc[j])^2/wt[j], j=1..nobs):
sfo := add(yc[j]*sqrt(wt[j]), j=1..nobs):
sfo2 := add(yc[j]*yc[j]*wt[j], j=1..nobs):
sw := add(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)
local gradf,nobs,npar,k,yck,xk,r,a,drv,niter,j,n,corr,alpha,
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;

gradf := linalg[grad](f,P);
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;
yck,drv := gry(f,gradf,X,P,npar,xk,p);
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];
yck,drv := gry(f,gradf,X,P,npar,xk,p1);
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;
yc[k],drv := gry(f,gradf,X,P,npar,xk,p)
od;
sdf2w := add((y[j] - yc[j])^2/wt[j], j=1..nobs):
sfo := add(yc[j]*sqrt(wt[j]), j=1..nobs):
sfo2 := add(yc[j]*yc[j]*wt[j], j=1..nobs):
sw := add(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;
yc[j],drv := gry(f,gradf,X,P,npar,xk,p);
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;
yc[j],drv := gry(f,gradf,X,P,npar,xk,p);
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:

Warning, `genfit` is implicitly declared local to procedure `fitstats`

>