TL;DR

共役性を持つ事前分布を使うことで、事後分布の計算が非常に楽になる。
二項分布、多項分布においては事前分布に共役性のあるベータ分布、ディリクレ分布を使うことで
コイン、サイコロの観測を行う前（事前）に
コインの表、サイコロが$k$の目である回数がいくらか得られているように
補正をすることに相当する

用語の整理

密度推定 (density estimation)
観測値の集合$\{x_1, \cdots, x_N\}$が与えられたときに確率変数の従う確率分布を求めること。
パラメトリックなモデル
少数のパラメータで確率分布（関数）が一意に決まるようなモデルのこと。1章で出てきたガウス分布は、平均と分散という２つのパラメータが決まると一意に分布が決まるので、パラメトリックである。
共役性 (conjugacy)
尤度と事前分布を掛けて事後分布を求める際、事後分布が事前分布と同じような関数の形になる
性質のこと。これにより、数学的に非常に簡単に計算ができるようになる。

2.1 二値変数

ベルヌーイ分布

0か1の二値しかとらない確率変数を考える。$x=1$が表、$x=0$が裏と思うと、この確率変数は
コイン投げの結果を表している。
$x=1$となる確率をパラメータ$\mu$で表す（$\mu$は確率なのでもちろん$0 \leq \mu \leq 1）$。
式で書くと
$$
p (x=1 | \mu ) = \mu \tag{2.1}
$$
コインが表 ($x=1$) の確率が$\mu$なんだから、当然$p(x=0 | \mu) = 1 - \mu$。
$x=0$のときと$x=1$のときの両方を1つの式にまとめて表現したのが、ベルヌーイ分布と呼ばれる。
$$
\text{Bern}(x | \mu) = \mu^x (1 - \mu)^{1-x} \tag{2.2}
$$
実際に$x=0$, $x=1$の両方を代入してみるとすぐ分かる。
期待値と分散を求めてみよう。期待値は

\begin{align}

\mathbb{E}[x] &= (x=0の確率)×(x=0のときのxの値) + (x=1の確率)×(x=1のときのxの値) \\
&= (1 - \mu) \cdot 0 + \mu \cdot 1 \\
&= \mu

\end{align}

分散の計算には式(1.39)を使ってみよう。まず$\mathbb{E}[x^2]$は

\begin{align}

\mathbb{E}[x^2] &= (x=0の確率)×(x=0のときのx^2の値) + (x=1の確率)×(x=1のときのx^2の値) \\
&= (1 - \mu) \cdot 0 + \mu \cdot 1
&= \mu
\end{align}

なので、

\begin{align}

\text{var}[x] &= \mathbb{E}[x^2] - \mathbb{E}[x]^2 \\
&= \mu - \mu^2 \\
&= \mu (1 - \mu)
\end{align}

観測値$ D = \{x_1, \cdots, x_N\}$が得られたときの尤度関数¹は
$$
p(D | \mu) = \prod_{n=1}^N \mu^{x_n} (1 - \mu)^{ (1 - x_n) }
$$
最尤推定で$\mu$を求めよう。前章同様、対数の尤度をとってから$\mu$で偏微分=0を解く。

\begin{align}

\ln p = \sum_{n=1}^N \left\{ x_n \ln \mu + (1 - x_n) \ln (1 - \mu) \right\} \\

\frac{ \partial \ln p }{ \partial \mu } = \frac{1}{\mu} \sum x_n - 
\frac{1}{1 - \mu} \sum (1 - x_n) = 0 \\
両辺に\mu(1 - \mu)をかけて \\
(1 - \mu) \sum x_n - \mu \sum (1 - x_n) = 0 \\
\sum x_n - \mu \sum x_n - \mu N + \mu \sum x_n = 0 \\
\mu_{ML} = \frac{1}{N} \sum x_n \tag{2.7}
\end{align}

