PRML３章の線形回帰を実装した

PRML３章の線形回帰を実装した。ベイズ線形回帰は未実装。

中谷さんのgihyo.jp連載記事「機械学習　はじめよう」
第11回 "線形回帰を実装してみよう"が非常に参考になりました。

http://gihyo.jp/dev/serial/01/machine-learning/0011?page=1

実装は以下の通り。

import System.Random
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple
import Control.Applicative

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

-- 座標を文字列にする
pointsToString :: (Show a) => [a] -> [a] -> String
pointsToString xs ys =
  concat $
  zipWith (\x y -> show x ++ " " ++ show y ++ "\n") xs ys

-- 座標を文字列にする
to3dPosString xs ts = concat $
                      zipWith (\(x,y) z ->    (show x) ++ " "
                                           ++ (show y) ++ " "
                                           ++ (show z) ++ "\n")
                      xs ts

-- ガウス基底関数
gaussianBasis s mu x = exp $ (- (x - mu) ** 2.0) / s

-- 多項式基底関数
polynomialBasis n x = x ** n

-- 二次元用多項式基底関数
-- 多次元の場合にこのような基底関数でよいのかどうかは確認していない。
polynomialBasis' n (x,y) = x ** n + y ** n

-- 基底関数の配列と入力変数xから計画行列Φを計算する
designMatrix basiss xs = (n >< m) [f x | x <- xs, f <- basiss]
  where
    n = length xs
    m = length basiss

-- 正則化係数、計画行列と目標変数の配列から
-- 基底関数の線形結合で用いるパラメータ w を計算する
parameters lambda phi ts = w
  where
    -- 配列を縦ベクトルに変換する
    t       = (length ts >< 1) ts 
    m       = cols phi
    -- (lambda * 単位行列) を計算する
    lambdaI = fromLists $ map (map (*lambda)) (toLists $ ident m)
    -- 計画行列の転置行列を計算する
    phiT = trans phi
    -- w = ((λ I + (Φ^T)Φ)^(-1))(Φ^T)t <= PRML P.143 (3.28)式より
    w = (inv (lambdaI + phiT `multiply` phi))
        `multiply` phiT `multiply` t

-- 線形回帰
-- 正則化係数, パラメータの要素数, 基底関数, 入力値xs,目標値tsから、
-- 線形回帰によって推定された関数fを返す。
linearRegression lambda m basis xs ts = f
  where
    -- 基底関数ベクトルを作る
    basiss = (\_ -> 1.0) : (map (basis . fromIntegral) [1..m])
    -- 計画行列を計算する
    dm = designMatrix basiss xs
    -- パラメータベクトルwを計算してパラメータのリストにする
    ws = concat . toLists $ parameters lambda dm ts
    -- 訓練データから推測される関数を返す
    f  = \x -> sum $ zipWith (*) ws $ map (\g -> g x) basiss

-- xが一次元の場合の計算
test2d = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 2*pi)
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map sin xs) gen
      f  = linearRegression 0.2 100 (gaussianBasis 0.2) xs ts
      f' = linearRegression 0.2 10 polynomialBasis xs ts
      ys  = map f xs
      ys' = map f' xs
      xtstr  = pointsToString xs ts
      xystr  = pointsToString xs ys
      xystr' = pointsToString xs ys'
  writeFile "train.dat"      xtstr
  writeFile "gaussian.dat"   xystr
  writeFile "polynomial.dat" xystr'


-- xが二次元の場合の計算
test3d = do
  gen <- getStdGen
  let xs = linearScale 50 (0, 2*pi)
      ys = linearScale 50 (0, 2*pi)
      xs' = [(x,y) | x <- xs, y <- ys]
      org = \(x, y) -> cos(0.5 * (x - pi)) * cos(0.5 * (y -pi))
      ts = addNoise (-0.2::Double, 0.2::Double) (map org xs') gen
      f = linearRegression 0.2 10 polynomialBasis' xs' ts
      zs  = map f xs'
      xtstr  = to3dPosString xs' ts
      xzstr  = to3dPosString xs' zs
  writeFile "train3d.dat"      xtstr
  writeFile "gaussian3d.dat"   xzstr
  return ()

main = do
  test2d
  test3d

xが1次元の場合の線形回帰の計算結果は以下のとおり。
ガウス基底関数の分散を大きくすると、もっとスムーズな関数を得る事は出来る。

xが2次元の場合の線形回帰の計算結果は以下の通り。
2次元での多項式基底関数のとり方が正しいかどうか要確認。