【これなら分かる!】ベータ分布とは

zuka

こんにちは。
zuka(@beginaid)です。

本記事は「これなら分かる!はじめての数理統計学」シリーズに含まれます。

もし不適切な内容があれば,記事下のコメント欄又はお問い合わせフォームよりご連絡下さい。

ベータ分布

ベルヌーイ分布の共役事前分布として導入される分布をベータ分布と呼び,$\Be (a, b)$と表します。ベータ分布はガンマ分布に従う確率変数の変数変換を利用しても導くことができます。ベータ分布には再生性はありません。ロードマップ中では,ベータ分布は二項分布のベイズ共役に相当しています。ベータ分布にはモーメント母関数が存在することはしますが,あまり利用されていないためここでは割愛することにします。ベータ分布は$0 < X < 1$で定義され,$a$と$b$は正定数です。

\begin{align}
f_{X}(x) &= \frac{1}{B(a, b)}x^{a-1}(1-x)^{b-1} \\[0.7em]
E[X] &= \frac{a}{a+b} \\[0.7em]
V[X] &= \frac{ab}{(a+b)^2(a+b+1)}
\end{align}

ただし,$B(\cdot)$はベータ関数を表す。

\begin{align}
B(a, b) &= \int_0^1 x^{a-1}(1-x)^{b-1}dx
\end{align}

確率密度関数

ベータ分布は二項分布の共役事前分布として導入する方法と, ガンマ分布に従う確率変数の変数変換を利用して導入する方法があります。

二項分布の共役事前分布

本シリーズではここまでベイズを扱ってこなかっため,困惑している方が多いと思います。簡単におさらいしていきましょう。

ベイズを説明する王道の例はコイントスです。今から行いたいのは,コイントスを行った結果からコインの形状を推定することです。数学的に言えば,$n$回の試行のうち$k$回表が出たという結果から,コインの形状$\mu$を決めるという話です。ここで,コインの形状$\mu$というのはコイントス1回の試行を決める値,すなわちベルヌーイ分布のパラメータです。例えば,$\mu=0.7$であれば,表の出る確率が$0.7$であるような形状のコインを表しています。このパラメータ$\mu$は,複数回の試行が独立であると仮定した場合には,二項分布のパラメータと等価になります。

パラメータ推定の最も基本的な手法は最尤推定です。これは「$n$回の試行のうち$k$回表が出た」という結果のみに依存してパラメータ$\mu$を1つの値に決定する手法です。具体的には,コイントスを複数回行った結果を二項分布で表し,二項分布を定めるパラメータを変数とみて偏微分を用いて推定します。蛇足ですが,最尤推定では偏微分に加えてラグランジュの未定乗数法もよく利用されます。さて,実際に偏微分を用いてコインの形状$\mu$を推定していきましょう。まず,$n$回の試行のうち$k$回表が出るようなコインの形状を$\mu$とおけば,二項分布は以下のように表されます。

\begin{align}
{}_n C_{k}\;\mu^{k} (1-\mu)^{(n-k)}
\end{align}

こいつを尤度なんて言ったりします。この二項分布の尤度が最大となるようにパラメータ$\mu$を求めるのが最尤推定です。二項分布を$\mu$を変数とする関数$f(\mu)$とし,$f'(\mu)=0$を満たすような$\mu$が$f(\mu)$を最大にする$\mu$の候補ですから,そいつを求めに行きます。

\begin{align}
f'(\mu) &= \frac{d}{d\mu}\left(\frac{n!}{k!(n-k)!}\mu^{k} (1-\mu)^{(n-k)} \right)\\[0.7em]
&= \frac{n!}{k!(n-k)!}\left(k\mu^{(k-1)}(1-\mu)^{(n-k)} – \mu^k (n-k)(1-\mu)^{(n-k-1)} \right) \\[0.7em]
&= \frac{n!\mu^{(k-1)}(1-\mu)^{(n-k-1)}}{k!(n-k)!}\left(k(1-\mu) – \mu (n-k)\right) \\[0.7em]
&= 0 \\[0.7em]
\therefore \mu &= \frac{k}{n}
\end{align}

