入力誤差モデルを用いた最尤推定による直線当てはめ

gistfile1.js

// 入力誤差モデルを用いた最尤推定による直線当てはめ
// 入力 データ点 points = [[x0,y0],[x1,y1],...,[xn-1,yn-1]]
// y = ax + b の a, b を格納した配列 [a,b] を返す
function lineFittingInputErrorModel(points) {

	// 点の数
	var N=points.length;	

	// xの重心
	var mx=0;	
	for(var i=0; i<N; ++i) {
		mx+=points[i][0];
	}
	mx/=N;

	// yの重心
	var my=0;	
	for(var i=0; i<N; ++i) {
		my+=points[i][1];
	}
	my/=N;

	var M=numeric.rep([2, 2], 0);
	for(var i=0; i<N; ++i) {
		M[0][0]+=(points[i][0]-mx)*(points[i][0]-mx);
		M[0][1]+=(points[i][0]-mx)*(points[i][1]-my);
		M[1][1]+=(points[i][1]-my)*(points[i][1]-my);
	}
	M[1][0]=M[0][1];

	// get minimum eigen value
	var lambda=0.5*(M[0][0]+M[1][1]-Math.sqrt((M[0][0]+M[1][1])*(M[0][0]+M[1][1])-4*(M[0][0]*M[1][1]-M[0][1]*M[1][0])));
	var vec=[lambda-M[1][1], M[1][0]];
	vec=numeric.div(vec, numeric.norm2(vec));

	// 直線 Ax + By + C = 0 の係数
	var A=vec[0];
	var B=vec[1];
	var C=-A*mx-B*my;

	// y = ax + b の a, b を返す
	// （注意）B=0の場合を考慮していない
	return [-A/B,-C/B];
}