パラメータの事後分布計算を実装した

Haskell 機械学習 PRML


PRML3.3.1のパラメータの分布を読んで、パラメータの事後分布計算を実装した。

事前分布、事後分布、尤度関数の関係


観測データ \bf{D}が得られた時に、そのデータがパラメータ \bf{w}に依存していたとする。
この時、観測データ \bf{D}があるパラメータ \bf{w}の条件下で得られる確率はP(\bf{D}|\bf{w})である。
P(\bf{D}|\bf{w})は \bf{w}の値によって変わる関数なので、尤度関数と呼ばれる。
また、パラメータ \bf{w}の確率分布 P(\bf{w})が事前分布であり、
観測データ \bf{D}が得られた後でわかるパラメータ \bf{w}の確からしさは事後分布 P(\bf{w}|\bf{D})である。
これら事前分布、事後分布、尤度関数は、ベイズの定理より以下の関係にある。(PRML 1.43式参照)
P(\bf{w}|\bf{D}) = \frac {P(\bf{D}|\bf{w}) P(\bf{w})} {p(\bf{D})}
つまり、事後分布 P(\bf{w}|\bf{D}) \propto 尤度関数 P(\bf{D}|\bf{w}) \times 事前分布 P(\bf{w})である。

尤度関数


観測データ \bf{D}のそれぞれの値が独立同分布である時、観測データ \bf{D}の得られる確率はは各データが得られる確率の積となる。
例えば、独立同分布で得られる値 y_1, y_2の同時確率は P(y_1,y_2) = P(y_1)P(y_2)である。
同様に、パラメータ \bf{w}の条件下で観測データ \bf{D}の得られる確率 P(\bf{D}|\bf{w})は、
それぞれの観測値が得られる確率の積となるので以下のように表現できる。
P(\bf{D}|\bf{w}) = \prod_{i=1}^{N} P(t_i|x_i,\bf{w})


特に、PRML3.1で導入した線形基底関数モデルを採用した場合に、入力 \bf{X} = \{x_1, ..., x_N\}と対応する目標値 \bf{t} = {t_1, ... , t_N}が得られた場合の尤度関数は、
P(\bf{t}|\bf{X},\bf{w},\beta) = \prod_{n=1}^{N}N(t_n|\bf{w}^{T}\bf{\phi}(x_n),\beta^{-1}) ... (3.10)
なお、入力変数のモデル化はしないので、以降 \bf{X}は省略する。

パラメータの分布


式(3.10)の尤度関数 P(\bf{t}|\bf{w})の指数部分が \bf{w}の2次形式であるため、
この尤度関数に対応する共役事前分布が以下のようなガウス分布で与えられる事に注目する。
P(\bf{w}) = N(\bf{w}|\bf{m_0}, \bf{S_0})
事後分布は尤度関数と事前分布の積に比例するので、以下のように書ける。
P(\bf{w}|\bf{t}) \propto \exp \Bigl(-\frac {1}{2} (\bf{w} - \bf{m_0})^{T}\bf{S}_{0}^{-1}(\bf{w} - \bf{m_0}) - \sum_{n=1}^{N}(t_{n}-\bf{w}^{T}\bf{\phi}(x_n))^2 \Bigr)
この式の指数部分を平方完成して、観測データ取得後の事後分布 P(\bf{w}|\bf{t}) = N(\bf{w}|\bf{m_N},\bf{S_N})も指数部分を展開して平方完成して比較すると、
以下のような期待値 \bf{m_N}と共分散 \bf{S_{N}^{-1}}が得られる。
\bf{m_N} = \bf{S_N}(S_{0}^{-1}\bf{m_{0}} + \beta \bf{\Phi}^{T}\bf{t})
\bf{S_{N}^{-1}} = \bf{S_{0}^{-1}} + \beta \bf{\Phi}^{T} \bf{\Phi}


ここで、 P(\bf{w}|\alpha) = N(\bf{w}|\bf{0}, \alpha^{-1}\bf{I})を事前分布とした場合に、
目標値 \bf{t}と入力値 \bf{X}が与えられた場合の、事後分布 P(\bf{w}|\bf{t})を計算した結果を示す。
目標値 \bf{t}は、 t = 3.14 * x + 1.23で計算したtに-1から1の範囲で得られるランダムな値を加えて作成した。
以下のグラフは、横軸に w_0、縦軸に[w_1]を取り、色で座標上の確率を表したものである。

