最小二乗法

一個前のガウスの消去法をつかって最小二乗法も求められるようにしよう・・。

実験などで例えば理論値f=x+6に対して以下の実験結果が得られたとする

x , f
-4 , 2.15
-3 , 3.32
-2 , 3.78
-1 , 4.99
0 , 6.23
1 , 6.89
2 , 7.12
3 , 9.18
4 , 9.73

これの数値部分がdata.txtに保存されているとして、これを読み込んで実行すれば

lsm = LeastSquaresMethod.new
lsm.input_data("data.txt",",")
p lsm.polynomial(1) #=>[0.941333333333333, 5.93222222222222]

と、近似直線はf=0.94x+5.93と求まり、確かに近似直線として良さそうである。

ここで、input_dataの第二引数markはデータの区切り文字列で、
空白とかタブとかでxとfを区切ってあるファイルでも大丈夫なようにしてある。
ploynomialは読み込んだデータ点から近似多項式を求めるメソッドで、
引数mは多項式の次数である。
引数mを2にすれば二次関数の多項式の近似を出力する。

ソースファイル least_squared_method.zip

require 'gaussian_elimination'
class LeastSquaresMethod
  def initialize
    @data = []
  end
  
  def input_data(filename, mark = ",")
    @data = []
    open(filename,"r") do |file|
      file.each_line do |line|
        @data.push(line.split(mark).map{|v| v.to_f})
      end
    end
  end
  
  def polynomial(m = 1)
    return "Dimension Error" if m < 1
    return "No Data" if @data.empty?
    invs_x,invs_f = involutions(m)
    equations = Array.new(m + 1){[]}
    equations.each_with_index do |eqation, n|
      (m+2).times do |k|
        if k == m + 1
          p = m - n
          eqation.push(invs_f[p])
        else
          p = 2*m - n - k
          eqation.push(invs_x[p])
        end
      end
    end
    return GaussianElimination.new.solve(*equations)
  end
  
  def involutions(m)
    data_x, data_f = @data.transpose
    invs_x = [data_x.size]
    invs_f = [data_f.inject(:+)]
    1.upto(2*m) do |n|
      invs_x.push(data_x.inject(0){|sum,x| sum + x**n})
      invs_f.push(@data.inject(0){|sum,d| sum + d[1]*(d[0]**n)}) if n <= m
    end
    return invs_x,invs_f
  end
end

lsm = LeastSquaresMethod.new
lsm.input_data("data.txt",",")
p lsm.polynomial(1)