Hatena::ブログ(Diary)

mmitouの日記

2012-09-20

多項式曲線フィッティングを実装してみた

| 18:43

PRMLの"1.1 例:多項式曲線フィッティング"を読んだ。

実際に実装して、科書通りにフィッティング出来ることを確認した。

多項式曲線フィッティングの実装

xを0~1の間で10個選び、それぞれのxに対してsin(2πx)を計算した値にノイズを加えて訓練集合tを作成する。

そして、この訓練集合tを以下の式でフィッティングする。

y¥left( x,w ¥right) = w_0 + w_1 x + w_2 x^2 + ¥ldots + w_M x^M = ¥sum_{j=0}^{M} w_j x^j ¥hspace{50}  (1.1)

フィッティングは、訓練集合tに対して二乗和誤差が最小になるwを見つける事で行う。

二乗和誤差は以下の式で表される。

 E(w) = ¥frac{1}{2} ¥sum_{n=1}^N ¥{ y(x_n,w)-t_n¥}^2 ¥hspace{50} (1.2)

(1.2)に(1.1)を代入して、各wで偏微分した結果が0になる時に二乗和誤差が最小である。

 ¥frac{¥partial E(w)}{¥partial w_i}=¥sum_{n=1}^{N}¥{¥sum_{j=0}^{M}w_j x_{n}^{j} -t_n¥}x_n^i ¥¥ = ¥sum_{n=1}^{N}¥{¥sum_{j=0}^{M}w_j x_{n}^{i+j} -t_n x_n^i¥} ¥¥ = ¥sum_{n=1}^{N} ¥sum_{j=0}^{M} w_j x_n^{i+j} - ¥sum_{n=1}^{N} t_n x_n^i ¥¥ = ¥sum_{j=0}^{M}w_j ¥sum_{n=1}^{N}x_n^{i+j} - ¥sum_{n=1}^{N}t_n x_n^i ¥¥ = ¥sum_{j=0}^{M} A_{ij} w_j - T_i = 0

したがって、

 ¥sum_{j=0}^{M} A_{ij} w_j = T_i

 A_{ij} = ¥sum_{n=1}^{N} x_n^{i+j} ¥¥ T_i = ¥sum_{n=1}^{N} t_n x_n^i

のように書く事が出来て(演習1.1)、M+1個の一次方程式から連立一次方程式が得られる。

その連立一次方程式はAw = Tの行列形式にして解く事が出来る。

ただし、このままでは過学習が発生するので、以下のように罰金項を加える。

 ¥widetilde{E(w)} = ¥frac{1}{2} ¥sum_{n=1}^{N} ¥{y(x_n, w) - t_n¥}^2 + ¥frac{¥lambda}{2}¥| w¥| ^2 ¥¥ ¥| w ¥| = w_0^2 + w_1^2 + w_2^2 + ¥ldots + w_M^2  ¥hspace{50} (1.4)

(1.4)も、偏微分して0になる条件から最適なwを求める。

(1.4)に(1.1)を代入して各wで偏微分すると、

 ¥widetilde{E(w)} = ¥frac{1}{2} ¥sum_{n=1}^{N} ¥{ ¥sum_{j=0}^{M} w_j x^j -t_n¥}^2 + ¥frac{¥lambda}{2} ¥sum_{m=0}^{M} w_m^2

 ¥frac{¥partial ¥widetilde{E(w)}}{¥partial w_i} = ¥sum_{j=0}^M A_{ij} w_j - T_i + ¥lambda w_i = 0

となり、

 A_{ij}’ = ¥begin{cases} ¥sum_{n=1}^{N} x_n^{i+j} & (i ¥neq j) ¥¥ ¥lambda + ¥sum_{n=1}^{N} x_n^(i+i) & (i = j) ¥end{cases}

とおくと、

 ¥sum_{j=0}^{M} A_{ij}’ w_j = T_i

が得られ(演習1.2)、ここから連立一次方程式を解くことでwを得る。

この一連の計算を実装したコードは以下の通りである。

import System.Random
import Control.Monad
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Data.List

trainingDataFileName = "training.dat"