w_0 = 1.23, w_1 = 3.14の座標で確立が集中的に高くなっていることがわかる。

この計算を行うプログラムを以下に示す。

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)

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

-- 座標を文字列にする
from3dToString :: (Show a1, Show a2, Show a3)
                  => [(a1, a2)] -> [a3] -> String
from3dToString xs ts = concat $
                      zipWith (\(x,y) z ->    (show x) ++ " "
                                           ++ (show y) ++ " "
                                           ++ (show z) ++ "\n")
                      xs ts

-- 多変量ガウス分布
-- mu    : D次元平均ベクトル 実装上はD×1の行列
-- sigma : D×Dの共分散行列
-- x     : D次元ベクトルx   実装上はD×1の行列
-- gaussDist :: Matrix a -> Matrix a -> Matrix a -> a
gaussDist mu sigma x = c * exp a
  where
    d = fromIntegral $ cols x
    c = 1.0 / (((2 * pi) ** (d / 0.5)) * ((det sigma) ** 0.5))
    y = x - mu
    a = -0.5 * mahalanobisDistance mu sigma x
    mahalanobisDistance avrg covariance x
      = ((trans $ x - avrg)
         `multiply` (inv covariance)
         `multiply` y) @@> (0,0)
      where
        y = x - avrg

-- 多項式基底関数
polynomialBasis :: Floating a => a -> a -> a
polynomialBasis n x = x ** n

-- 基底関数のリストと入力変数xからベクトル?(x)を計算する。
-- 実装上は M×1の行列
basisFuncVector basiss x = (m >< 1) $ [f x | f <- basiss]
  where
    m = length basiss

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

-- 事後分布
postDist alpha beta xs ts = gaussDist m s
  where
    phi = designMatrix (map polynomialBasis [0,1]) xs
    phiTphi = (trans phi) `multiply` phi
    n = (cols phiTphi)
    alphaI = diag . fromList $ take n $ repeat alpha
    betaI  = diag . fromList $ take n $ repeat beta
    invS = alphaI + betaI `multiply` phiTphi
    s = inv invS
    m = betaI `multiply` s `multiply` (trans phi)
        `multiply` ((length ts >< 1) ts)

-- 事後分布を受け取り、関数内で生成した各座標での確率を計算する。
-- 座標と確率を文字列にして返す。
calcDist dist = str
  where
    xy  = [(x,y) :: (Double, Double)
           | x <- linearScale 100 (-5, 5)
           , y <- linearScale 100 (-5, 5)]
    xy' = map ((2 >< 1) . (\(x,y) -> [x,y])) xy
    z   = map dist xy'
    str = from3dToString xy z

-- 事後分布を計算してファイルに出力する。
-- xの個数が少ないと1を超えた確率が出る場合がある
runPostDistCalc = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 10)
      f  = \x -> 3.14 * x + 1.23
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map f xs) gen
      alpha = 1.0
      beta  = 1.0
      g = postDist alpha beta xs ts
      str = calcDist g
  writeFile "gaussDist.dat" str
  return ()

目標値tと入力値xを一つずつ入力して事後分布を更新していく計算方法を考える事も出来る。
以下のアニメーションでは、更新された事後分布 P(\bf{w})が徐々に w_0 = 1.23, w_1 = 3.14に収束していく事がわかる。

このアニメーションの確率分布は以下のプログラムで計算した。

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

-- ガウス分布の共分散と期待値を保持するデータ構造
data GaussDist = GaussDist {
    mu    :: Matrix Double
  , sigma :: Matrix Double
} deriving (Show, Eq)

-- パラメータWとその確率を保持するためのデータ構造
data Pos = Pos Double Double Double deriving Eq 
instance Show Pos where
  show (Pos x y z) = show x ++ " " ++
                     show y ++ " " ++
                     show z ++ "\n"

-- 多変量ガウス分布
-- mu    : D次元平均ベクトル 実装上はD×1の行列
-- sigma : D×Dの共分散行列
-- x     : D次元ベクトルx   実装上はD×1の行列
-- gaussDist :: Matrix a -> Matrix a -> Matrix a -> a
gaussDist mu sigma x = c * exp a
  where
    d = fromIntegral $ cols x
    c = 1.0 / (((2 * pi) ** (d / 0.5)) * ((det sigma) ** 0.5))
    y = x - mu
    a = -0.5 * mahalanobisDistance mu sigma x
    mahalanobisDistance avrg covariance x
      = ((trans $ x - avrg)
         `multiply` (inv covariance)
         `multiply` y) @@> (0,0)
      where
        y = x - avrg

