最小二乗近似曲線の式を求める - コンピュータを楽しもう！！

ポケコンを目指しているLuaridaですが、Luaridaを作って、先ずやってみたかったのが、実はこれでした。最小二乗近似やりたいよね。ということで、作ってみました。
プログラムはここからダウンロードしてください。(LeastSquare)

グローバル宣言

関数は下記のようなものを用意しました。元データはcsvファイルに用意します。それをValue.XとValue.Yの配列に読み込みます。

------------------------------------------
--最小二乗近似曲線を求めるプログラム
------------------------------------------
--関数宣言--------------------------------
main={}      --mainメソッド
printscrl={} --スクロール有り文字表示
gauss={}     --ガウス消去を行う
split={}     --文字の分解
readdata={}  --CSVファイルを読み込みます
kinji={}     --近似計算

--グローバル変数宣言----------------------
Path = "/sdcard/luarida/"  --luaファイルを保存しているPath
Fname = "leastsquare.csv"  --データCSVファイル
Value ={ X={}, Y={} }      --Value.X[]とValue.Y[]の配列にデータの値を読み込む

元データ(leastsquare.csv)

csvファイルの読み込み

filenameで指定されたファイルからcsvデータを読み込みます。読み込んだデータは画面に表示されます。

------------------------------------------
--CSVファイルを読み込みます
------------------------------------------
function readdata( filename )
local fp
local msg
local str
local t={}
local i = 1
 --ファイルを開きます
 fp,msg = io.open( Path..filename, "r")
 if( not(fp) )then
   dialog( Path..Fname.."がオープンできません","プログラムを終了します", 0 )
   return
 end
 --CSVデータを読み込みます
 while(true)do
   str = fp:read("*l")  --1行読み込み
   if( str==nil )then break end  --読込むデータが無ければ終了
   str = string.gsub( str,"\r","" ) --改行コードを外す
   t = split( str, "," )
   if( t[1]~=nil and t[2]~=nil )then
     Value.X[i] = t[1]
     Value.Y[i] = t[2]
     printscrl( "X="..Value.X[i].." Y="..Value.Y[i], 20, color(0,0,0), color(255,255,255) )
     i = i + 1
   end
 end
 io.close(fp)
end

近似式を求めるために一次方程式を作る

読み込んだデータから近似式の係数を求める連立一次方程式を作ります。まず、何次の回帰曲線を求めるかを選びます。プログラムでは1〜6次式まで求めることができるようにしました。

ラジオボタンで指定します。

 item.clear()
 item.add( "1次近似", 0 )
 item.add( "2次近似", 0 )
 item.add( "3次近似", 0 )
 item.add( "4次近似", 0 )
 item.add( "5次近似", 0 )
 item.add( "6次近似", 0 )
 k = item.radio("何次近似を行いますか", 1 )

3次の回帰曲線を求めてみます。
それでは方程式を立てます。詳しい説明は省きますが、昔に書いた資料がこちらのサイトに残っていますのでこちらを見てください。

------------------------------------------
--近似計算
-- k次近似
------------------------------------------
function kinji( k )
local gn
local a={
        {0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        ,{0,0,0,0,0,0,0,0}
        }
local n = #Value.X
local px
local i,j
 gn = k + 1
 a[1][1] = n
 for i=1, n do
   px = 1
   for j=1,gn-2 do
     px = px*Value.X[i]
     a[1][j+1] = a[1][j+1] + px
   end
   px = 1
   for j=3,gn do
     px = px*Value.X[i]
   end
   for j=1,gn do
     px = px*Value.X[i]
     a[j][gn] = a[j][gn] + px
   end
   px = 1
   for j=0,gn-1 do
     a[j+1][gn+1] = a[j+1][gn+1] + Value.Y[i]*px
     px = px * Value.X[i]
   end
 end
 for i=2,gn do
   for j=1,gn-1 do
     a[i][j] = a[i-1][j+1]
   end
 end
 --for i=1,gn do
 --  printscrl( a[i][1].." "..a[i][2].." "..a[i][3].." "..a[i][4].." "..a[i][5], 20, color(0,0,0), color(255,255,255) )
 --end
 return( gauss( a, gn ) )
end

これで、aに求める連立方程式が設定できました。

ガウス消去

連立方程式を解くためにガウス消去法を使います。aに対して行列を入れ替えていきます。下記のプログラムはみその計算物理学のサイトに説明しているガウスの消去法の数値計算例(C言語)をLua言語に書き直したものです。

------------------------------------------
--ガウス消去を行います
-- a[n][n+1]配列です
------------------------------------------
function gauss( a, n )
local i,j,k,l
local m
local pivot
local b
local p,q
local x={}
 for i=1,n do
   m = 0
   pivot = i
   for l=i,n do
     if( math.abs( a[l][i] )>m )then
       m = math.abs( a[l][i] )
       pivot = l
     end
   end
   if( pivot~=i )then
     for j=1,n+1 do
       b = a[i][j]
       a[i][j] = a[pivot][j]
       a[pivot][j] = b
     end
   end
 end
 for k=1,n do
   p = a[k][k]
   a[k][k] = 1
   for j=k+1,n+1 do
     a[k][j] = a[k][j] / p
   end
   for i=k+1,n do
     q = a[i][k]
     for j=k+1,n+1 do
       a[i][j] = a[i][j] - q*a[k][j]
     end
     a[i][k] = 0
   end
 end
 for k=1,n do
   i = n - k + 1
   x[i]=a[i][n+1]
   for l=i+1,n do
     j = n - l + i +1
     x[i] = x[i] - a[i][j]*x[j]
   end
 end
 return( x )
end

結果表示

最後に求めた係数を表示します。
(追記:係数の表示が逆です。間違っています。最小二乗近似曲線の式を求める(2)に追記しています。)リンクはこちらです。

 for i=1,k+1 do
   printscrl( "X^ "..(i-1).."の係数: "..x[i], 20, color(0,0,0), color(255,255,255) )
 end

以上です。
グラフを描きたい場合は、グラフを描くのプログラムを持ってくれば、グラフが描けると思います。