この結果は何を意味しているのでしょうか。試行回数における成功回数の割合がコインの形状を定めているということですね。直感にも合います。例えば4回の試行のうち3回表が出れば,0.75だけ表が出やすいコインを用いているということです。なお,通常最尤推定では計算簡略化のため対数尤度を最大化することが多いです。対数をとっても(対数関数は単調増加であることから)最大となるパラメータの値は変わらないからです。今回は対数を取らなくてもそこまで複雑にならないため生の尤度関数を用いて求めてしまいました。

さて,ここからいよいよベイズの話に突入していきます。最尤推定はコインの形状を決め打ちしてしまいます。これは点推定と呼ばれていて,求めたいパラメータを1つの値に定めることを意味しています。しかし,先ほどの例でいえば,4回の試行のうち4回表が出れば,必ず表が出るコインを用いているということになります。

…本当でしょうか。この世にはたくさんの形状のコインが存在しますが,100%表が出るようなコインなど存在するのでしょうか。というより,もはやそれってコインとは言えないですよね。そこで,得られた結果からコインの形状の「それっぽさ」を求めようとするのがベイズ推定です。最尤推定では「表の出る確率が100%のコインを使っている!」と断言していたのに対し,ベイズでは「たしかに表の出る確率が100%っぽいけど,それは60%くらいの確信度で,表の出る確率はもっと低いかもしれない。さすがに4回の試行だけでは100%とは言い切れないよな。」というような推定を行います。

それっぽさというのは,確率分布を用いて表現することができます。つまり,ベイズ推定では,ある確率分布の形状を定めるパラメータを求めることが目標になります。これを数学的に表現すると,得られた結果(尤度関数)に対して事前知識(事前分布)を導入することで,曖昧性を加味した確率分布(事後分布)のパラメータを求めるということです。この操作は,\hyperlink{定理:ベイズの定理}{ベイズの定理}によって数学的に正当性が裏付けられています。分かりやすい形で再掲しておきます。

\begin{align}
(\text{事後分布}) \propto (\text{尤度関数}) \times (\text{事前分布})
\end{align}

本題に戻りましょう。ベータ分布は二項分布の共役事前分布として導入されます。共役事前分布というのは,事後分布と形が同じになるような事前分布のことを指します。この場合,尤度関数である二項分布にかけ合わせることで同じ形が出現するような分布であることを意味しています。なぜ共役事前分布だと嬉しいのかというと,計算が楽になるからです。ベイズ推定では,事前知識を導入して事後分布を計算し,その事後分布を事前知識として改めて導入して事後分布を導入し…という更新作業を繰り返し行うことになります。そのため,事前分布と事後分布の形が同じであれば非常に都合が良いのです。さて,二項分布の共役事前分布としてのベータ分布を数式で表すと以下のようになります。

\begin{align}
(\text{ベータ分布}) \propto (\text{二項分布}) \times (\text{ベータ分布})
\end{align}

ここで,二項分布の形は分かっていますので,求めたいパラメータ$\mu$に関する部分のみ抽出して代入してみます。

\begin{align}
(\text{ベータ分布}) \propto \mu^k (1-\mu)^{(n-k)} \times (\text{ベータ分布})
\end{align}

したがって,ベータ分布は以下のような形になっていれば都合が良いことが推測されます。

\begin{align}
(\text{ベータ分布}) &\propto \mu^a (1-\mu)^b\\[0.7em]
(\text{ベータ分布}) &= C \mu^a (1-\mu)^b\quad (\because C\text{は正規化項})
\end{align}

あとは,ベータ分布が確率密度関数の定義を満たすように,積分して1になるような正規化項$C$を求めてあげるだけです。さて,まずはベータ分布の形を決める部分を$\mu$に関して積分したものを$B(a+1, b+1)$とおきます。なぜ$a+1$,$b+1$となっているかというと,その方が後々きれいな形が出現するからです。ここは定義として受け入れてください。$B(a,b)$はベータ関数と呼ばれています。$\mu$の定義域が$[0,1]$なので積分範囲も$[0,1]$になります。

\begin{align}
B(a+1, b+1) &= \int_{0}^{1}\mu^a (1-\mu)^b d\mu
\end{align}

改めてベータ分布の積分を書き直すと以下のようになり,正規化項$C$は$1/B(a+1, b+1)$に設定してあげればよいことになります。

\begin{align}
\int_{0}^{1}(ベータ分布)d\mu &= C\int_{0}^{1}\mu^a (1-\mu)^b d\mu \\[0.7em]
&= C B(a+1,b+1) \\[0.7em]
\therefore C &= \frac{1}{B(a+1, b+1)}
\end{align}

