β分布

Qiitaに書いてどーすんだ説はありますが。
ディリクレ分布の理解を進めるにはβ分布をガツンと理解しておく事に限る、と思います。
その足跡記事です。

結果的に、[黄色い本](https://www.amazon.co.jp/dp/4621061224)の2.1をネチネチと説明した、という感じの文章になってしまいました。

# 1回のコイントス
ここに、{表1, 裏0}がある確率分布に従って出るコインがある。
次にコイントスをした際に、表裏どちらが出るかで**私の人生**が決まる。
やばいっすね。
計算しましょう。


コインのどちらが出るか、という確率は、何か`μ`に支配されているはずです。

```math
\begin{eqnarray}
p(x=1|μ)&=&μ,&\ where\ 0≤μ≤1 \\
p(x=0|μ)&=&1-μ&
\end{eqnarray}
```

この`μ`のもとで`x`がどんな値をとるかという確率は、

```math
μ^x(1-μ)^{1-x} \equiv Bern(x|μ)
```

これは、**独立した1回の試行で**コインがとる値の確率分布です。
（これをベルヌーイ分布 Bernoulli distributionと呼びます。）

例えば`μ=0.5`の時、`x`は`Bern(x|μ)`に従うので、

```math
p(x = 1|μ = 0.5) = 0.5^1\times(1-0.5)^{1-1} =0.5\\
p(x = 0|μ = 0.5) = 0.5^0\times(1-0.5)^{1-0} =0.5
```

と計算できます。

# 最尤推定
ベルヌーイ分布を支配している`μ`。
我々は、その正体を明らかにせねばなるまい。
シスの暗黒卿は、パルパティーン議員だったのだ。

手元のデータ

```math
D \{ x_1, x_2, ...\ x_i\ ... , x_N \}
```
があるとして、`μ`を最尤推定してみる。

`D`を**独立**なベルヌーイ試行の集合とすると、

尤度は、`μ`のもとでそれぞれのデータが得られる確率を単純に掛ければよいから、

```math
p(D|μ)=\prod_{i=1}^N p(x_i|μ)=\prod_{i=1}^N μ^{x_i}(1-μ)^{1-x_i}
```

尤度の最大化は対数尤度の最大化と等価なので、

```math
\begin{eqnarray}
ln\ p(D|μ)
&=& ln\ \prod_{i=1}^N p(x_i|μ) \\
&=& \sum_{i=1}^N ln\ p(x_i|μ) \\
&=& \sum_{i=1}^N \big( x_i ln\ μ+(1-x_i)ln(1-μ) \big)
\end{eqnarray}
```

これが最大となるのは、微分=0の時なので、

```math
\begin{eqnarray}
\frac{d}{dμ} \Big( ln\ p(D|μ) \Big) &=& 
\frac{d}{dμ}\sum_{i=1}^N \big( x_i ln\ μ+(1-x_i)ln(1-μ) \big) \\
0&=&\frac{1}{μ_{ML}}\sum x_i\ - \frac{1}{1-μ_{ML}} \sum (1-x_i)  \\
μ_{ML}&=& \frac{\sum x_i}{\sum x_i+\sum(1-x_i)} \\
&=& \frac{1}{N}\sum_{i=1}^N x_i
\end{eqnarray}
```

こ、これは…**平均**ですね(標本平均)。

つ、つまり、これまでの`N`回の試行`D`において、

```math
p(x=1|μ)=μ=\bar{x}
```

これまで100回コインを投げて40回表が出たら、これまでの実績から考えた最も尤もらしいコインが表になる確率は0.4。

うーむ、そらそうだ。


# 事前分布

ただですね、これまで3回の実績しかないとします。
それが全て表だった場合、平均は1なので、以降永遠に表が出続けるというのが最も尤もらしいということになります。

これは3回で**全て**が表されているという前提が怪しい。
その限定された実績で**過**学習が起きているという事です。

**今更**ですが、コインは世界にこれ1枚では無いわけです。
もちろん、完全に同一なコインは1枚だけですが、大雑把に見れば、似たようなコインは山ほどあり、Homo族の猿の一種はすぐコイントスをする性質がありますよね。で、概ね50%50%ぐらいの割合なんじゃないっすかね？

これを**前提**として、今起きた「3回連続表」という現象を評価すうると、やや表に偏っているかな、ぐらいの評価になるはずです。これは「以降、永遠に表が出続ける」よりもずいぶんと穏やかで妥当な推測ですよね。

こうした**前提**の事を、「事前分布」と言います。

という事で、**人生が掛かってます**からね。
大先輩[ケダリオン](https://ja.wikipedia.org/wiki/%E3%82%B1%E3%83%BC%E3%83%80%E3%83%AA%E3%82%AA%E3%83%BC%E3%83%B3)の如く、事前分布の肩の上に乗りましょう。

# 求めたいモノ

我々は既に**決意**を固めました。
もはや次の1回の(そしてこれまで行われてきたあらゆる)コイントスを、**独立な**ベルヌーイ試行とは**考えません**。

次の1回のコイントスで表がでる確率は、**これまでのデータ**に縛られる。


```math
p(x=1|D,μ)
```
`0 ≤ μ ≤ 1`より、

```math
\begin{eqnarray}
p(x=1|D,μ) &=& \int_{0}^{1} p(x=1|μ)\ p(μ|D)dμ \\
&=& \int_{0}^{1} μ\ p(μ|D)dμ
\end{eqnarray}
```
右辺の形は、「データに基づいてμが発生する確率」と「μの実現値」の積の総和ですわな。
これは、**期待値**と呼ばれるやつです。

```math
p(x=1|D,μ)=\mathbb{E}[μ|D]
```

データが与える`μ`の期待値を求めればよろしい。
データとは、「N回の試行の結果、k回の{1}が得られた」という情報です。

（あとでここに戻ります）

# N回試行の尤度
今、我々は、[B系](https://ja.wikipedia.org/wiki/%E3%83%99%E3%82%A4%E3%82%BA%E7%A2%BA%E7%8E%87)のモノの考え方を採用しています。

```math
事後分布 \propto 尤度関数 \times 事前分布
```

推定したいのは`μ`なので、事前分布・事後分布は`p(μ|D)`という形になるはずです。

従って、尤度関数は、`p(D|μ)`の形になる。
この`D`には、「N回投げたらk回の{1}が得られた」という情報が入るはずです。

`k`の分布を考えると、`N`は変数になるので、

```math
p(k|N,μ)= {}_NC_k\ μ^k(1-μ)^{N-k} \equiv Bin(k|N,μ)
```
これを二項分布(Binomial distribution)とモウします。

ここで確率の左側に来ているのが、コインを何度か投げたうちの{1}が出た**回数**、という事に注意してください。


さて、`μ`に依存する要素だけで書き直すと、

```math
\begin{eqnarray}
p(k,l|μ)&=& {}_NC_k\ μ^k(1-μ)^{l},\ where\ N=k+l \\
& \propto & μ^k(1-μ)^{l}
\end{eqnarray}
```

# 事前分布、どーする

再掲。

```math
\begin{eqnarray}
事後分布 &\propto& 尤度関数 &\times 事前分布 \\
& \propto & μ^k(1-μ)^{l} &\times 事前分布
\end{eqnarray}
```

もう一度回り道をしますが、ベイズ的にやりたい事は、事前分布を元に事後分布を計算し、それを事前分布として事後分布を計算し、...という情報のアップデートを積み重ねてより正確な`μ`の事後分布を知りたい、という事です(逐次計算)。

そのためには、正の実数`a`,`b`を用いて、

```math
事前分布 \propto μ^{a}(1-μ)^{b}
```

という格好になっているとスゴく便利ですよね。


```math
\begin{eqnarray}
事後分布 &\propto& 尤度関数 \times 事前分布 \\
& \propto & μ^k(1-μ)^{l} \times μ^{a}(1-μ)^{b}\\
& \propto & μ^{k+a}(1-μ)^{l+b}
\end{eqnarray}
```

となって、これはすぐに次の事前分布として利用できます。

正規化項を加えて、

```math
p(μ|a,b)=\frac{\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}μ^{a}(1-μ)^{b} 
```

(20160819. [7of9](http://qiita.com/7of9)様よりご指摘をいただき修正しました。)


としますが、ちょっとお行儀を良く書き直してβ分布(beta distribution)と呼びます。

```math
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \equiv Beta(μ|a,b) \\
```
ただし、

```math
\Gamma(x)\equiv \int_0^\infty u^{x-1}e^{-μ}du
```


なんでワザワザaをa-1に書き直したのかというと、期待値をすっきり書けるからです。

```math
\mathbb{E}[Beta(μ|a,b)]= \frac{a}{a+b}
```

Γ関数や期待値の詳細をAppendixにすっ飛ばして先に進みます。

事後分布は、

```math
p(μ|k,l,a,b) \propto  μ^{k+a-1}(1-μ)^{l+b-1}
```

正規化項をしっかり計算すると、

```math
p(μ|k,l,a,b) = \frac{\Gamma(k+a+l+b)}{\Gamma(k+a)\Gamma(l+b)} μ^{k+a-1}(1-μ)^{l+b-1}
```

これ自体がBeta分布になっています。（共役事前と言います）
ので、

```math
\mathbb{E}[Beta(μ|k,l,a,b)]= \frac{k+a}{k+a+l+b}
```

（あとでここに）を再掲

```math
p(x=1|D,μ)=\mathbb{E}[μ|D]
```

試合終了。


```math
p(x=1|D)=\frac{k+a}{k+a+l+b}
```

# Beta分布とは何者だったのか。

[祭りのあと](http://www.kasi-time.com/item-5917.html)に残った感慨は、結局、一体全体、Beta分布とは何か、という話です。

**コインを1回振って**出た結果の分布が、ベルヌーイ分布。
**コインをN回振って**表が出る回数の分布が、二項分布。

そして、**コインそのものの性質を決めるのが、Beta分布**です。
（これが**だけ**が言いたかった）

変数a, bを様々に取る事で、完璧なコインからイカサマコインまで、あらゆるコインの`μ`の事前分布＆事後分布を与える事ができる(ようにデザインされている)。

こうした、実際の試行に先立って与えられたコインの性質を規定している変数a, bの事を、ハイパーパラメータと呼んだりします。

# Appendix
## 1. Γ関数

```math
\Gamma(x)\equiv \int_0^\infty u^{x-1}e^{-μ}du
```

β分布の正規化項で暗闇を破る雷鳴の如く登場したこの関数。

復習：部分積分

```math
\int f(x)g'(x) dx = f(x)g(x)-\int f'(x)g(x)dx
```

ですよな。

```math
f(x)=u^{x-1}\\
g(x)=-e^{-u}
```

として、見くらべて頂きながら、

```math
\begin{eqnarray}
\Gamma(x)&=& \int_0^\infty u^{x-1}e^{-μ}du \\
&=& \Big( -u^{x-1}e^{-u} \Big)_0^\infty + (x-1)\int_0^\infty u^{x-2}e^{-u}du \\
&=& (x-1)\int_0^\infty u^{x-2}e^{-u}du \\
&=& (x-1)\Gamma(x-1)
\end{eqnarray}
```

という事で、これの正体は一般化された階乗関数です。

Nが自然数であれば、

```math
\Gamma(N+1) = N!
```

## 2. β分布の正規性

本文より、


```math
事前分布 \propto μ^{a}(1-μ)^{b}
```

となって、いるとスゴく便利で、

正規化項を加えて、行儀良く書き直してβ分布(beta distribution)と呼びます。

```math
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \equiv Beta(μ|a,b) \\
```

え〜と、ですね。正規化項がこの形なのは、立派な理由がありまして。
立派すぎて、書くのを躊躇われます。

```math
\int_0^1μ^{a-1}(1-μ)^{b-1}dμ\ \equiv\ B(a,b)
```

と置きます。β関数と呼びます。

```math
μ=sin^2θ
```
と置きます。

```math
\begin{eqnarray}
B(a,b)&=&\int_0^1μ^{a-1}(1-μ)^{b-1}dμ\\
&=&\int_0^{\frac{\pi}{2}} \big(sin^2θ \big)^{a-1} \big( 1-sin^2θ \big)^{b-1}\ 2sinθcosθdθ \\
&=&2\int_0^{\frac{\pi}{2}} sin^{2a-1}θ\ cos^{2b-1}θ\ dθ \\
\end{eqnarray}
```
となりますね。

一方、Γ関数において`u=s^2`と置いて、

```math
\begin{eqnarray}
\Gamma(x) &=& \int_0^\infty u^{x-1}e^{-μ}du
&=&2\int_0^\infty s^{2x-2}e^{-2s^2}ds
\end{eqnarray}
```

さぁて、うまくいくかな…!?

```math
\begin{eqnarray}
Γ(a)Γ(b)&=&
\Big( 2\int_0^\infty  x^{2a-1}e^{-x^2}dx \Big)
\Big( 2\int_0^\infty  y^{2b-1}e^{-y^2}dy \Big) \\
&=& 4\int_0^\infty\int_0^\infty x^{2a-1} y^{2b-1}e^{-(x^2+y^2)}dxdy \\
&=& 4\int_0^\infty\int_0^{\frac{\pi}{2}}
 (r\ cosθ)^{2a-1} (r\ sinθ)^{2b-1}e^{-r^2}dθdr \\
&=& 4\int_0^\infty\int_0^{\frac{\pi}{2}}
 r^{2(a+b)-2}e^{-r^2}cos^{2a-1}θ\ sin^{2b-1}θ\ dθdr \\
&=&
\Big( 2\int_0^\infty  r^{2(a+b)-2}e^{-r^2} dr \Big)
\Big( 2\int_0^{\frac{\pi}{2}}  cos^{2a-1}θ\ sin^{2b-1}θ\ dθ \Big) \\
&=&Γ(a+b)B(a,b) \\
\end{eqnarray}
```

ぶふー。極座標変換ですね。`x=r cosθ`, `y=r sinθ`
xもyも0~∞つまり正の値だから、θの積分範囲は0~pi/2になる。
で、要するに

```math
B(a,b)=\frac{Γ(a)Γ(b)}{Γ(a+b)}
```

つまり、


```math
\begin{eqnarray}
Beta(μ|a,b)&=&\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \\
&=&\frac{1}{B(a,b)}μ^{a-1}(1-μ)^{b-1}
\end{eqnarray}
```

と書けますね。

そして、

```math
\begin{eqnarray}
\int_0^1 Beta(μ|a,b)dμ 
&=&\int_0^1  \frac{1}{B(a,b)}μ^{a-1}(1-μ)^{b-1}dμ \\
&=&\frac{1}{B(a,b)}\int_0^1 μ^{a-1}(1-μ)^{b-1}dμ \\
&=&\frac{B(a,b)}{B(a,b)}=1
\end{eqnarray}
```

おお、正規化されていました。

## 3. β分布の期待値

息切れしてきました。
Texのトレーニングはもう…。

まず、`B(a-1,b)`を考えます。

```math
\begin{eqnarray}
B(a-1,b)&=&\frac{Γ(a-1)Γ(b)}{Γ(a+b-1)} \\
&=& \frac{Γ(a)}{a-1} \frac{Γ(b)}{Γ(a+b-1)} \\
&=& \frac{a+b-1}{a-1} \frac{Γ(a)Γ(b)}{Γ(a+b)} \\
&=& \frac{a+b-1}{a-1} B(a,b) \\
\end{eqnarray}
```


いざ、期待値は、実現値と確率の掛け算の総和っと叫んで突撃。

```math
\begin{eqnarray}
\mathbb{E}[Beta(μ|a,b)] &=& \int_0^1 Beta(μ|a,b)\ μ\ dμ  \\
&=&\int_0^1 \frac{μ^{a-1}(1-μ)^{b-1}}{B(a,b)}\ μ\ dμ \\
&=&\int_0^1 \frac{μ^{a}(1-μ)^{b-1}}{B(a,b)}\ dμ \\
\end{eqnarray}
```

ここで、`a=α-1`と置くと、

```math
\begin{eqnarray}
\mathbb{E}[Beta(μ|a,b)]
&=&\int_0^1 \frac{μ^{α-1}(1-μ)^{b-1}}{B(α-1,b)}\ dμ \\
&=&\frac{α-1}{α+b-1}
\int_0^1 \frac{μ^{α-1}(1-μ)^{b-1}}{B(α,b)}\ dμ \\
\end{eqnarray}
```

正規化されているから、積分の中身＝１。
`α=a+1`に戻して、

```math
\mathbb{E}[Beta(μ|a,b)]=\frac{a}{a+b}
```


ぶ、分散は…気が向いたら追加します。(永久に向かない可能性大)

# 参考文献
書籍
[黄色い本](https://www.amazon.co.jp/dp/4621061224)

webサイト
http://bio-info.biz/statistics/distribution_beta_distribution.html
https://lecture.ecc.u-tokyo.ac.jp/~nkiyono/2006/miya-gamma.pdf