表が出る確率$\mu$の推定値はコインを投げた回数分の表が出た回数² というそのまんまの結果が導出できた。

二項分布 (binomial distribution)

ベルヌーイ分布が1回のコイン投げで表が出る確率を表していたのに対し、
二項分布はN回のコイン投げでm回表が出る確率を表す。
表が出る確率が$\mu$なので
表が$m$回、裏が$(N-m)$回でる確率は$\mu^m (1-\mu)^{(N-m)}$に比例することは
イメージがつくだろう。このままだと$m$が0の確率から$m$が$N$の確率までを
足し合わせると1にならない（N=2の場合で計算してみると簡単に分かる）。
これでは確率分布の定義を満たしていない（足し合わせると1にならないとダメ）ので
係数をかける。その係数は$N$個のうち$m$個の組み合わせの数で与えられるので、
式で書くと
$$
\text{Bin}(m| N, \mu) = \frac{N!}{(N-m)!m!} \mu^m (1-\mu)^{N-m}
\tag{2.9}
$$

2.1.1 ベータ分布

二項分布のパラメータ$\mu$を（最尤推定でなく）MAP推定で解くには
先に述べたように$\mu$の事前確率分布が必要になる。しかし、適当な事前分布に
設定してしまうと式が複雑になり、解くのが大変になる。
そこで、式を簡単にするために、尤度と事前確率の積に比例する事後確率が
尤度と同じ関数形式（これについてはこの後すぐ説明）になるような事前確率を設定したい。このような、事後分布と事前分布が同じ関数形式になる性質を
共役性 (conjugacy) という。
ここでいう関数形式は式(2.9)のように$\mu$のべき乗と$(1-\mu)$のべき乗で
表される形とのこと。尤度がべき乗の積で書けるので、事前分布もべき乗の積で
書けるものを使えば、事後分布もべき乗の積で書けるものになる。
このべき乗の積で書ける分布の一つがベータ分布であり、二項分布の事前分布として
よく用いられる。
$$
\text{Beta}( \mu | a,b) =\frac{ \Gamma(a+b) }{ \Gamma(a)\Gamma(b) }
\mu^{a-1} (1 - \mu)^{b-1} \tag{2.13}
$$
ベータ分布にガンマ関数からなる複雑な係数がついているが、こうしておくとテキストの図のように
$\mu$がどの値をとりそうかという分布の形をいろいろ変えることができ、
なおかつ足し合わせて1になるという確率分布の定義を満たすためにある、とぐらいに
理解しておけばいい
ベータ関数を事前分布としたとき、事後分布は式(2.9)と(2.13)の積なので$l$が
$l = N - m $、すなわちコイン投げの裏の回数とすると
$$
p( \mu | m,l,a,b) \propto \mu^{m+a-1} (1-\mu)^{l+b-1} \tag{2.17}
$$
となり、事前分布であるベータ関数$\mu^{a-1} (1 - \mu)^{b-1}$のべき乗の部分に
表の回数$m$と裏の回数$l$を足しただけの形となっており、計算がすごく簡単になっている。
事後分布についてもう少し詳しく見てみよう。
観測データ集合$D$が与えられ、次の観測が$x=1$ (表) となる確率$p(x=1 | D)$を予測しよう。
- 式変形にあたり、$p(x=1 | D)$に$x=1$になる確率$\mu$を入れ込む。
  $p(x=1, \mu | D)$と書き、$\mu$について確率のsum rule (式1.10) を使って$\mu$が消えるようにする。
  　これで、
  $$
  p(x=1 | D) = \int_0^1 p(x=1, \mu | D) d\mu
  $$
  と書ける。
