{MCMCpack} イカサマサイコロを振ろう。

# おさらい
直近で二つの記事を書きました。
・[β分布](http://qiita.com/kilometer/items/3780f172bd314e934281)
・[ディリクレ分布](http://qiita.com/kilometer/items/6b6c9f451285c809167e#_reference-3eb75c16051be3a9bf06)

随時、参照してください。

# 確率分布からのサンプリング

これからサイコロを投げようと思うのですが、最初に確率分布からのサンプリングのお話。

例えば、

```math
X[x_1,...,x_{1000}]\sim N(μ=0,δ=1)
```
というのは、1000個のデータを平均0分散1の確率分布からサンプリングする事を示します。

至極当たり前の事を書きますが、これは個々のxが0に等しい事を保証するのではなく、
個々のxが大体0だよ、という意味であります。

で、このμを別の確率分布からサンプリングしてくる事を考えます。

```math
X[x_1,...,x_{1000}]\sim N(μ, 1) \\
μ\sim N(0, 1)
```

これは、Xはだいたい「だいたい0の分布からサンプルされた値」だ、と言っているわけです。

こうしたサンプリングを、概念図で書くとこんな感じになります。
(この書き方をgraphical modelと呼んだりします。)
<img width="300" alt="スクリーンショット 2016-09-16 11.04.08.png" src="https://qiita-image-store.s3.amazonaws.com/0/92401/5ef2aa42-8788-92c2-6bbf-7f9f9f425639.png">

実際のXの分布がどうなるかRで描いてみます。
とりあえずXを10回ぐらい生成してみます。

```R
set.seed(1)
N <- 10
myu <- rnorm(N, 0,1)

dat <- NULL
for(i in 1:N){
  dat <- rbind(dat, cbind(i, rnorm(1000, rnorm(1, myu[i], 1), 1)))
}
dat <- data.frame(dat)
```
<img width="111" alt="スクリーンショット 2016-09-16 10.42.17.png" src="https://qiita-image-store.s3.amazonaws.com/0/92401/a1b7e69b-2fab-5cba-7fdf-8d05e9fcdc8f.png">

で、描画

```R
library(ggplot2)

dat$i <- factor(dat$i)

ggplot(dat, aes(x=V2, color=i))+
  geom_density(size=1.5)
```
<img width=500 src=https://qiita-image-store.s3.amazonaws.com/0/92401/e878e11a-56bb-ea7d-6235-b195377e69fb.png>

与えた条件は、Xの平均の平均が0、Xの平均の分散が1、Xの分散が1の3つです。
まぁ概ね予想通りなのではないでしょうか。緑のぶっ飛んだヤツが出てくれて美味しい。

# サイコロを振ろう

正六面体の理想サイコロを振るとします。
結果は多項分布になりますよね。

```math
\frac{N!}{m_1!m_2!...m_K!}\prod_k μ_k^{m_k} \equiv Mult(m|N,μ)
```

μは、理想サイコロの場合、μ1~6が全て等しく1/6。
Rでとりあえず10回投げると、

```r
set.seed(1)
rmultinom(10, 1, c(1,1,1,1,1,1))
```

<img width=350, src=https://qiita-image-store.s3.amazonaws.com/0/92401/6e00ff73-9f64-eac9-1eec-4f8fad9c1de0.png>

`rmultinom`の2番目の引数は、同時に投げるサイコロの数、です。
3番目の引数が、各面の出る偏りを示す割合を入れますが、勝手に確率にスケールしてくれます。

とりあえず1000回投げて、結果を見てみると、

```r
set.seed(1)
dat <- rmultinom(1000, 1, c(1,1,1,1,1,1))
apply(dat, 1, sum)
```

![スクリーンショット 2016-09-16 11.26.47.png](https://qiita-image-store.s3.amazonaws.com/0/92401/d3d2706e-dab6-6141-ecdf-7e837c6823c3.png)

`1000/6=166.66..`なので、まぁ概ね均一ですね。

では、1が他の面よりも3倍も出やすいイカサマサイコロを振ってみましょう。

```r
set.seed(1)
dat <- rmultinom(1000, 1, c(3,1,1,1,1,1))
apply(dat, 1, sum)
```

![スクリーンショット 2016-09-16 11.31.32.png](https://qiita-image-store.s3.amazonaws.com/0/92401/04e527a5-199a-492e-8881-553c49dc601a.png)

これは、理想のイカサマサイコロですよね。
「1が他の面よりも3倍も出やすい」事が「真」だ、という意味です。
でも実際のイカサマサイコロは、この「3倍」が厳密ではなく、揺らぎを持っているわけです。

つまり、`c(3,1,1,1,1,1)`で与えていたμを、確率分布で与えたい。
μの分布？あ、それはディリクレ分布ですな。

```math
\frac{Γ(a_0)}{Γ(a_1)...Γ(a_K)}\prod_k μ_k^{a_k-1} \equiv Dir(μ|a)
```

「1が3倍出やすい**ような**」サイコロを100個ほど用意し、それぞれ1000回投げてみましょう。

グラフィカルモデルを書くと、こんな感じ。

<img width="300" alt="スクリーンショット 2016-09-16 11.50.54.png" src="https://qiita-image-store.s3.amazonaws.com/0/92401/81c70474-d7ac-7b7d-0876-bb47b17a0321.png">

ディリクレ分布の乱数は、`MCMCpack::rdirichlet`で発生させる。

```R
library("MCMCpack")
library("ggplot2")
library("tidyr")
set.seed(1)

N <- 100
a <- rdirichlet(N, c(3,1,1,1,1,1))

dat <- NULL
for(i in 1:N){
dat <- rbind(dat, apply(rmultinom(1000, 1, a[i,]), 1, sum))
}
dat <- data.frame(dat)

ggplot(gather(dat), aes(x=value, color=key))+
  geom_density(size=2)+
  xlim(0, 1000)
```

![スクリーンショット 2016-09-16 11.52.35.png](https://qiita-image-store.s3.amazonaws.com/0/92401/7017ce27-029c-6b40-02f8-06a6218d1432.png)

<img width="555" alt="スクリーンショット 2016-09-16 18.22.01.png" src="https://qiita-image-store.s3.amazonaws.com/0/92401/f04ada53-8b48-da7c-9aa5-bf0842750419.png">


サイコロNo.4なんかは、面1が95/1000回しか出ていない。
が、100個のサイコロで面kが出ていた回数の分布はグラフのようになり、
全体としては面1が極端に多く出るように偏っているのが分かる。

なんというか、いい感じに偏ったイカサマサイコロを選ぶ、というのが既にギャンブルですよね。