Practical Probabilistic Programming with Monads を読む (1) : とりあえずサイコロを振る

Practical Probabilistic Programming with Monads (PDF) という論文を紹介して頂いたので少し読んでみます．関連のコードが adscib/monad-bayes: A library for probabilistic programming in Haskell. にあり，今回は準備編ということで背景の紹介など全部すっ飛ばしてサイコロを振ってみましょう．統計界隈の読者を想定してステップバイステップでいきます．環境は Ubuntu です．殆どの部分は他のディストリビューションや Mac とかでも動くと思います．apt だけ適宜読み替えてください．

準備

• Haskell と stack 導入
  • 色々紹介記事があるのでよしなにやってください．
• 必要なライブラリ導入
  • 覚えてる範囲では $ sudo apt install liblapack-dev libblas-dev
• $ git clone https://github.com/adscib/monad-bayes
  • 手許に持たなくとも動かすこと自体はできるのですが，色々と便利なのでローカルに落としておきます．
  • 「論文内のコードについては haskell2015 ブランチを見てね」と書いてあるのですが，ざっと見た感じ随分粗いコードな上論文と仕様がちがうところもあり，とりあえず master を使います．
  • ついでに tag とかも作っとくと便利．
• stack.yaml を編集して $ stack haddock
  • 取ってきたmonad-bayes のディレクトリ内で実行してください．時間がかかると思います．
  • stack.yaml の編集は必須ではありませんが，resolver が結構古くなっているので自分の環境に合わせて，例えば今なら resolver: lts-8.13 というふうに書き換えるとやりやすいと思います（問題なくビルドは通ります）．
  • これで stack haddock --open でドキュメンテーションが見られるようになると思います．

雛形

好きな名前，例えば mbplay をつけて $ stack new mbplay しましょう．以降このディレクトリ内で作業します．

stack と cabal の準備

この monad-bayes は stackage になく， まだ hackageにもあがっていない¹ ので我々が自分で stack さんに探しに行く場所を教えてあげる必要があります．さっき出来た stack.yaml のpackages: のところに次の下三行を追記します：

stack.yaml

packages:
- '.'
- location:
    "../monad-bayes" # さっき clone したディレクトリ
  extra-dep: true    # 依存先（本体の一部じゃなくて）として扱ってもらう

続いて，そこにある monad-bayes というパッケージを使うよ，ということを cabal に書き込みます．

mbplay.cabal

library  
  hs-source-dirs:      src
  exposed-modules:     Lib
  build-depends:       base >= 4.7 && < 5
                     , monad-bayes

こんな感じ． executable mbpy-exe の方にも同様に monad-bayes を書き足しましょう²．

実際にインポートしてみて

Lib.hs

module Lib (someFunc) where
import Control.Monad.Bayes.Simple

stack build が通るのを確認してみてください．

Dice の用意

あとはサイコロを振って終わりです． 例は monad-bayes の models/ にあるのでここから Dice.hs を拝借します．

terminal

# mbplay に今いるとして
cp ../monad-bayes/models/Dice.hs src/

.cabal も編集して Dice っていうモジュールを export してるで，ということを書きます

mbplay.cabal

library
  hs-source-dirs:      src
  exposed-modules:     Dice

振る

とりあえず 6面ダイスを振りたい．関連する部分は以下です．

Dice.hs

-- | A toss of a six-sided die.
die :: MonadDist d => d Int
die = categorical $ zip [1..6] (repeat (1/6))

-- | A sum of outcomes of n independent tosses of six-sided dice.
dice :: MonadDist d => Int -> d Int
dice 1 = die
dice n = liftM2 (+) die (dice (n-1))

die は論文では面の数を受けて Dist Int を返していましたがここでは 6面固定，MonadDist³ という型クラスの d に Int を包んで返しています． 1 から 6 まで確率が 1/6 ずつの categorical distribution として作られています． dice で n 回振れますが，これは素直な再帰．

この distribution からサンプルを取りたいのですが，元論文には sample :: Samplable a => StdGen -> d a -> a という分かりやすい関数があったのが色々と変更が入っています（この辺については次回以降にできれば書きたい）．結論だけいえばSamplerIO を使うのが今回は正解で，初期化もして乱数を使ってとって来てくれる sampleIO :: SamplerIO a -> IO a が使えます． instance MonadDist SamplerIO なので die :: d Int が直接いけまして

Main.hs

main = sampleIO die >>= print

でさいころを振れます．

折角だから 「20回振る」のを 5000 回やってみます．

app/Main.hs

module Main where

import Dice
import Control.Monad.Bayes.Sampler

import Data.List (group, sort)
import Control.Monad (replicateM)

main :: IO ()
main = do
    results <- sampleIO . replicateM 5000 $ dice 20
    mapM_ putStrLn . hist . map length . group . sort $ results

-- create a very basic histogram. 
-- `div` 2 . (+1) is just there to make the bars shorter.
hist :: [Int] -> [String]
hist = map (flip replicate '#' . (`div` 2) . (+1))

$ stack build && stack exec mbplay-exeの結果は

やりました．

---

これでは全然パッケージのことがわからないので続く…予定．

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

1. Readme にある程度安定したらアップするよって書いてありますがまだのようです ↩

2. もちろんライブラリ側か Main.hs のどっちかでしか使わないならそちらの方だけでいいです． ↩

3. class (...) => MonadDist m where categorical :: [(a,CustomReal m)] -> m a ... ↩