2017-05-26

事前分布、事後分布、MCMC

事前分布は、パラメタに分布を仮定し、その尤度が計算できることとして決まる
- a,b,sdに一様分布を設定し、何か具体的な値に対して、その総合的な対数尤度をlogpriorが返す

:{
let logprior :: Double -> Double -> Double -> Double;logprior a b sd = sum [ap, bp, sdp]
       where
              ap  = logDensity (uniformDistr 0 10) a
              bp  = logDensity (normalDistr 0 5) b
              sdp = logDensity (uniformDistr 0 30) sd
:}

事後確率は、a,b,sdに対して、事前尤度に観察尤度を掛けあわせればよい。対数尤度なら加算
- 事前分布も尤度関数もa,b,sdの関数として作ってある

let logposterior :: Double -> Double -> Double -> Double; logposterior a b sd = loglikelihood a b sd + logprior a b sd

let nIters :: Int; nIters = 100000

- seedをとってMT擬似欄数列発生器を作り、ランダムなevalStateをさせる。平均0,指定sdでの正規乱数の長さnItersの列を取り出す

:{
let normalisedProposals :: Int -> Double -> Int -> [Double]; normalisedProposals seed sigma nIters =
       evalState (replicateM nIters (sample (Normal 0.0 sigma)))
       (pureMT $ fromIntegral seed)
:}

- Accept/Rejectの判別のための一様乱数列の発生

:{
let acceptOrRejects :: Int -> Int -> [Double]; acceptOrRejects seed nIters =
       evalState (replicateM nIters (sample stdUniform))
       (pureMT $ fromIntegral seed)
:}

- 採否の閾値を計算。新たな係数セットを現在の係数セットと比べる。比べるにあたって、事後分布としての尤度で返す。次のoneStepではこの尤度を確率的に採否するべく、上で発生した一様乱数の値と比較する

:{
let acceptanceProb :: Double -> Double -> Double -- ^ Proposal parameters
       -> Double -> Double -> Double -- ^ Old parameters 
       -> Double; acceptanceProb a' b' sd' a b sd =
       exp ((logposterior a' b' sd') - (logposterior a b sd))
:}

- 上で説明した通り

:{
let oneStep :: (Double, Double, Double, Int)    -- ^ Current location
       -> (Double, Double, Double, Double) -- ^ Delta location and acceptance probability
       -> (Double, Double, Double, Int) ; oneStep (a, b, sd, nAccs) (da, db, dsd, acceptOrReject) =
       if acceptOrReject < acceptanceProb (a + da) (b + db) (sd + dsd) a b sd
       then (a + da, b + db, sd + dsd, nAccs + 1)
       else (a, b, sd, nAccs)
:}

- 更新するための諸関数ができたら、変数初期値のセットを作り、採否をかませた更新関数を使いつつ、scanl関数にて、連鎖・更新しながらそれを全部記録していく(scanlはこちら)
以下はただのコード記録

import Control.Monad
import Control.Monad.State (evalState)
import Data.Histogram (asList)
import Data.Histogram.Fill
import Data.Histogram.Generic (Histogram)
import Data.List
import Data.Random
import Data.Random.Source.PureMT
import Statistics.Distribution hiding (mean)
import Statistics.Distribution.Normal
import Statistics.Distribution.Uniform
import Text.JSON
import qualified Data.Histogram as H
import qualified Data.Vector.Unboxed as V

let predict :: Double -> Double-> Double -> Double; predict a b x = a * x + b
    
let trueA, trueB, trueSd :: Double; trueA = 5;  trueB = 0; trueSd = 10

let nSamples :: Int; nSamples = 31

let xs :: [Double]; xs = take nSamples $ [-15.0 ..]

let arbSeed :: Int; arbSeed = 4022 -- Arbitrary seed for generating random numbers

:{
let ys :: [Double]; ys = zipWith (+) noise (map (predict trueA trueB) xs)
       where noise = evalState (replicateM nSamples (sample (Normal 0.0 trueSd)))
              (pureMT $ fromIntegral arbSeed)  
:}


let exportTestData :: IO ();exportTestData = writeFile "test-data.json" . encode . toJSObject $ [("x", xs), ("y", ys)]

exportTestData

:{
let loglikelihood :: Double -> Double -> Double -> Double; loglikelihood a b sd = sum $ zipWith lh xs ys
       where lh x = logDensity (normalDistr (predict a b x) sd)
:}

let exportLikelihoodA :: IO (); exportLikelihoodA = writeFile "log-likelihood.json" . encode . toJSObject $ (\as -> [ ("x",as), ("loglikelihood", map (\a -> loglikelihood a trueB trueSd) as )]) [3, 3.05 .. 7::Double]

:{
let logprior :: Double -> Double -> Double -> Double;logprior a b sd = sum [ap, bp, sdp]
       where
              ap  = logDensity (uniformDistr 0 10) a
              bp  = logDensity (normalDistr 0 5) b
              sdp = logDensity (uniformDistr 0 30) sd
:}
let logposterior :: Double -> Double -> Double -> Double; logposterior a b sd = loglikelihood a b sd + logprior a b sd

let nIters :: Int; nIters = 100000

:{
let normalisedProposals :: Int -> Double -> Int -> [Double]; normalisedProposals seed sigma nIters =
       evalState (replicateM nIters (sample (Normal 0.0 sigma)))
       (pureMT $ fromIntegral seed)
:}

:{
let acceptOrRejects :: Int -> Int -> [Double]; acceptOrRejects seed nIters =
       evalState (replicateM nIters (sample stdUniform))
       (pureMT $ fromIntegral seed)
:}
:{
let acceptanceProb :: Double -> Double -> Double -- ^ Proposal parameters
       -> Double -> Double -> Double -- ^ Old parameters 
       -> Double; acceptanceProb a' b' sd' a b sd =
       exp ((logposterior a' b' sd') - (logposterior a b sd))
:}

:{
let oneStep :: (Double, Double, Double, Int)    -- ^ Current location
       -> (Double, Double, Double, Double) -- ^ Delta location and acceptance probability
       -> (Double, Double, Double, Int) ; oneStep (a, b, sd, nAccs) (da, db, dsd, acceptOrReject) =
       if acceptOrReject < acceptanceProb (a + da) (b + db) (sd + dsd) a b sd
       then (a + da, b + db, sd + dsd, nAccs + 1)
       else (a, b, sd, nAccs)
:}