遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

MCMCをやってみる

Haskell 確率的プログラミング
  • こちらのページをやってみる
  • やるために、stackで環境を整える
  • stack new hogeして、そのcabalファイルのmainの必要パッケージを以下のようにする
executable pfptest-exe
  hs-source-dirs:      app
  main-is:             Main.hs
  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
  build-depends:       base
                     , pfptest
                     , probability
                     , histogram-fill
                     , random-fu
                     , statistics
                     , json
                     , vector
                     , random-source
                     , mtl
> 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
    • Haskell Control.Monadで検索するとHackageページが見つかり、その一番上の「バー」にbase-4.9.1.0: Basic librariesと出る。このbaseはわざわざcabalに書かなくても、importできる
    • 次にHaskell Data.Histogramとすると、histogram-fill-0.8.5.0: Library for histograms creation.と出る。このhistogram-fillまでがパッケージ名なのでこれを書く。以下同様
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
  • こうやると、test-data.jsonという名前のファイルができて、その中身は
{"x":[-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15],"y":[-73.70187524627536,-53.32544108228022,-85.5561593826951,-53.2923266140698,-59.30169338011609,-34.71713105945278,-64.23269931716732,-37.70895778207112,-32.49602734293523,-12.5987398418178,-32.78105748127745,-17.02959946540419,-12.709949190288553,-24.007282280667503,0.8279906165058932,-9.244400639134525,17.09474416823518,23.228482221669115,14.690049579435886,29.41188167743224,39.42186361486318,31.413353262653644,33.189437506871066,31.840776644320762,42.815737959293465,54.51844255645821,34.5842642875685,65.82943043766038,70.11853692362467,56.332662810976046,67.91804611706799]}
  • 乱雑項を含むデータxs,ysができたので、回帰直線の傾きの尤度を計算してみる
    • a,b,xが与えられた時に、観察値yの生起確率を乱雑項がsd標準偏差の正規分布と仮定したときの対数尤度を計算し、それをすべてのxs,ysペアについて行って足し合わせる関数が loglikelihood
    • それを固定のb,固定のsdのもとでaだけ3から7まで振ってjsonファイルに書き出させるのがexportLikelihoodA
:{
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]
exportLilekiooodA
  • logDensityという関数は以下のように、連続分布と、そのサポート上の値とをとってその対数確率密度・対数尤度を返す関数
> :type logDensity
logDensity :: ContDistr d => d -> Double -> Double
>||
-normalDistrは平均とSDをとって、具体的な正規分布を返す関数
>|haskell|
> :type logDensity
logDensity :: ContDistr d => d -> Double -> Double
  • precictは、２つの係数と台の値とをとって、その予測値を返す関数
> :type logDensity
logDensity :: ContDistr d => d -> Double -> Double

{"x":[3,3.05,3.0999999999999996,3.1499999999999995,3.1999999999999993,3.249999999999999,3.299999999999999,3.3499999999999988,3.3999999999999986,3.4499999999999984,3.4999999999999982,3.549999999999998,3.599999999999998,3.6499999999999977,3.6999999999999975,3.7499999999999973,3.799999999999997,3.849999999999997,3.899999999999997,3.9499999999999966,3.9999999999999964,4.049999999999996,4.099999999999996,4.149999999999996,4.199999999999996,4.249999999999996,4.299999999999995,4.349999999999995,4.399999999999995,4.449999999999995,4.499999999999995,4.5499999999999945,4.599999999999994,4.649999999999994,4.699999999999994,4.749999999999994,4.799999999999994,4.849999999999993,4.899999999999993,4.949999999999993,4.999999999999993,5.049999999999993,5.0999999999999925,5.149999999999992,5.199999999999992,5.249999999999992,5.299999999999992,5.349999999999992,5.3999999999999915,5.449999999999991,5.499999999999991,5.549999999999991,5.599999999999991,5.649999999999991,5.69999999999999,5.74999999999999,5.79999999999999,5.84999999999999,5.89999999999999,5.9499999999999895,5.999999999999989,6.049999999999989,6.099999999999989,6.149999999999989,6.199999999999989,6.2499999999999885,6.299999999999988,6.349999999999988,6.399999999999988,6.449999999999988,6.499999999999988,6.549999999999987,6.599999999999987,6.649999999999987,6.699999999999987,6.749999999999987,6.7999999999999865,6.849999999999986,6.899999999999986,6.949999999999986,6.999999999999986],"loglikelihood":[-159.40621882876707,-157.14658427153978,-154.94894971431242,-152.81331515708507,-150.73968059985776,-148.72804604263047,-146.77841148540307,-144.89077692817577,-143.06514237094842,-141.30150781372106,-139.59987325649377,-137.9602386992664,-136.3826041420391,-134.86696958481176,-133.41333502758442,-132.02170047035708,-130.69206591312974,-129.4244313559024,-128.21879679867504,-127.07516224144769,-125.99352768422038,-124.97389312699303,-124.01625856976571,-123.12062401253839,-122.28698945531103,-121.51535489808369,-120.80572034085631,-120.158085783629,-119.57245122640164,-119.04881666917427,-118.58718211194697,-118.18754755471959,-117.8499129974923,-117.5742784402649,-117.3606438830376,-117.20900932581023,-117.1193747685829,-117.09174021135556,-117.12610565412821,-117.22247109690086,-117.38083653967347,-117.60120198244616,-117.8835674252188,-118.22793286799147,-118.63429831076412,-119.10266375353675,-119.63302919630945,-120.22539463908204,-120.87976008185471,-121.59612552462738,-122.37449096739998,-123.21485641017267,-124.11722185294529,-125.08158729571792,-126.10795273849061,-127.19631818126324,-128.3466836240359,-129.5590490668085,-130.83341450958116,-132.16977995235382,-133.56814539512644,-135.02851083789912,-136.55087628067176,-138.13524172344435,-139.78160716621701,-141.48997260898966,-143.2603380517623,-145.09270349453496,-146.98706893730758,-148.9434343800802,-150.96179982285287,-153.04216526562553,-155.1845307083981,-157.3888961511708,-159.6552615939434,-161.98362703671606,-164.37399247948872,-166.82635792226137,-169.34072336503394,-171.91708880780664,-174.5554542505793]}
|