- 次にproduct rule (式1.11)を使い、
  $$
  \int_0^1 p(x=1, \mu | D) d\mu = \int_0^1 p(x=1 | \mu) p(\mu | D) d\mu
  $$
  と変形する。$p(x=1 | \mu)$は表の確率$\mu$が与えられたときに表 ($x=1$) が出る確率の意味であり、
  $p(x=1 | \mu)$は$\mu$そのもの。なので、
  $$
  \int_0^1 p(x=1 | \mu) p(\mu | D) d\mu =
  \int_0^1 \mu p(\mu | D) = \mathbb{E}[ \mu | D ] \tag{2.19}
  $$
  　となる。
$p(x=1 | D)$は$D$が与えられたとき（事後分布）の$\mu$の期待値$\mathbb{E}[ \mu | D ]$であるということになった。事後分布は式2.17の形のベータ分布になることとベータ分布の期待値の計算式2.15より
$$
p(x=1 | D) = \frac{m+a}{m+a+l+b} \tag{2.20}
$$
となる。事前分布を考えず、（最尤推定で）コインが表になる確率を推定すれば、
コインを振った回数分の表が出た回数、すなわち$\frac{m}{m+l}$となる。
これと式2.20を比較すると、実際の表と裏の観測値$m, l$に仮想の観測値$a, b$を合わせている
ことが分かる。データを観測する前（事前）の状態から、$a, b$を与えることで
観測後（事後）の表になる確率の推定値に補正をかけているのだ。

2.2 多値変数

多値でのxの分布

0か1かの二値しか扱えなかった二値変数を多値が扱えるように拡張したものをこれから扱う。
例としては目が1から6まであるサイコロを考えればよい。
多値変数を扱う上で便利な表現である1-of-K符号化法（one-hotベクトルともいう）を見てみよう。
$K$次元ベクトル$\mathbf{x}$の$k$番目の値が出たときは$x_k$を1とし、その他は0とするベクトルで
表現する。
$$
\mathbf{x} = (0,0,1,0,0,0)^T \tag{2.25}
$$
サイコロを例にすれば、$K$は6となり、目が3の場合は上式(2.25)のように
ベクトルの3番目の要素を1とし、他を0として表現する。
1の目が出る確率から6の目が出る確率までを表した$K$次元のベクトルを$\mathbf{\mu}$とおく。
すると、$x_l=1$となる確率を１つの式にまとめて次のように書ける。
$$
p( \mathbf{x} | \mathbf{\mu} ) = \prod_{k=1}^K \mu_k^{x_k}
$$
$l=k$のとき確率は$\mu_l^1 = \mu_l$となり、$l \neq k$のときは$\mu_l^0 = 1$となるので
これらを掛け合わせると結局$\mu_l$となるため。
N個の観測値$\mathbf{x_1}, \cdots, \mathbf{x_N}$が与えられた際の尤度関数は
$$
p(D | \mathbf{\mu} ) = \prod_{n=1}^N \prod_{k=1}^K \mu_k^{x_{nk}}
= \prod_{k=1}^K \mu_k^{m_k} \tag{2.29}
$$
ここで${x_{nk}}$は$n$番目の観測値$\mathbf{x_n}$の$k$番目の要素、
$m_k$は$m_k = \sum_n x_{nk}$であり、N回サイコロを振った中で目が$k$となった合計回数を表す。
尤度関数が書けたので、最尤推定でパラメータ$\mathbf{\mu}$を推定しよう。
$\mathbf{\mu}$は確率であるため、$\sum_k \mu_k = 1$を満たさないといけない。
そこで、最尤推定の際には単純に尤度関数を最大化するだけではなく、
$\sum_k \mu_k = 1$という制約が満たされる範囲で尤度が最大となる$\mathbf{\mu}$を求めたい。
このように制約の下、関数を最大化するにはラグランジュの未定乗数法と呼ばれる手法を使う。
　　- ラグランジュの未定乗数法では(最大化したい式)+$\lambda$ (制約式) = 0と式をおき、
　　この式をパラメータについて偏微分した式、$\lambda$について偏微分した式の
　　連立方程式を解けばよい。
- 式2.29の対数、つまり対数尤度は
  $$
  \sum_{k=1}^K m_k \ln \mu_k
  $$
  と書けるので、ラグランジュの未定乗数法で解く式はつまり、
  $$
  \sum_{k=1}^K m_k \ln \mu_k + \lambda \left( \sum_{k=1}^K \mu_k -1 \right) = 0 \tag{2.31}
  $$
  パラメータ$\mu_k$について偏微分すると
  $$
  \frac{m_k}{\mu_k} + \lambda = 0
  $$
  $$
  \mu_k = - \frac{m_k}{\lambda} \tag{a}
  $$
  $\lambda$について偏微分すると
  　 $$
  \sum_{k=1}^K \mu_k -1 = 0
  $$
  $\mu_k$を代入して
  $$
  \sum_{k=1}^K ( - \frac{m_k}{\lambda} ) = 1
  $$
  $\lambda$は$k$によらないので$\sum$の外側に出して