そこで,以下では$B(a+1, b+1)$を計算していきます。結論から言うと,ベータ関数はガンマ関数の積を用いて表すことができます。ガンマ関数は積分でしたね。積分の積に相性が良いのは,極座標変換です。そこで,まず,ベータ分布を極座標で表そうというモチベーションから,$\mu=\sin^2\theta$という変数変換を行います。なぜ2乗で変数変換するのかというと,ベータ関数に$1-\mu$という形が表れており,二乗で定義した方が$\sin$と$\cos$の相互変換が上手くいくからです。$\mu$の定義域が$[0,1]$なので,$\theta$の定義域は$[0, \pi/2]$になります。

\begin{align}
B(a,b) &= \int_{0}^{\frac{\pi}{2}} \sin^{2a}\theta\cos^{2b}\theta\cdot 2\sin\theta \cos\theta d\theta\\[0.7em]
&= 2\int_{0}^{\frac{\pi}{2}} \sin^{2a+1}\theta\cos^{2b+1}\theta d\theta \label{正弦波導入後のベータ関数}
\end{align}

ここで天下り的になってしまうのですが,ガンマ関数の積を考えることでベータ関数の形を出現させることができます。ガンマ関数は以下の形をしていました。

\begin{align}
\Gamma(a) &= \int_{0}^{\infty} t^{a-1}e^{-t}dt
\end{align}

ここで,極座標変換を用いて式(\ref{正弦波導入後のベータ関数})に近づけるために,$t=x^2$の変数変換を行います。

\begin{align}
\Gamma(a) &= \int_{0}^{\infty} x^{2a-2}e^{-x^2}2xdx\\[0.7em]
&= 2\int_{0}^{\infty} x^{2a-1}e^{-x^2}dx
\end{align}

同様に,$\Gamma(b)$に対して$t=y^2$の変数変換を行います。

\begin{align}
\Gamma(b) &= \int_{0}^{\infty} y^{2a-2}e^{-y^2}2ydy\\[0.7em]
&= 2\int_{0}^{\infty} y^{2a-1}e^{-y^2}dy
\end{align}

ここで,2つのガンマ分布の積を計算してみます。

\begin{align}
\Gamma(a)\Gamma(b) &=
4\int_{0}^{\infty} x^{2a-1}e^{-x^2}dx
\int_{0}^{\infty} y^{2a-1}e^{-y^2}dy\\[0.7em]
&= 4\int_{0}^{\infty} \int_{0}^{\infty}
x^{2a-1}y^{2b-1}e^{-(x^2+y^2)} \label{極座標変換前のガンマ分布の積}
dxdy
\end{align}

$x^2+y^2$という極座標変換において出てきて嬉しい項が出現しました。さて,満を持して極座標変換を行います。$x=r\cos\theta, y=r\sin\theta$とおきます。このとき,以下のようにヤコビアン$J$は$r$となります。

\begin{align}
|J| &=
\begin{vmatrix}
\partial x / \partial r & \partial x / \partial \theta \\
\partial y / \partial r & \partial y / \partial \theta \\
\end{vmatrix}\\[0.7em]
&=
\begin{vmatrix}
\cos \theta & -r\sin\theta \\
\sin \theta & r\cos \theta \\
\end{vmatrix}\\[0.7em]
&= r \cos^2\theta + r \sin^2\theta
= r
\end{align}

したがって,式(\ref{極座標変換前のガンマ分布の積})は以下のように変換されます。

\begin{align}
\Gamma(a)\Gamma(b)
&= 4\int_{0}^{\infty} \int_{0}^{\frac{\pi}{2}}
r^{2a-1}\cos^{2a-1}\theta \; r^{2b-1}\sin^{2b-1}\theta
e^{-r^2}rdrd\theta \\[0.7em]
&= 2\int_{0}^{\infty} r^{2(a+b)-1} e^{-r^2}dr \cdot
2\int_{0}^{\frac{\pi}{2}}\cos^{2b-1}\theta \sin^{2a-1}\theta d\theta\\[0.7em]
&= \Gamma(a+b)B(a, b)
\end{align}

はい。ここできれいな形を得るために,冒頭で無理矢理$B(a+1, b+1)$と定義したのです。

