実験系では非線形の関数を得られた実験データにフィッティングさせることはよくある。簡単な関数だったExcelでできるし、複雑でも有名な関数だったら、専門のソフトを使えばできる。でも、複雑な関数だったり、多量のデータを連続的に、色々な関数に対してフィッティングさせようとしたら、自分でプログラムを組む必要が出てくる。ここではフリーソフトであるScilabやOctaveを使って、タダかつ正確かつ高速に非線形フィッティングすることを目指す。

1. Octaveの場合

leasqrを使うとできる。この関数は、そもそも回帰計算用に作られているので、わかりやすい。

Levenberg-Marquardt algorithm（レーベンバーグ・マルカート法）という方法を使っているらしい。ちなみに今のとこ一番いい非線形回帰法らしい。詳しいことは全く分かりません(^^;

[f, p] = leasqr(x, y, pin, F)

x, yデータに対して関数Fをフィッティングし、最適なパラメータを返す、という関数。超便利♪

x, yは測定データベクトルで同じ行数、1列のベクトル（別に1列である必要はないのかもしれない。。。）

pinは最適化したいパラメータの初期値を入れたベクトル。

Fはy=f(x,p)という形の関数で、フィッティングさせる式。y1, y2...ymだとm行1列の関数ベクトルになる。

で、何が返ってくるかというと…

まずpは最適化されたパラメータベクトル、pinからスタートしてpに至るわけだ。

fは最適化されたときの関数の値のベクトル。Fと同じくm行1列。

東京大学の関先生が、この関数を使って、"SWRC Fit" という土壌水分特性曲線の回帰プログラムを書かれています。

2. Scilabの場合

lsqrsolve関数を使うと、leasqrと同様にLevenberg-Marquardt algorithm（レーベンバーグ・マルカート法）を使った非線形回帰が行える。ただし、回帰用だけに作られた関数ではないので、フィッティングに利用するには多少の工夫が必要。

※leastsq関数やdatafit関数を使えば、Octaveのleasqrと同じように、簡単に記述できることに気づく（Oct. 29, 2007）。ただし非線形回帰手法にLevenberg-Marquardtを使いたい場合は、lsqrsolveを使うしかないようだ。

[x, v] = lsqrsolve(x0, fct, m)

非線形関数fctの2乗和（関数は一つでも複数でもよい）を最小にするxを求める関数。汎用性高し☆

x0は最適化したいパラメータの初期値を入れたベクトル、leasqrのpinに相当している。

fctはv=fct(x, m)という形で、v1, v2...vmだとm行1列の関数ベクトルになる。mは関数の個数を表すことになる。leasqrと違って、フィッティングさせる場合にはfctは誤差の関数になるわけだ。

mは前述したように、関数の個数。

で、何が返ってくるかというと…

まずxは最適化されたパラメータベクトル、xoからスタートしてxに至るわけだ。

vは最適化されたときの関数の値のベクトル。Fと同じくm行1列。

データフィットプログラムの例、2次関数の例

//まずフィッティングさせるデータを作る。
a=34; b=12; c=14; //2次関数の係数を定義

function y = FF(x)
y = a * (x-b)+ c * x .* x;
endfunction //2次関数を定義

X=(0:.1:3)';
Y=FF(X)+100*(rand()-.5); //乱数を使って誤差を加えてやる。データは0から3まで0.1毎に定義

//誤差を加えたデータに対して2次関数をフィッティングさせる
function e = f1(abc, m)
a=abc(1); b=abc(2); c=abc(3);
e = Y - (a * (X-b) + c * X .* X);
endfunction //最小化させる誤差関数を定義

[abc, v] = lsqrsolve([10; 10; 10], f1, size(X,1)); //abcは最適化されたパラメータ, vはそれを
//代入したときの関数の値、[10; 10; 10]はパラメータ初期値、f1は2乗和を最小化する関数、size(X,1)で
//Xの行数つまり、関数の個数が求まる。ちなみにsize(X,2)だと列数を求める。
abc
a=abc(1); b=abc(2); c=abc(3);
y=FF(X);
plot(X, Y, '+', X, y, '-')