Skip to content

Instantly share code, notes, and snippets.

@sasekazu
Created January 12, 2015 07:31
Show Gist options
  • Save sasekazu/f158372b03a43d5f7011 to your computer and use it in GitHub Desktop.
Save sasekazu/f158372b03a43d5f7011 to your computer and use it in GitHub Desktop.
入力誤差モデルを用いた最尤推定による直線当てはめ
// 入力誤差モデルを用いた最尤推定による直線当てはめ
// 入力 データ点 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];
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment