Haskellでメトロポリスヘイスティング

  • こちらには、前の記事で勉強したDistのことが書いてあるのだが、それを使ってメトロポリスヘイスティングを実装している。
  • 以下の通り
  • このHaskellコードを読んでいくことにする
  • メトロポリスヘイスティングでは(関数 mh として作成されている)
    • 事象の確率は計算できる確率分布があるときに、モンテカルロで事象を発生させ、その採否を確率的に決めることで、分布からの事象の集合を発生させる方法
    • 以下では、そのような確率を計算できる分布を第三引数Dist aで与えている
    • 第1、第2引数は、初期値依存性が高いと想定される初期の発生事象として捨てる回数(burn-in回数)と、その後の、蓄積して分布乱数として採用する回数とを与えている
    • 返り値は、事象集合をリストで持つDist [a] である
mh :: Int      -- burn-in 期間
   -> Int      -- サンプリング数
   -> Dist a   -- `Conditional`を含む事前分布
   -> Dist [a] -- 事後分布からサンプリングした点列上の確率分布
mh dc tc d = fmap (map fst . drop dc) $ proposal >>= iterate (dc + tc)
  where
    proposal = prior d
    iterate 0 _ = pure []
    iterate n (x,s) = do
      (y,r) <- proposal
      accept <- bernoulli $ min 1 (r / s)
      let next = if accept then (y,r) else (x,s)
      rest <- iterate (n-1) next
      return $ next:rest
    • "proposal >>= iterate (dc + tc)"
      • proposalという、「確率計算できる確率分布」が、Dist a になっているので、そのaをモナドであるDistの箱から出して
      • iterate (dc + tc) a とする。これは、バーンインとその後の集計との総合計の回数の繰り返しをすることを言っている
      • iterateという関数はwhere以下に定義されている
      • こうして発生した(dc + tc)個の乱数のうち、前dc個を捨てて(drop dc)、その上で、乱数には、事象の値とその生起確率とがペアで入っているので、その事象値の方(fst 第一要素)を取り出すという処理を、dc個すててtc個になったペアにmap適用する
      • ただし、そうするためには、$以下の部分がDistでくるまれているので、外して適用して、くるみ直すためにfmapを使っている
    • where以下は?
      • 全dc + tc回の繰り返しにあたって、繰り返し0回の状態は、空リストに相当する状態であって、それをDistとしてハンドリングするために pure をかぶせている( iterate0 _ pure [])
      • 1つインスタンスを採用するたびに、そのインスタンス(事象と確率のペア)とをnextとして、返却リストDistの1要素として追加している(next:rest)
      • 追加したら、iterate (n-1) nextとして、そのときのインスタンスを受け取りつつ、1回少ない繰り返しをする状態に移行する
      • どうしてそのときのインスタンスを渡すかというと、次のインスタンス候補を、オリジナルの入力分布から乱択し、その新しいインスタンスの採否を、一つ前に採用されたインスタンスと比較して決めるから。その採否は、両者の生起確率 r s を比較し、より高確率なインスタンスが新しい方ならば、必ず採択し、そうでなければ、確率 r/s で採択し、採択しない時は、そのときの値を再度採用する、となっている
      • 採択したなら、次に回に回す「そのとき」は更新し、採択しないなら更新しない
    • このmhは「メトロポリスヘイスティングサンプリングできるよ」という状態を作る
    • 以下にあるように、そこにRandom.sampleを適用することで、実際にバーンイン回だけ回して捨てて、その後の値ベクトルを包んだIOモナドを作る
  td <- trainingData -- trainingDataIOモナドとして作るのがtrainingData、そこから中身を取り出して、乱点リストをtdとする