$$
- \frac{1}{\lambda} \sum_{k=1}^ Km_k = 1
  $$
  $$
- \sum_{k=1}^ Km_k = \lambda
  $$
  $\sum_{k=1}^ Km_k$はサイコロの目が$k$となった回数$m_k$をすべての目に対して足し合わせている
  ことを意味する。つまり、1の目の回数+…+6の目の回数＝サイコロを振った回数であり、
  $\sum_{k=1}^ Km_k = N$である。よって、
  $$
  \lambda = - \frac{1}{N}
  $$
  これを式(a)に代入することで
  $$
  \mu_k^{\text{ML}} = \frac{m_k}{N} \tag{2.33}
  $$
  が得られる。
  結局、小難しい式を解いてきたが
  目が$k$となる確率の推定値はサイコロを振った回数$N$ 分の $k$の目が出た回数$m_k$という
  直観と一致した結果が導出できた、ということになる。

多項分布 (multinomial distribution)

それぞれの目が出る確率$\mathbf{\mu}$と観測値の総数$N$が与えられている条件で
目が1の回数が$m_1$、目が2の回数が$m_2$、…となる同時確率は
多項分布で表される。多項分布は(2.29)より$\prod_{k=1}^K \mu_k^{m_k}$に比例した形となる。
足し合わせて1になるように、係数は$N$個の対象を$m_1, \cdots, m_k$の$K$個のグループに分ける
数を使い、
$$
\text{Mult}(m_1, \cdots, m_k | \mathbf{\mu}, N) =
\frac{N!}{m_1! m_2! \cdots m_k!} \prod_{k=1}^K \mu_k^{m_k} \tag{2.34}
$$
と書ける。

2.2.1 ディリクレ分布

二項分布に対するベータ分布と同様、多項分布に対する便利な事前分布を考えよう。
多項分布と共役性を持つには、$\prod_{k=1}^K \mu_k$のべき乗で書ける関数であればよい。
そのように書ける関数、かつ足し合わせて1になるように正規化された確率分布が
ディリクレ分布である。
$$
\text{Dir}( \mathbf{\mu} | \mathbf{\alpha} ) =
\frac{ \Gamma(\alpha_0) }{ \Gamma(\alpha_1) \cdots \Gamma(\alpha_k) }
\prod_{k=1}^K \mu_k^{ \alpha_k - 1}
$$

$\alpha$は図2.5のようにディリクレ分布の形をいろいろと変えるパラメータだと
理解しておけばよい。

尤度と事前分布を書けた事後分布は
$$
\prod_{k=1}^K \mu_k^{ \alpha_k + m_k -1 }
$$
の形となり、ベータ分布の議論の際と同様、$k$の目を観測した回数$m_k$に
仮想の観測値$\alpha_k$を足して補正をかけている。

尤度関数はパラメータ$\mu$を固定値と思ったときに観測値が得られる確率。1個目の観測値が得られる確率×2個目の観測値が得られる確率×…で求まる。 ↩
表が$x=1$としていることに注意 ↩