\begin{align}
B(a-1,b-1) &= \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
\end{align}

以上をまとめると,ベータ分布の確率密度関数$f(\mu)$は以下のようになります。

\begin{align}
f(\mu) &= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}
\mu^{a-1} (1-\mu)^{b-1}
\end{align}

ガンマ分布に従う確率変数の変数変換

ベータ分布はガンマ分布に従う確率変数の変数変換からも求めることができます。この手法はディリクレ分布の導出にも利用されるために,拡張性が高いです。

天下り的ですが,$U$と$V$が独立かつ

\begin{align}
U &\sim \mathrm{Ga}(a, c) \\[0.7em]
V &\sim \mathrm{Ga}(b, c)
\end{align}

とします。ここで,

\begin{align}
X &= \frac{U}{U+V} \\[0.7em]
Y &= U+V
\end{align}

を考えると,実は

\begin{align}
X &\sim B(a,b) \\[0.7em]
Y &\sim \mathrm{Ga}(a+b, c)
\end{align}

となります。実際に,同時確率密度関数$f(x, y)$を求めることで示していきます。今回は変数変換の逆関数である

\begin{align}
u &= xy \\[0.7em]
v &= (1-x)y
\end{align}

が存在しますので,変数変換を使います。$f_{U}(u)$と$f_{V}(v)$が独立であることに注意しましょう。

\begin{align}
f_{X, Y}(x, y) &= f_{U, V}(u, v) |J(u,v)|\\[0.7em]
&= f_{U}(u)f_{V}(v) |J(u,v)|
\end{align}

ヤコビアンを求めましょう。

\begin{align}
|J(u, v)| &=
\begin{vmatrix}
\partial u / \partial x & \partial u / \partial y \\
\partial v / \partial x & \partial v / \partial y \\
\end{vmatrix}\\[0.7em]
&=
\begin{vmatrix}
y & x \\
-y & 1-x \\
\end{vmatrix}\\[0.7em]
&= y
\end{align}

したがって,同時確率分布は以下のように計算できます。

\begin{align}
f_{X, Y}(x, y)
&= \frac{c^a}{\Gamma(a)} (xy)^{a-1}e^{-cxy}\cdot
\frac{c^a}{\Gamma(b)} ((1-x)y)^{b-1}e^{-c(1-x)y}\cdot
y\\[0.7em]
&= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1}(1-x)^{b-1}\cdot
\frac{c^{a+b}}{\Gamma(a+b)}y^{a+b-1}e^{-cy}
\end{align}

これを噛み砕くと,以下のような形になっています。

\begin{align}
f_{X, Y}(x, y) &= (\text{Xを変数とするベータ分布}) \times (\text{Yを変数とするガンマ分布})
\end{align}

これは,$X, Y$が独立で$X\sim B(a, b)$,$Y\sim \mathrm{Ga}(a+b, c)$を表しています。

ベータ分布の上側確率

さらに,ベータ分布の上側確率は二項分布の確率関数の部分和と等しいという性質があります。これは,ガンマ分布のときと同様に部分積分を用いて示すことができます。

\begin{align}
&\int_{p}^{1}\frac{z^{k-1}(1-z)^{n-k}}{B(k,n-k+1)}dz \notag\\[0.7em]
&= \left[\frac{-1}{n-k+1}\frac{z^{k-1}(1-z)^{n-k+1}}{B(k, n-k+1)} \right]_{p}^{1} + \frac{k-1}{n-k+1} \int_{p}^{1}\frac{z^{k-2}(1-z)^{n-k+1}}{B(k, n-k+1)} dz \\[0.7em]
&= \frac{n!}{(k-1)!(n-k)!} p^{k-1}(1-p)^{n-k}\frac{(k-1)}{(n-k+1)}\cdot \frac{\Gamma(n+1)}{\Gamma(k)\Gamma(n-k+1)} \int_{p}^{1} z^{k-2}(1-z)^{n-k+1}dz \\[0.7em]
&= {}_n C_{n-k}\; p^{k-1}(1-p)^{n-k}
+ \frac{(k-1)}{(n-k+1)}\cdot \frac{\Gamma(n+1)(n-k+1)}{(k-1)\Gamma(k-1)\Gamma(n-k+2)}
\int_{p}^{1} z^{k-2}(1-z)^{n-k+1}dz \\[0.7em]
&= {}_n C_{n-k}\; p^{k-1}(1-p)^{n-k}
+ \frac{\Gamma(n+1)}{\Gamma(k-1)\Gamma(n-k+2)}
\int_{p}^{1} z^{k-2}(1-z)^{n-k+1}dz \\[0.7em]
&= {}_n C_{n-k}\;p^{k-1}(1-p)^{n-k+1} +
\int_{p}^{1}\frac{z^{k-2}(1-z)^{n-k+1}}{B(k-1, n-k+2)}dz
\end{align}