-- 座標wの確率を計算する
-- wはパラメータベクトル 実装上はD×1の行列
probability :: GaussDist -> Matrix Double -> Double
probability dist w = gaussDist m s w
  where
    m = mu dist
    s = sigma dist

-- 多項式基底関数
polynomialBasis :: Floating a => a -> a -> a
polynomialBasis n x = x ** n

-- 確率分布を更新する。演習3.8により共分散と期待値を更新する
-- 第一引数 係数α
-- 第二引数 係数β
-- 第三引数 基底関数にx_n+1を適用した?_n+1縦ベクトル(実装上はM行1列の行列)
-- 第四引数 新しい目標値t_n+1
-- 第五引数 古いガウス分布
update :: Double -> Double -> Matrix Double ->
          Double -> GaussDist -> GaussDist
update alpha beta phi t dist = dist'
  where
    m = mu    dist
    s = sigma dist
    invS' = (inv s) + ((1><1) [beta])*(phi `multiply` (trans phi))
    lhs   = (inv s) `multiply` m
    rhs   = ((1><1) [beta * t]) * phi
    m'    = (inv invS') `multiply` (lhs + rhs)
    dist' = GaussDist m' (inv invS')

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

-- 訓練データを作成する
-- 関数t(x) = 3.14 x + 1.23から 0 <= x < 5の範囲で生成した値に
-- -1.0 から 1.0の範囲の値をランダムに加えた(x,t)のリストを返す。
makeTrainingData :: IO [(Double, Double)]
makeTrainingData = do
  gen <- getStdGen
  let xs = linearScale 1000 (0, 5)
      f  = \x -> 3.14 * x + 1.23
      ts = addNoise (-1.0 ::Double, 1.0::Double) (map f xs) gen
      xts =  zip xs ts
  return xts

-- 事後分布の逐次更新計算を行う
-- 入力された目標値のリストから目標値を一つずつ取り出して事後分布を更新していく。
sequentialUpdateCalc :: Double -> Double ->
                        [(Double -> Double)] ->
                        GaussDist -> [(Double, Double)] -> [[Pos]]
sequentialUpdateCalc alpha beta basiss priorDist xts =
  seqCalc priorDist xts
  where
    update' = update alpha beta
    -- -5 < w0 < 5 , -5 < w1 <5の範囲にある確率P(w1,w2)を全て計算する
    w = [(2 >< 1) [w0, w1] :: Matrix Double
        | w0 <- linearScale 100 (-5, 5)
        , w1 <- linearScale 100 (-5, 5)]
    -- 戻り値で座標wを返す際に行列のままだと表示しづらいので、リストに変換する
    wl = map (concat . toLists) w
    -- 実際の逐次計算は以下の関数で行う。
    -- (x,t)を一つずつ
    seqCalc _ [] = []
    seqCalc priorDist ((x,t):xts) = poss:(seqCalc postDist xts)
      where
        -- 基底関数にxを適用したベクトル?nを計算する
        phi  = (2 >< 1) [f x | f <- basiss]
        -- 全ての(w0,w1)に対して、その座標での確率を計算する
        ps   = map (probability priorDist) w
        -- w0,w1,P(w0,w1)からなるリストを作成する
        poss = zipWith (\ (w0:w1:[]) p -> Pos w0 w1 p) wl ps
        -- 事前分布から事後分布を計算する(演習3.8)
        postDist = update' phi t priorDist

-- 文字列のリストの各要素をファイルに出力する。        
writeSequential name strs = writeSeq 0 strs
  where
    writeSeq _ [] = return ()
    writeSeq n (str:strs) = do
      let fileName = name ++ show n ++ ".dat"
      writeFile fileName str
      writeSeq (n + 1) strs

-- 逐次更新処理を実行する      
runSequentialUpdate = do
  -- 訓練データを作成する
  xts <- makeTrainingData
  let basiss  = map polynomialBasis ([0,1] :: [Double])
      -- 事前分布 P(w) = N(w|m0,s0)
      m0 = (2 >< 1) [0,0]
      s0 = diag . fromList $ [1,1]
      priorDist = GaussDist m0 s0
      -- 更新されていく確率分布を計算する
      posss = sequentialUpdateCalc 1 1 basiss priorDist xts
      -- 計算した確率分布を文字列にする
      strs = map (concat . (map show)) posss
  writeSequential "update" strs
  return ()