最小２乗法

２次元座標データの直線近似

　 n 個の２次元座標点（ x i , y i ）を直線 ( y - y 0 ) = a ( x - x 0 ) により 近似する．ベクトル X = ( x 1 , x 2 , ..., x n)T ，および， Y = ( y 1 , y 2 , ..., y n)T により， ( Y - y 0 ) ≒ a ( X - x 0 ) とする． さらに， Y' = Y - y 0 , X' = X - x 0 とすると， Y' ≒ a X' により， a = ( Y' X'T ) / ( X' X'T ) が求められる．
　ここで，始点 （ x 0 , y 0 ） を適当に決めると，そのときの２乗誤差を決定数することができる． 始点の決定方法にはいろいろな方法が考えられるが，簡単なのは，与えられた２次元座標点（ x i , y i ）の １つを用いるものである．それぞれの座標点を始点として２乗誤差を求め，それが最小となるものを計算結果とする．


プログラムリスト

上記の文字に対応していない場合があります．コメントを参照してください．


//最小２乗法処理ルーチン
#include 
//  n : 入力はデータの総数，出力は始点のインデックス
double CFractalDim::Minimize2error(double *x, double *y, int& n)
{
	int loop, center, P;
	double Y, Z, DIV, error, p;

	double mse = DBL_MAX, DIVMAX = 0.0;
	for(center = 0; center < n; center++){						// 始点の移動
		Y = Z = 0.0;
		for(loop = 0; loop < n; loop++){
			p = x[loop] - x[center];
			Y += p * p, Z += p * (y[loop] - y[center]);
		}
		DIV = Z / Y;									//傾きを決定する

		//近似直線の式    y = DIV * (x - x[center]) + y[center];

		//誤差を計算
		error = 0.0;
		for(loop = 0; loop < n; loop++){
			p = DIV * (x[loop] - x[center]) - (y[loop] - y[center]);
			error += p * p;
		}

		//最小かどうか調べる
		if(error < mse){
			mse = error, DIVMAX = DIV, P = center;
		}
	}
	n = P;

	return DIVMAX;										// 傾きを戻り値とする
}

言語は一応C++です．（C言語も混ざっていますが）．
バグがありましたら連絡してください．お願いします．
戻る
2002/4
by Hiroyuki Koga