対数の最小二乗

◆対数へのフィッティング

　何かの現象を捕らえたデータがあって、その傾向を知りたいときや予測がしたい場合に関数へのフィッティングを行いたくなります。Excelの機能を使うと簡単に求めることができます。いったいどのように計算しているのでしょうか？今回は対数へのフィッティングを行ってみたいと思います。

　以下のJavaアプレットは、マウスで5点クリックすると対数のフィッティング結果を表示するサンプルです。

◆対数にフィッティングさせる簡単な方法

　このプログラムはマウスでクリックされた座標を、y = a ln(x) + b の式に最小二乗法を使ってフィッティングさせています。(lnは自然対数の意味です)

　最小二乗法を使った簡単な例は、y = ax + b の１次関数へのフィッティングです。目的の関数と入力データの二乗誤差が最小になるように a と b の値を求めます。式で書くと次のようになります。

J = Σ( y - (ax + b))^2 → min

　この問題を解くために、Jの式で a と b の編導関数がゼロになることを利用します。途中の計算を飛ばして結論を書いてしまうと、以下の連立方程式を解けばよいことになります。

　ここで、y = ax + b へのフィッティング方法が示されました。では、y = a ln(x) + b への拡張を行いたいと思います。といっても x を ln(x) に置き換えるだけです。

　機械的に置き換えれば終了です。

◆プログラミングはどうするの

　以下の手順が主な内容です。

　手順１．マウスでクリックされた座標の配列を入力データとして、先に導いた連立方程式を作ります。

　手順２．連立方程式を解いて a と b を求めます。

　手順３．a と b を使ってフィッティング結果を描画します。

◆手順１

　連立方程式を2次元配列にして、値を足しこむだけでOK簡単ですね。

for(int i=0;i<_x.length;i++)
{
	double lx = Math.log(_x[i]);
	mtx[0][0] += lx*lx;
	mtx[0][1] += lx;
	mtx[0][2] += lx*_y[i];
	mtx[1][2] += _y[i];
	mtx[1][1] += 1;
}
mtx[1][0]=mtx[0][1];

　こんな感じに配列に値を入れるだけです。

◆手順２

　連立方程式を解くわけですが、簡単なので直接 a b を求めることもできます。しかし、今後の拡張性も考えて汎用的なロジックを実装します。ガウスジョルダン法で連立方程式を解く関数を作成しました。ガウスジョルダン法の詳しい説明は沢山ありますので割愛します。

if( !GaussJordan(mtx,2) )
	return null;

　作成した関数を呼び出せば終わりです。

◆手順３

　アプレットのpaint(Graphics g)メソッドで描画します。

for(double x=minx;x<maxx;x+=2)
{
	double y = ab[0] * Math.log((double)x) + ab[1];
	g.drawRect((int)x, (int)y, 1, 1);
}

　ポリラインで書けば綺麗と思いますが、ちょっと手抜きをして小さな四角を並べました。

◆ソースのダウンロードはこちら
Latest update at 2007/4/24