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

Practical Probabilistic Programming with Monads ([PDF](http://mlg.eng.cam.ac.uk/pub/pdf/SciGhaGor15.pdf)) という論文を紹介して頂いたので少し読んでみます．関連のコードが [adscib/monad-bayes: A library for probabilistic programming in Haskell.](https://github.com/adscib/monad-bayes) にあり，今回は準備編ということで背景の紹介など全部すっ飛ばしてサイコロを振ってみましょう．統計界隈の読者を想定してステップバイステップでいきます．環境は 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](https://www.stackage.org/lts-8.13) になく， まだ [hackage](http://hackage.haskell.org/)にもあがっていない[^hack] ので我々が自分で stack さんに探しに行く場所を教えてあげる必要があります．さっき出来た `stack.yaml` の`packages:` のところに次の下三行を追記します：

```yaml:stack.yaml
packages:
- '.'
- location:
    "../monad-bayes" # さっき clone したディレクトリ
  extra-dep: true    # 依存先（本体の一部じゃなくて）として扱ってもらう
```
続いて，そこにある `monad-bayes` というパッケージを使うよ，ということを cabal に書き込みます．

```cabal:mbplay.cabal
library  
  hs-source-dirs:      src
  exposed-modules:     Lib
  build-depends:       base >= 4.7 && < 5
                     , monad-bayes
```
こんな感じ． `executable mbpy-exe` の方にも同様に `monad-bayes` を書き足しましょう[^nec]．

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

```haskell:Lib.hs
module Lib (someFunc) where
import Control.Monad.Bayes.Simple
```

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

## Dice の用意

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

```bash:terminal
# mbplay に今いるとして
cp ../monad-bayes/models/Dice.hs src/
```
`.cabal` も編集して Dice っていうモジュールを export してるで，ということを書きます

```mbplay.cabal
library
  hs-source-dirs:      src
  exposed-modules:     Dice
```

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

```haskell: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`[^md] という型クラスの `d` に `Int` を包んで返しています． 1 から 6 まで確率が `1/6` ずつの categorical distribution として作られています． `dice` で `n` 回振れますが，これは素直な再帰．

この distribution からサンプルを取りたいのですが，元論文には `sample :: Samplable a => StdGen -> d a -> a` という分かりやすい関数があったのが色々と変更が入っています（この辺については次回以降にできれば書きたい）．結論だけいえば[`SamplerIO`](https://github.com/adscib/monad-bayes/blob/1a9fc60cf9f8529b336667e3878e7d853be291cc/src/Control/Monad/Bayes/Sampler.hs#L40-L46) を使うのが今回は正解で，初期化もして乱数を使ってとって来てくれる `sampleIO :: SamplerIO a -> IO a` が使えます． [`instance MonadDist SamplerIO`](https://github.com/adscib/monad-bayes/blob/1a9fc60cf9f8529b336667e3878e7d853be291cc/src/Control/Monad/Bayes/Sampler.hs#L82-L86) なので `die :: d Int` が直接いけまして

```Main.hs
main = sampleIO die >>= print
```

でさいころを振れます．


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

```haskell: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`の結果は

![out.png](https://qiita-image-store.s3.amazonaws.com/0/63937/d254c9e0-e7d7-1500-31f1-229d904fcaa9.png)

やりました．

---

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

[^hack]: [Readme](https://github.com/adscib/monad-bayes/commit/2f42955a14a6d5064631fa6f027d0e0333cddee1) にある程度安定したらアップするよって書いてありますがまだのようです
[^nec]: もちろんライブラリ側か `Main.hs` のどっちかでしか使わないならそちらの方だけでいいです．
[^md]: [`class (...) => MonadDist m where categorical :: [(a,CustomReal m)] -> m a ...`](https://github.com/adscib/monad-bayes/blob/1a9fc60cf9f8529b336667e3878e7d853be291cc/src/Control/Monad/Bayes/Simple.hs#L93-L105)

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a>.