確率とモナドと確率的プログラミング

この記事は[ADVANCED BEGINNER](http://qiita.com/lotz/items/0d68c8440d1f362d0c32#advanced-beginner)から[COMPETENT](http://qiita.com/lotz/items/0d68c8440d1f362d0c32#competent)の方を対象読者として書かれています。

コインの裏表やサイコロの出目はよく[確率変数](https://www.wikiwand.com/ja/%E7%A2%BA%E7%8E%87%E5%A4%89%E6%95%B0#)によって表されます。確率変数が互いに依存しているようなモデルを記述する手法として[グラフィカルモデル](https://www.wikiwand.com/ja/%E3%82%B0%E3%83%A9%E3%83%95%E3%82%A3%E3%82%AB%E3%83%AB%E3%83%A2%E3%83%87%E3%83%AB)と云うものがあります。例えば、ある分布に従って表が出る確率が偏ったコインが選ばれた後、そのコインを投げて表裏が決まるような実験を考えた場合、コインの確率変数を $X$, コインの表裏の確率変数を $Y$ とすると、この系を記述するグラフィカルモデルは

![](https://i.gyazo.com/99aedd3a5e23f059190aba06e4746955.png)

このようになります。ところで $X$ は表が出る確率`Double`上の確率変数`RVar Double`で、 $Y$ は $X$ の結果に依存したコインの裏表`Bool`上の確率変数`Double -> RVar Bool`であると考えるとします。今コインがランダムに選ばれると言う構造を **忘れて** コインの表裏が出る確率変数`RVar Bool`を知りたいと思うと

```hs
RVar Double -> (Double -> RVar Bool) -> RVar Bool
```

このような関数が欲しくなるのではないでしょうか。お気付きの通りこれはモナドの`(>>=)`を思い出させるような形をしていて、実際`RVar`はモナドのインスタンスとなるように作ることができます。

確率分布はモナドだ
----------------
確率変数`RVar`を考える前に確率分布`Dist`を考えることにしましょう。さらに実装の簡便さからここでは離散確率分布を考えます。

```hs
type Prob = Double
newtype Dist a = Dist [(a, Prob)]
```

確率分布は確率変数のとる値から $[0,1]$ への関数ですが、`Dist`はそのグラフとなるように定義しています。`Prob`は単なる`Double`のエイリアスですが、値が $[0,1]$ を取ることを期待されるものとします。今、`f :: a -> b`なる関数があって`a`上の確率分布が与えられている時に、`b`上の確率分布`Dist b`を`Dist a`から生成した値に`f`を適用して得られた`b`の分布と考えると

```hs
instance Functor Dist where
  fmap f (Dist as) = Dist [(f a, p) | (a, p) <- as]
```

この様に`Dist`が`Functor`のインスタンスになることがわかります。

次に`f :: a -> b`が確率分布`Dist (a -> b)`から確率的に生成されていたとして同時に`a`上の確率分布`Dist a`が与えられていた時に、それぞれから生成された値を適用した`f a :: b`上の分布を考えると

```hs
instance Applicative Dist where
  pure x = Dist [(x, 1.0)]
  (Dist fs) <*> (Dist as) = Dist [(f a, p * q) | (f, p) <- fs, (a, q) <- as]
```

この様に`Applicative`のインスタンスも定義できると思います。ちなみにこれを使って

```hs
joint :: Dist a -> Dist b -> Dist (a, b)
joint da db = (,) <$> da <*> db
```

このような`a`と`b`の同時分布を求める関数を書くことができます。

最後に確率分布`Dist a`上の確率分布、すなわち`Dist (Dist a)`を考えてみましょう。`Dist (Dist a)`から値（確率分布）を生成し、さらに生成された確率分布から`a`の値を生成することを考えると

```hs
join :: Dist (Dist a) -> Dist a
join ass = Dist [(a, p * q) | (Dist as, p) <- ass, (a, q) <- as]
```

この様に`a`上の確率分布上の確率分布から`a`上の確率分布を得る関数を作ることができます。これを使えば

```hs
instance Monad Dist where
  m >>= k = join (fmap k m)
```

`Dist`を`Monad`のインスタンスにすることができました。

`Dist`を使った簡単な例を見てみましょう。サイコロの出目の確率分布`dice`を以下のように定義します。

```hs
dice :: Dist Int
dice = Dist $ map (\x -> (x, 1/6)) [1..6]
```

これを使って2つのサイコロの目の和を表す確率分布がどのようになるのかを計算してみます。その前に確率分布を表示する関数を作っておきましょう。

```hs
showDist :: (Show a, Ord a) => Dist a -> String
showDist (Dist as) =
  show .
    map (\as -> (fst (head as), sum (map snd as))) .
      groupBy (\a b -> fst a == fst b) .
        sort $ as
```

ここで作った`Dist`の実装では`a`の同じ値が中のリストに複数個現れてしまうのでうまくグループ化して確率を足して重複を無くした後に表示するようにしています。これらを使えば

```hs:Main.hs
main :: IO ()
main = do
  putStrLn $ showDist $ do
    x <- dice
    y <- dice
    pure (x + y)
```

```bash
$ runhaskell Main.hs
[(2,2.7777777777777776e-2),(3,5.555555555555555e-2),(4,8.333333333333333e-2),(5,0.1111111111111111),(6,0.1388888888888889),(7,0.16666666666666669),(8,0.1388888888888889),(9,0.1111111111111111),(10,8.333333333333333e-2),(11,5.555555555555555e-2),(12,2.7777777777777776e-2)]
```

![output.png](https://qiita-image-store.s3.amazonaws.com/0/17134/93fef080-807f-d078-aedf-38d9fa2a2da4.png)

この様に２つのサイコロの目の和を表す確率分布を作ることができました。

さて、確率変数の話に戻りましょう。`RVar`という型は実際[rvar](http://hackage.haskell.org/package/rvar)というHackageの中で定義されていて, `RVar`がモナドのインスタンスになることは確率分布がモナドのインスタンスになることを使って定義されています。[random-fu](http://hackage.haskell.org/package/random-fu)というHackageではよく知られた確率分布とそれに従う確率変数が定義されており、当然連続なものもたくさんあります。例えばこれを使って混合ガウス分布に従うような確率変数は

```hs
import Data.Random
import Data.Random.Distribution.Dirichlet
import Data.Random.Distribution.Exponential
import Data.Random.Distribution.Categorical

gaussianMixtureModel :: [Double]           -- 各正規分布の係数
                     -> [(Double, Double)] -- 各正規分布のパラメータ
                     -> RVar Double
gaussianMixtureModel weights params = do
  let gs = map (\(mu, sigma) -> normal mu sigma) params
  selected <- categorical $ zip weights gs
  selected

main :: IO ()
main = do
  results <- sample $ do
    weights <- dirichlet [1, 1]
    params <- sequence $ replicate (length weights) $ do
      mu <- normal 0 10
      sigma <- exponential 1
      pure (mu, sigma)
    sequence $ replicate 1000 $ gaussianMixtureModel weights params
  print results
```

![output.png](https://qiita-image-store.s3.amazonaws.com/0/17134/f4c09323-4115-34c6-116f-45621a4c06e2.png)

このように定義することができます。グラフィカルモデルがもつ[DAG](https://www.wikiwand.com/ja/%E6%9C%89%E5%90%91%E9%9D%9E%E5%B7%A1%E5%9B%9E%E3%82%B0%E3%83%A9%E3%83%95)の構造をモナドを使って綺麗に記述できているのが分かると思います。

### 参考
* [A Probabilistic Functional Programming Library for Haskell] (http://web.engr.oregonstate.edu/~erwig/pfp/)
* [プログラミングのための確率統計 in Haskell](http://note.golden-lucky.net/2010/12/1-2-3-4-5-6-16-16-16-16-16-16-246-135.html)
* [FUNCTIONAL PEARLS: Probabilistic Functional Programming in Haskell](https://web.engr.oregonstate.edu/~erwig/papers/PFP_JFP06.pdf)
* [Giry monad in nLab](https://ncatlab.org/nlab/show/Giry+monad)

確率的プログラミング言語
-----------------------

グラフィカルモデル（特にここではDAGとなるものを扱うので[ベイジアンネットワーク](https://www.wikiwand.com/ja/%E3%83%99%E3%82%A4%E3%82%B8%E3%82%A2%E3%83%B3%E3%83%8D%E3%83%83%E3%83%88%E3%83%AF%E3%83%BC%E3%82%AF)）を考えた時に実際に生成された値からモデルを構成する確率変数の確率分布を推定するのは重要な問題です。簡単な例で言えばあるデータが正規分布に従っていると仮定した時に実際に得られたデータからその平均と分散の値を求めるといったような問題です。それぞれの確率モデルごとに最適な推定アルゴリズムは違っていますが、問題に合わせて逐一実装するのは少々面倒でしょう。さらに複雑な確率モデルになっても最適な推定アルゴリズムが存在するとは限りません。そこでデータが得られた後の事後分布をサンプリングされた点列で近似することを考えます。この方法ならモデルがある程度複雑になっても同じ方法で近似的に事後分布を得ることができます。

このような確率モデルを構築してサンプリングによって事後確率分布を推定するというプロセスは色々な問題に共通して使うことができます。[確率的プログラミング](https://www.wikiwand.com/en/Probabilistic_programming_language)はプログラムによって確率モデルを記述し、実行としてサンプリングを行って推定結果が出力されるというアプローチでこれらの問題を統一的に扱います。確率的プログラミング言語はたくさん存在していて、古くからある[BUGS](http://www.openbugs.net/w/FrontPage)や[Stan](http://mc-stan.org/)は実務でも利用されています。これらはモデルを記述する専用の言語を用いて直感的にモデルを組み立てることができます。またLispを拡張した[Church](https://probmods.org/)という言語はモデルを組み立てるのにホスト言語の機能がそのまま使えると言うメリットが有ります。また[Anglican](http://www.robots.ox.ac.uk/~fwood/anglican/)や[PyMC3](http://pymc-devs.github.io/pymc3/)はそれぞれClojureとPythonのライブラリとして提供されており既存の開発環境をそのまま使えるというメリットが有ります。Haskellを使って作られた言語としては[Hakaru](http://hakaru-dev.github.io/)や[BAli-Phy](http://www.bali-phy.org/)といったものがあります。

確率的プログラミング言語は新しい言語を作ることで確率変数を第一級オブジェクトとして使えるという利点があります。その結果、宣言が簡潔になったりifなどの文法と組み合わせて使うことができるようになります。ところでHaskellでは確率変数はモナドのインスタンスになることを見てきたと思います。この性質を使えば元の言語に手を加えることなくあたかも確率変数が第一級オブジェクトかのようにHaskellでも確率的プログラミングを扱うことができるようになるのではないでしょうか。

Haskellで実装してみる
-------------------
前に作った`Dist`では確率モデルを記述することができますが得られたデータから確率分布を推定する事はできませんでした。なのでここでは得られたデータを使って事後分布をサンプリングする機能を作っていきたいと思います。なおここでの話は [Practical Probabilistic Programming with Monads](http://mlg.eng.cam.ac.uk/pub/pdf/SciGhaGor15.pdf) を大いに参考にして書かれています。

今回作る`Dist`は

* `random-fu`で定義されている確率変数を埋め込めるように
* データからサンプリングが出来るように

してみたいと思います。

```hs
{-# LANGUAGE GADTs #-}

type Prob = Double

data Dist a where
  Pure :: a -> Dist a
  Bind :: Dist b -> (b -> Dist a) -> Dist a
  Primitive :: RVar a -> Dist a
  Conditional :: (a -> Prob) -> Dist a -> Dist a

condition :: (a -> Prob) -> Dist a -> Dist a
condition = Conditional
```

前の`Dist`とはぜんぜん違うものになってしまいました。それもそのはず`Dist`はもはや直接確率分布を表すものではなくGADTｓを使った単なるDSLになってしまったのです。気持ちとしては`Pure`と`Bind`は`Monad`のインスタンスにするために、`Primitive`は`random-fu`の`RVar`を埋め込むために、`Conditional`はデータから得られた条件を反映するために用意しています。これらのDSLに意味を持たせるにはインタプリタを実装する必要がありますが、その前に一旦`Monad`までのインスタンスを定義してしまいましょう。

```hs
instance Functor Dist where
  fmap f d = d >>= (pure . f)

instance Applicative Dist where
  pure = Pure
  f <*> a = a >>= \x -> fmap ($x) f

instance Monad Dist where
  (>>=) = Bind
```

`Pure`と`Bind`はモナドになることを意識して用意しているので`Monad`のインスタンスはすんなり作ることができます。`Monad`であれば`Applicative`と`Functor`のインスタンスは自動的に出来上がりますね。
さてインタプリタですが、

```hs
instance Sampleable Dist IO a where
  sampleFrom _ (Pure x) = pure x
  sampleFrom s (Bind d f) = do
    x <- sampleFrom s d
    sampleFrom s (f x)
  sampleFrom s (Primitive d) = sampleFrom s d
  sampleFrom _ (Conditional _ _) = undefined
```

この様に`random-fu`の`Sampleable`の実装に押し付ける形で作ります。これで`Dist a`から`IO a`に変換できるようになるはずです。ここで`Conditional`の場合の返り値が`undefined`になってて危険そうですが、`sampleFrom`は`Conditional`を全て取り除かれた`Dist`に対して適用されるという前提で話を進めていきます。

`Conditional`が全て取り除かれた状態とはどういうものかというと与えられた条件に対してベイズ更新が全て終わった後の事後確率分布です。事後確率分布を求めるためにいよいよサンプリングの関数を作っていきましょう。

```hs
import qualified Data.Random.Distribution.Bernoulli as Random

bernoulli a = Primitive $ Random.bernoulli a

prior :: Dist a -> Dist (a, Prob)
prior (Conditional c d) = do
  (x,s) <- prior d
  return (x, s * c x)
prior (Bind d f) = do
  (x,s) <- prior d
  y <- f x
  return (y,s)
prior d = do
  x <- d
  return (x,1)

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
```

`mh`は[メトロポリスヘイスティング法](https://www.wikiwand.com/ja/%E3%83%A1%E3%83%88%E3%83%AD%E3%83%9D%E3%83%AA%E3%82%B9%E3%83%BB%E3%83%98%E3%82%A4%E3%82%B9%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0%E3%82%B9%E6%B3%95)を用いて事後分布を求める関数です。これは設定値を受け取って`Conditional`を含む事前分布から事後分布からサンプリングした点列上の確率分布を求める関数になっています。返り値の`Dist`には`Conditional`は含まれませんので安心して結果を取り出すことができます。 `prior`は`mh`の中で提案分布として使用するために尤度を引き回しながら事前分布を求める関数です。

それではこれらを使って簡単なベイズ線形回帰を試してみましょう。

```hs
import qualified Data.Random as Random

normal a b = Primitive $ Random.normal a b

type Point = (Double, Double)
type Param = (Double, Double)

linear :: Dist Param
linear = do
  a <- normal 0 5
  b <- normal 0 5
  return (a,b)

point :: Point -> Dist Param -> Dist Param
point (x,y) = condition (\(a,b) -> pdf (Random.Normal (a*x + b) 1) y)

points :: [Point] -> Dist Param -> Dist Param
points ps d = foldl (flip point) d ps

trainingData :: IO [(Double, Double)]
trainingData =
  let points = [-5..5]
      equation = \x -> 2 * x + 1
      addNoise x = sample $ Random.Normal x 1
   in sequence . map (\x -> (,) x <$> addNoise (equation x)) $ points

main :: IO ()
main = do
  td <- trainingData

  putStrLn "Sampling..."
  xs <- sample $ mh 100 100000 (points td linear)

  putStrLn "Drawing..."
  r2AxisMain $ myaxis $ nub xs
    where
      myaxis :: [(Double, Double)] -> Axis B V2 Double
      myaxis mydata1 = r2Axis &~ do
        scatterPlot mydata1 $ key "data 1"
```

学習データとして $y = 2x + 1$ に標準正規分布に従うノイズを加えたものを使っています。グラフのプロットには[plots](http://hackage.haskell.org/package/plots)というライブラリを使っています。

```bash
$ runhaskell Main.hs
Sampling...
Results:
a: 1.9233401839331328
b: 0.9046290119555563
Drawing...
```

プロットした図はこんな感じ。

![output.png](https://qiita-image-store.s3.amazonaws.com/0/17134/11c4779f-0769-8932-0708-f852042897c9.png)

まとめ
------
確率分布はモナドになり、それを使ってモデルの記述やパラメータの推定が直感的で簡潔に書けることを紹介してきました。本当は確率的プログラミングの型理論の話から量子コンピュータの話まで書きたかったけど知識も時間も足りずに挫折…。

この記事は[第37回Haskellもくもく会@朝日ネット](https://haskellmokumoku.connpass.com/event/45821/)のもくもく時間を利用して書かれました。