二次元空間上に分布する2つ以上の座標集合を直線で近似する際、最小二乗法という手法を使います。
定義
直線の式を一次関数
[math]f(x)=ax+b\,[/math]
とおくと、直線の傾き a と切片 b は次の式で求めることができます。
[math]a=\frac{\displaystyle n\sum_{k=1}^n x_ky_k-\sum_{k=1}^n x_k\sum_{k=1}^n y_k}{\displaystyle n\sum_{k=1}^n x^2_k-\left( \sum_{k=1}^n x_k \right)^2}[/math]
[math]b=\frac{\displaystyle \sum_{k=1}^n x^2_k\sum_{k=1}^n y_k-\sum_{k=1}^n x_ky_k\sum_{k=1}^n x_k}{\displaystyle n\sum_{k=1}^n x^2_k-\left( \sum_{k=1}^n x_k \right)^2}[/math]
Ruby スクリプト
この式を Ruby で書いていきたいと思います。
まず、座標値の集合を配列で定義します。
s = [[ 32.0, 24.0],
[131.0, 62.0],
[159.0, 110.0],
[245.0, 146.0]]
今回は4つの座標値をサンプルとして用います。配列の中は、
[[x0, y0], [x1, y1], ... ]
という構造として扱います。 次に、a と b を求めなければいけないのですが、上の式で頻出している数列の総和[math]\sum_{k=1}^n[/math]をシングルトンメソッドとして座標値集合 s に突っ込みます。
def s.sum(&b)
self.inject(0.0) do |t, pnt|
t + b.call(pnt[0], pnt[1])
end
end
この sum メソッドは、自身の全ての要素について与えられたブロック引数を処理し、その結果の総和を返すメソッドです。Enumerable モジュールに組み込まれている inject メソッドもそもそも総和を求めるメソッドなのでそのままaとbを求める式に利用することも可能ですが、呼び出し側で最低限の記述をしたいために sum メソッドでラッピングしています。
つまり、
[math]\sum_{k=1}^n x_k+y_k[/math]
という式は、
s.sum {|x, y| x + y}
というシンプルなブロック構文付きのメソッド呼び出しに置き替えることができます。今回の場合、k は 1 から総数 n まで、つまり全ての要素についての処理のみしかありませんので、sum メソッドにはこれらを指定する引数は取りません。
この sum メソッドを用いて n, a, b をそれぞれ記述すると、次のようになります。
n = s.size.to_f
a = (n * s.sum{|x,y|x*y} - s.sum{|x,y|x} * s.sum{|x,y|y}) /
(n * s.sum{|x,y|x**2} - (s.sum{|x,y|x})**2)
b = (s.sum{|x,y|x**2} * s.sum{|x,y|y} - s.sum{|x,y|x*y} * s.sum{|x,y|x}) /
(n * s.sum{|x,y|x**2} - (s.sum{|x,y|x})**2)
これを実行すると
p a => 0.5913598269802649
p b => 1.6747445255474454
という結果になりました。これだけ見ると正直、合っているのかどうか分かりません。
グラフにプロット
そんなわけで gnuplot でプロットしてみます。 まずは座標値群ファイルを準備します。
$ cat ./ten
32.0 24.0
131.0 62.0
159.0 110.0
245.0 146.0
次に先ほど求めた直線の式から、[math]x = 0[/math] の時と [math]x = 300.0[/math] の時の x, y それぞれの値を求めます。
if b.nan?
# 切片の解なし(Y軸に対して平行)
else
x = 0.0
puts "#{x} #{a*x+b}"
x = 3000.0
puts "#{x} #{a*x+b}"
end
結果は
0.0 1.6747445255474454
300.0 179.08269261962693
となります。 この結果を ./sen に保存しておいて、
$ gnuplot4-qt
gnuplot> plot "./ten", "./sen" w lp
すると次のように表示されました。
んー、なかなか合ってるような気がします！
どうもちゃんと動いているようなので、以下に Ruby スクリプトの全文を掲載しておきます。