この操作を次々適用していくことで,最終的に以下を得ます。

\begin{align}
\int_{p}^{1}\frac{z^{k-1}(1-z)^{n-k}}{B(k,n-k+1)}dz
&= \sum_{x=0}^{k-1} {}_n C_{x}p^x (1-p)^{n-x}
\end{align}

モーメント母関数

通常,モーメント母関数はモーメント母関数の定義に従って計算していきます。しかし,ベータ分布のモーメント母関数は複雑な計算を要するため,あまり利用されていません。本シリーズでも割愛することにします。

平均・分散

連続分布の平均と分散を求める際には「モーメント母関数の性質」(離散分布の平均と分散を求める際には「確率母関数の性質」)を利用します。

しかし,ベータ分布のモーメント母関数は複雑な形をしているため,ここでは平均と分散の定義に従って愚直に計算していきます。これはパズル問題的な感覚でスルスルと求めていけますので心配はいりません。利用するのは以下の関係式のみです。

\begin{align}
B(a,b) &= \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} \\[0.7em]
\Gamma(a+1) &= a\Gamma(a)
\end{align}

実際に計算していきます。

\begin{align}
E[X] &= \frac{1}{B(a, b)}\int_{0}^{1}x\cdot x^{a-1}(1-x)^{b-1} dx \\[0.7em]
&= \frac{1}{B(a, b)}\cdot B(a+1, b) \\[0.7em]
&= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \cdot
\frac{\Gamma(a+1)\Gamma(b)}{\Gamma(a+b+1)} \\[0.7em]
&= \frac{\Gamma(a+b)}{\Gamma(a)}\cdot
\frac{a\Gamma(a)}{(a+b)\Gamma(a+b)} \\[0.7em]
&= \frac{a}{a+b}
\end{align}

まさに,ガンマ関数を階乗的な感覚で次数を落としていき,約分していくパズルです。同様に$E[X^2]$も求めていきます。

\begin{align}
E[X^2] &= \frac{1}{B(a, b)}\int_{0}^{1}x^2\cdot x^{a-1}(1-x)^{b-1} dx \\[0.7em]
&= \frac{1}{B(a, b)}\cdot B(a+2, b) \\[0.7em]
&= \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \cdot
\frac{\Gamma(a+2)\Gamma(b)}{\Gamma(a+b+2)} \\[0.7em]
&= \frac{\Gamma(a+b)}{\Gamma(a)}\cdot
\frac{(a+1)a\Gamma(a)}{(a+b+1)(a+b)\Gamma(a+b)} \\[0.7em]
&= \frac{(a+1)a}{(a+b+1)(a+b)}
\end{align}

したがって,以下のように分散を計算できます。

\begin{align}
V[X] &= E[X^2] – E[X]^2 \\[0.7em]
&= \frac{(a+1)a}{(a+b+1)(a+b)} – \left(\frac{a}{a+b}\right)^2 \\[0.7em]
&= \frac{ab}{(a+b)^2(a+b+1)}
\end{align}

再生性

再生性を示すためには,再生性を示したい分布に従う独立な二つの確率変数を考え,その和のモーメント母関数(離散分布の場合はモーメント母関数)を計算したときに,パラメータが和の形になっていることを示します。

ベータ分布のモーメント母関数は複雑な形をしており,積をとっても同じ形が出現したいため,再生性はありません。

ロードマップ

確率分布のロードマップ

さて,ロードマップに戻りましょう。 ベータ分布は二項分布のベイズ共役(共役事前分布)として導出できましたね。また,ガンマ分布に従う確率変数を利用した変数変換を通じても導出することができました。また,ベータ分布の上側確率は二項分布の確率関数の部分和とも関係がありましたね。

シェアはこちらからお願いします!
URLをコピーする
URLをコピーしました!

コメント

コメントする