-- 多項式の係数を計算 mは次数
polynomialCoefficients :: Field a => Bool -> a -> Int -> [a] -> [a] -> [a]
polynomialCoefficients isRegularized lambda m xn tn = [mtrxW @@> (i, 0) | i <- [0 .. m]]
  where
    m'    = m + 1
    -- Ai,j と Ti は 演習1.1を参照せよ
    a i j = sum $ map (^ (i + j)) xn
    -- 式1.4のように正則化した場合、Ai,jは i == jで係数λが加えられる
    -- 正則化した場合のAi,jの導出は演習1.2を参照せよ
    a' i j | i /= j    =  sum $ map (^ (i + j)) xn
           | otherwise =  lambda + (sum $ map (^ (i + j)) xn)
    t i   = sum $ zipWith (*) tn $ map (^ i) xn
    -- Aw = T の行列として扱い、連立方程式を解く
    fA | isRegularized = a'
       | otherwise     = a
    mtrxA = (m' >< m') [fA i j | i <- [0 .. m], j <- [0 .. m]]
    mtrxT = (m' >< 1 ) [t i| i <- [0 .. m]]
    mtrxW = linearSolve mtrxA mtrxT

-- ノイズを付与する
addNoise :: (RandomGen g, Random a, Num a) => (a, a) -> [a] -> g -> [a]
addNoise range xs g = zipWith (+) xs (randomRs range g)

-- gnuplotでフィッティング曲線, 訓練データ, sin(x)をpng形式で出力する
toFittingCurveScript :: Show a => [[a]] -> String -> String
toFittingCurveScript wss attr = "set xrange [-0.5: 2*pi + 0.5]\n"
                           ++ "set yrange [-1.5: 1.5]\n"
                           ++ "set term png\n"
                           ++ (concat $ map order (zip [0..] wss))
  where
    order (m, ws) = "set output \"order" ++ attr ++ show m ++ ".png\"\n"
                    ++ "set title " ++ "\"order " ++ attr ++ show m ++ "\"\n"
                    ++ "plot " ++ (toPolynomial ws)
                           ++ " w l lt 1 t \"fitting\", " 
                           ++ "sin(x) lt 2, "
                           ++ "\"" ++ trainingDataFileName ++ "\" w p pt 6 lt 3\n"
    toPolynomial xs = concat . map toTerm $ zip [0..] xs
    toTerm (i, coefficient) = "+ (" ++ show coefficient ++ ")" 
                              ++ (concat . take i $ repeat "*x")

erms :: Floating a => [a] -> [a] -> [a] -> a
erms xs ts ws = (2 * e / n) ** 0.5
  where
    f x = (sum . zipWith (*) ws $ scanl (*) 1 $ repeat x)
    e   = sum . map (^2) . zipWith (-) ts $ (map f xs)
    n   = genericLength xs

ermsList :: Floating a => [a] -> [a] -> [[a]] -> [a]
ermsList xs ts wss = map (erms xs ts) wss

valuesToString :: Show a => [a] -> String
valuesToString xs = concat $ map (\x -> show x ++ "\n") xs

main = do
  gen <- getStdGen
  let xs   = linearScale 10 (0,2*pi)
      ts   = addNoise (-0.5 :: Double, 0.5 :: Double) (fmap sin xs) gen
      wss  = map (\m -> polynomialCoefficients False 0   m xs ts) [0 .. 10]
      wss' = map (\m -> polynomialCoefficients True  1.5 m xs ts) [0 .. 10]
      testxs = linearScale 100 (0,2*pi)
      testys = addNoise (-0.5 :: Double, 0.5 :: Double) (fmap sin testxs) gen

  writeFile "plot.plt" $ toFittingCurveScript wss ""
  writeFile "plot_reg.plt" $ toFittingCurveScript wss' "_regularized"
  writeFile "erms_train.dat" . valuesToString $ ermsList xs ts wss
  writeFile "erms_test.dat"  . valuesToString $ ermsList testxs testys wss
  writeFile trainingDataFileName $ concat . map 
                (\(x,y) -> show x ++ " " ++ show y ++ "\n") $ zip xs ts

実行手順は以下の通り。

ghc test.hs
./test
gnuplot -p < plot.plt
gnuplot -p < plot_reg.plt

正則化無しで過学習する(1.1)式は以下のようになる。

f:id:mmitou:20120920182555p:image

f:id:mmitou:20120920182556p:image

f:id:mmitou:20120920182557p:image

f:id:mmitou:20120920182558p:image

正則化ありで、過学習を抑制した結果は以下の通り。

f:id:mmitou:20120920182559p:image

f:id:mmitou:20120920182600p:image

f:id:mmitou:20120920182602p:image

f:id:mmitou:20120920182601p:image

f:id:mmitou:20120920183610p:image

Ermsを計算した結果は、erms_train.datとerms_test.datに書きだすので、

gnuplotで以下のように実行して結果を表示する。

set term png
set output "erms.png"
plot "erms_train.dat" w lp lt 3 pt 6, "erms_test.dat" w lp lt 1 pt 6

f:id:mmitou:20120920183611p:image

トラックバック - http://d.hatena.ne.jp/mmitou/20120920/1348134209
リンク元