演習a 第１回発表

Python によるベイズ統計学入門2.1~3.1

2023-04-14

Kaoru Babasaki* and Keigo Takahashi*

*Keio University, Graduate school of Economics

ベイズ統計学による推定の概要

真の値が未知であるパラメータ $q$ を確率変数であると解釈し,その確率分布を使ってパラメータの値に関する統計的推論を行う.
具体的には,最初に推論の出発点となるパラメータの確率分布として事前分布を想定し,これをデータが持つパラメータに関する情報で更新することで事後分布を求めてパラメータの値に関する推論を行うという手順を踏む.

ベイズの定理・用語

データ $D=\{x_1,\cdots,x_n\}$ 下での パラメータ $q$ の確率分布 （事後分布） は

条件付き確率の定義から

$\begin{align*} p(q|D) &= \frac{p(q,D)}{p(D)} \\ &= \frac{p(D|q)p(q)}{p(D)} \\ &= \frac{p(D|q)P(q)}{\int p(D|q)p(q) dq} \end{align*}$

$p(q)$ を事前分布・ $p(D|q)$ を尤度・ $p(D)=\int p(D|q)p(q) dq$ を周辺尤度と呼ぶ.

※ パラメータ とは,分布の特性・形状を決定する変数のことである.

※ パラメータを確率変数をみなし,分布として推定するのがベイズ統計学の特徴の一つ.

尤度 $p(D|q)$ と周辺尤度 $p(D)$

パラメータ $q$ の具体的な値が与えられた下での,データ群 $D$ が観測される確率

要するにこれは独立を仮定して,各事象が観測される確率の積をとったもの.

今回はベルヌーイ分布を想定しているので以下の様な計算になる.

$\begin{align*} P(D|q) &= Pr\{X_1 = x_1 , \dots , X_n = x_n|q\} &= \Pi_{i-1}^{n} Pr\{X_i = x_i|q\} \\ &= \prod_{i=1}^n q^{x_i}(1-q)^{1-x_i} \\ &= q^y(1-q)^{n-y} , y = \Sigma_{i=1}^{n} x_i \end{align*}$

※ $p(D) = \int_{q=0}^{1} p(D,q) dq = \int_{q=0}^{1} p(D|q)p(q) dq = (\Sigma_{q=0}^{1} p(D|q)p(q) )$

$p(D)$ はパラメータ $q$ の値域で積分・足し合わせたものである.

これはデータ $D$ が観測される平均的な可能性を指しており $q$ の値に依らない

内閣支持率 $q$ の事後分布の導出

事前分布はベータ分布,パラメータは $\alpha_0,\beta_0$ とする.

$q$ の事後分布はパラメータ $q$ にまつわる部分(比例記号)を抜き出すと

$\begin{align*} p(q|D) &\propto p(D|q)p(q) \\ &\propto q^{y}(1-q)^{n-y} \times \frac{q^{\alpha_0-1}(1-q)^{\beta_0-1}}{B(\alpha_0,\beta_0)} \\ &\propto q^{y+\alpha_0-1}(1-q)^{n-y+\beta_0-1} , y = \Sigma x_i \end{align*}$

基準化定数 $K$ を, $\int_{0}^{1} p(q|D) dq = 1$ の条件から求めると $K = B(y+\alpha_0,n-y+\beta_0)$

よって具体的な事後分布 $p(q|D)$ の形状は ベータ分布 となる.

$p(q|D) = \frac{q^{y+\alpha_0-1}(1-q)^{n-y+\beta_0-1}}{B(y+\alpha_0,n-y+\beta_0)}$

ベータ分布に従うという意味で $q|D \sim Be(y+\alpha_0,n-y+\beta_0) , y = \Sigma x_i$

概要

ベイズ統計学は,パラメータ $(=q)$ を確率変数とみなし,分布を出す. そして分布の形状からあたりをつける.

「あたり」を定量的に推定・議論をしたい. またパラメータに対する仮説を検証したい.

点推定
区間推定
仮説検定

ベイズ統計学における以上の概念の分析をこの節で紹介している.

代表的な損失関数 $L(q,\delta)$ と点推定の値 $\delta$

ベイズ統計学で使われる代表的な損失関数と対応する点推定の値は以下の表の通り

損失関数 $L(q,\delta)$	点推定 $\delta$
$(q-\delta)^2$	事後分布の平均
$q-\delta$	事後分布の中央値
$1-\unicode{x1D7D9}_{q}(\delta)$	事後分布の最頻値

$\begin{align*} \unicode{x1D7D9}_{q}(\delta) = \begin{cases} 1 & (\delta = q)\\ 0 & (\delta \not= q) \end{cases} \end{align*}$

補足事後分布の最頻値 $Mode_q$

$Mode_q = argmax_{0 \leq q \leq 1} p(q|D)$

最頻値による点推定を MAP 推定 と呼ぶ.（機械学習の文脈で何回かみたことがある)

メリット

周辺尤度 $P(D)$ の計算が必要がなく,計算が楽

パラメータ $q$ が表れるのは分子部分 $p(D|q)p(q)$ だけなので, $argmax_{0 \leq q \leq 1} p(D|q)p(q)$ を考えればよい.

$argmax_{0 \leq q \leq 1} q^{y+\alpha_0-1}(1-q)^{n-y+\beta_0-1}$
$\frac{\partial q^{y+\alpha_0-1}(1-q)^{n-y+\beta_0-1}}{\partial q} = 0$
積の微分の公式と合成関数の微分の公式の 2 つを使って計算→次スライド.

デメリット

事後分布の最頻値が一つであるとは限らない.峰が複数ある分布の可能性もあるので,簡易的に採用するのは NG.

補足：微分計算

積の微分と合成関数の微分を使う.

$\begin{align*} \frac{\partial q^{\alpha-1}(1-q)^{\beta-1}}{\partial q} &= (q^{\alpha-1})^{\prime}(1-q)^{\beta-1} + q^{\alpha-1}((1-q)^{\beta-1})^{\prime} \\ &= (\alpha-1)q^{\alpha-2}(1-q)^{\beta-1} - q^{\alpha-1}(\beta-1)(1-q)^{\beta-2} \\ &= q^{\alpha-2}(1-q)^{\beta-2}((\alpha-1)(1-q) - (\beta-1)q ) \end{align*}$

$\frac{\partial q^{y+\alpha_0-1}(1-q)^{n-y+\beta_0-1}}{\partial q} = 0$ より, $(\alpha-1)(1-q) - (\beta-1)q = 0$ であれば良い.

$\begin{align*} (\alpha-1)(1-q) - (\beta-1)q &= -(\alpha-1)+(-\alpha-\beta+2)q \\ (\alpha+\beta-2)q &= (\alpha-1) \\ q &= \frac{\alpha-1}{\alpha + \beta+2} \end{align*}$

各損失関数における点推定 $\delta$ の導出

そもそもパラメータが推定対象であり真の値は知り得ない.従って $L(q,\delta)$ は計算できない.

なので, $L(q,\delta)$ の期待値を事後分布で評価したもの $E_q[L(q,\delta)|D] = \int L(q,\delta)p(q|D) dq$ を考え,それを最小化する $\delta$ を考える.

※ $E_q[L(q,\delta)|D]$ は,データ $D$ が与えられた下で $q$ について評価したもので条件付き期待値を指す.

2 乗誤差損失→事後分布の平均値

$L(q,\delta) = (q-\delta)^2$ より,

$\begin{align*} E_q[(q-\delta)^2|D] &= E_q[(q-E_q[q|D]+E_q[q|D]-\delta)^2|D] \\ &= E_q[(q-E_q[q|D])^2-2(q-E_q[q|D])(E_q[q|D]-\delta)+(E_q[q|D]-\delta)^2|D] \\ &= E_q[(q-E_q[q|D])^2|D] + (E_q[q|D]-\delta)^2 \end{align*}$

※ $E_q[-2(q-E_q[q|D])(E_q[q|D]-\delta)|D]$ は実は 0 になる. 計算してみる

$\delta$ の値が含まれる(今は期待損失を最小にするような $\delta$ を探している.)のは, $(E_q[q|D]-\delta)^2$ .

$(E_q[q|D]-\delta)^2$ を最小化すればよく,2 次関数であるからその最小点は

$\delta^* = E_q[q|D]$

※ $(q-\delta)^2 = (q-E_q[q|D]+E_q[q|D]-\delta)^2$ が肝

絶対誤差損失→事後分布の中央値

$L(q,\delta) = |q-\delta|$ , $p(q|D)$ の累積分布関数を $P(q|D)$ とする.

※ 累積分布関数 $P(\delta) = \int^{\delta}_{-\infty}p(x)dx$

$\begin{align*} E_q[|q-\delta||D] &= \int_0^1 |q-\delta|p(q|D)dq \\ &= \int_{0}^{\delta} (\delta-q)p(q|D)dq + \int_{\delta}^{1} (q-\delta)p(q|D)dq \\ &= \delta*\int_{0}^{\delta}p(q|D)dq - \int_{0}^{\delta}qp(q|D)dq + \int^1_{\delta} qp(q|D)dq + -\delta*\int^1_{\delta} p(q|D)dq \\ &= \delta P(\delta|D) - \int_{0}^{\delta}qp(q|D)dq + \int^1_{\delta} qp(q|D)dq -\delta(1-P(\delta|D)) \\ &= 2 \delta P(\delta|D) -\delta - \int_{0}^{\delta}qp(q|D)dq +(\int_{0}^{1}qp(q|D)dq - \int_{0}^{\delta}qp(q|D)dq) \\ \end{align*}$

絶対誤差損失→事後分布の中央値続き

式展開の続き

$\begin{align*} E_q[|q-\delta||D]&= 2 \delta P(\delta|D) -\delta - \int_{0}^{\delta}qp(q|D)dq +(\int_{0}^{1}qp(q|D)dq - \int_{0}^{\delta}qp(q|D)dq) \\ &= 2 \delta P(\delta|D) -\delta -2 \int_{0}^{\delta}qp(q|D)dq + \int_{0}^{1}qp(q|D)dq \\ &= 2 \delta P(\delta|D) -\delta -2\{\mid qP(q|D) \mid^{\delta}_{0} - \int_{0}^{\delta} P(q|D)dq\} + E_q[q|D] \\ &= 2 \int_{0}^{\delta} P(q|D)dq - \delta + E_q[q|D] \end{align*}$

最後の式の最小化を考えると,1 階の条件は $\nabla_{\delta}R(\delta^{*}|D) = 2P(\delta^{*}|D) -1 = 0$

$P(\delta^{*}|D) = \frac{1}{2}$ となり $\delta^{*}$ は事後分布の中央値となる.

※ 最適解なのかどうかをちゃんと議論する際,2 階の条件まで考える必要がある.

0-1 損失→事後分布の最頻値

$L(q,\delta) = 1 - \unicode{x1D7D9}_{q}(\delta)$

ただし $\delta$ の定義域で少し異なる損失関数を定義する.

$\begin{align*} L_{\epsilon}(q,\delta) &= 1 - \unicode{x1D7D9}_{[q-\epsilon,q+\epsilon]}(\delta) \\ &= \begin{cases} 0 & (q-\epsilon \leq \delta \leq q+\epsilon)\\ 1 & (\delta < q-\epsilon, q+\epsilon < \delta) \end{cases} \end{align*}$

0-1 損失→事後分布の最頻値続き

$\begin{align*} E_q[L_{\epsilon}(q,\delta)|D] &= \int_{0}^{1} L_{\epsilon}(q,\delta)p(q|D)dq \\ &= \int_{0}^{1} p(q|D)dq - \int_{0}^{1} \unicode{x1D7D9}_{[q-\epsilon,q+\epsilon]}(\delta) p(q|D)dq \\ &= 1 - \int_{\delta-\epsilon}^{\delta+\epsilon} p(q|D)dq \end{align*}$

絶対誤差を最小化するためには, $\int_{\delta-\epsilon}^{\delta+\epsilon} p(q|D)dq$ を最大にする $\delta$ を求める.

$Pr\{\delta-\epsilon \leq \delta+\epsilon|D\}$ という狭い区間の積分値を大きくすることになるため,最頻値になる.

区間推定

点推定は,パラメータの真の値が分からず,期待損失を最小化している以上不確実性が残る.

区間推定は,真の値が含まれる可能性の高い区間を提示する.

$Pr\{a \leq q \leq b|D\} = \int_{a}^{b}p(q|D)dq$

$Pr\{a \leq q \leq b|D\} = 0.95,0.90 \cdots$ となるような $(a,b)$ を求めれば良いが,この条件だけでは解の組が無数に存在するため,他の制約条件を考える必要がある.

一般に以下 2 つの区間が使用されることが多い.

信用区間
HPD(highest posterior density)区間

HPD(highest posterior density)区間

$c$ を区間の外側の確率とする.区間に真の値がある確率は $100(1-c)$

このとき, $100(1-c)$ HPD 区間 $[a_{hpd},b_{hpd}]$ は

$\begin{align*} Pr\{a_{hpd} \leq q \leq b_{hpd}|D\} & = 1-c \\ p(a_{hpd}|D) &= p(b_{hpd}|D) \end{align*}$

を満たす.

信用区間との違い

区間内の確率密度が区間外の確率密度より必ず高くなる
区間外の 2 つの面積が等しいとは限らない

仮説検定

パラメータの推定に加えて,パラメータに関する仮説を検討したいときもある.(例:内閣支持率 $q$ 3 割超えているか？)

仮説 $H_i$ はパラメータのとりうる範囲・特定の値か否かで定義される.

ある範囲内の値をとる (例 $\{q : 0.5 \leq q \leq 1\}$ )
特定の値に等しい (例 $q = 0.5$ )
特定の値に等しくない (例 $q \neq 0.5$ )

「ある範囲内の値をとる」という仮説については,範囲を $S_i$ とすると,仮説 $H_i : q \in S_i$ の妥当性は

$Pr\{H_i|D\} = Pr\{q \in S_i|D\} = \int_{S_i}p(q|D)dq$

範囲 $S_i$ の事後分布における確率を計算するだけである.

1 に近い値なら仮説が成り立ち,逆に 0 に近い値なのであればその仮説は成り立たないと判断.

事後オッズ比・ベイズファクター $B_{01}$

ここでは 2 つの仮説 $H_0,H_1$ の比較を前提とした計算を紹介する.(2 つ以上の仮説の比較は…?)

事後オッズ比は,範囲 $S_i$ の事後分布における確率の比をとる
ベイズファクター $B_{01}$ は,事前分布による仮説支持のバイアスを取り除くために,範囲 $S_i$ の事前分布における確率を事後オッズ比に組み入れて計算.

$\begin{align*} 事後オッズ比 &= \frac{Pr\{H_0|D\}}{Pr\{H_1 | D\}} = \frac{Pr\{q \in S_0 \mid D\}}{Pr\{q \in S_1 \mid D\}} \\ B_{01} &= \frac{\frac{Pr\{H_0|D\}}{Pr\{H_0\}}}{\frac{Pr\{H_1|D\}}{Pr\{H_1\}}} \\ &= \frac{Pr\{H_0|D\}}{Pr\{H_1|D\}} \div \frac{Pr\{H_0\}}{Pr\{H_1\}} \end{align*}$

ただ,ベイズファクターは事後分布がどれだけ値に貢献しているのか判断しにくい.

ベイズファクター $B_{01}$ の評価

常用対数値 $\log_{10}B_{01}$ の値を使って,Jeffreys が,事後分布においてどちらの仮説を支持するかの目安を出している.

等級	ベイズファクターの常用対数値	$H_1$ に対する支持
0	$0 < \log_{10}B_{01}$	$H_0$ が支持される
1	$-\frac{1}{2} < \log_{10}B_{01} < 0$	それほどではない
2	$-1 < \log_{10}B_{01} < -\frac{1}{2}$	相当なものである
3	$-\frac{3}{2} < \log_{10}B_{01} < -1$	強い
4	$-2 < \log_{10}B_{01} < -\frac{3}{2}$	かなり強い
5	$\log_{10}B_{01} < -2$	決定的である

特定の値に等しい・等しくないの仮説の検証

特定の区域内に存在する・ある値以上・以下の他に,特定の値に等しい・等しくないという仮説を検証したいときもある. 例）回帰係数 $\beta$ が 0 に等しい・等しくない.

ベイズ・ファクター $B_{01}$ を使って,等しい仮説 $H_0$ ・等しくない仮説 $H_1$ を比較・検証する.

$\begin{align*} \begin{cases} H_0 : & q = q_0 \\ H_1 : & q \neq q_0 \\ \end{cases} \end{align*}$

ただし, $Pr \{ q = q_0 \mid D \} , Pr \{ q = q_0\}$ は 0 であり,このままでは実はベイズファクターが計算できない.

※ 連続型確率分布では,特定の値における確率(密度)は 0 である.

SDDR(Savage-Dickey Density Ratio)

ディラックのデルタ関数 $\delta(x)$ を組み入れた $p(q) = p_0\delta(q-q_0) + (1-p_0)f(q)$ を事前分布に想定することで回避する.

$p_0$ は仮説 $H_0$ が正しい確率を表す.

ベイズ・ファクターの計算は,以上の事前分布を想定すると

$B_{01} = \frac{f(q_{0}|D)}{f(q_{0})} , f(q|D) = \frac{p(D|q)f(q)}{\int_{0}^{1}p(D|q)f(q)dq}$

$f(x)$ は,ディラックのデルタ関数を組み入れていない元の事前分布？(ここではベータ分布)
元の事前分布とそこから計算される事後分布における $q_0$ の確率密度を比較

SDDR(Savage-Dickey Density Ratio)の導出

事前分布 $p(q) = p_{0} \delta(q-q_{0}) + (1-p_{0})f(q)$

※ $\delta(q-q_{0})$ は $q \neq q_{0}$ のときは 0, $q=q_{0}$ のときは特に値を考えず, $q=q_{0}$ という事実が必要.

事前オッズ比

$\begin{align*} \frac{Pr\{q=q_0\}}{Pr\{q \neq q_0\}} = \frac{p_0}{1-p_0} \end{align*}$

※ $p_0$ は仮説 $H_0=\{q=q_0\}$ が正しい確率を表す.

SDDR(Savage-Dickey Density Ratio)の導出続き

事後オッズ比を求める.

$q$ の事後分布は

$\begin{align*} p(q|D) &= \frac{p(D|q)p(q)}{\int_{0}^{1}p(D|q)p(q)dq} \\ &= \frac{p(D|q)\{p_0\delta(q-q_0) + (1-p_0)f(q)\}}{\int_{0}^{1}p(D|q)\{p_0\delta(q-q_0) + (1-p_0)f(q)\}dq} \\ &= \frac{p_{0}p(D|q) \delta(q-q_0) + (1-p_{0})p(D|q)f(q)}{p_{0}\int_{0}^{1}p(D|q)\delta(q-q_{0})dq + (1-p_{0}) \int_{0}^{1}p(D|q)f(q)dq} \\ &= \frac{p_{0}p(D|q) \delta(q-q_0) + (1-p_{0})p(D|q)f(q)}{p_{0}p(D|q_0) + (1-p_0) \int_{0}^{1}p(D|q)f(q)dq} \\ \end{align*}$

3 行目から 4 行目にかけては,ディラックのデルタ関数の性質

任意の連続関数 $g(x)$ に対して, $\int_{-\infty}^{\infty}g(x)\delta(x)dx=g(0)$ を使用.

SDDR(Savage-Dickey Density Ratio)の導出続き

事後オッズ比は,分母の $p_{0}p(D|q_0) + (1-p_0) \int_{0}^{1}p(D|q)f(q)dq$ はキャンセルされるため

$\begin{align*} \frac{Pr\{q=q_{0}|D\}}{Pr\{q \neq q_{0}|D\}} &= \frac{p_{0}}{1-p_{0}} \times \frac{p(D|q_0)}{\int_{0}^{1}p(D|q)f(q)dq} \end{align*}$

ベイズ・ファクター $B_{01}$ は,今までの事前オッズ比・事後オッズ比より

$\begin{align*} B_{01} &= \frac{Pr\{q=q_{0}|D\}}{Pr\{q \neq q_{0}|D\}} \div \frac{Pr\{q=q_0\}}{Pr\{q \neq q_0\}} \\ &= \frac{p(D|q_{0})}{\int_{0}^{1}p(D|q)f(q)dq} \end{align*}$

仮説 $H_0$ が正しい確率 $p_{0}$ はベイズ・ファクターを求める際には不要!

SDDR(Savage-Dickey Density Ratio)の導出終

最後に,分子と分母に $f(q_0)$ をかけると SDDR の完成.

$\begin{align*} B_{01} &= \frac{p(D|q_{0})}{\int_{0}^{1}p(D|q)f(q)dq} \times \frac{f(q_{0})}{f(q_{0})} \\ &= \frac{p(D|q_{0})f(q_{0})}{\int_{0}^{1}p(D|q)f(q)dq} \times \frac{1}{f(q_{0})} \\ &= \frac{f(q_{0} | D)}{f(q_{0})} \end{align*}$

※ $f(q_{0} | D) = \frac{p(D|q)f(q)}{\int_{0}^{1}p(D|q)f(q)dq}$

この節のポイント

データを $D=(x_1, \dots, x_n)$ , 将来観測されるであろうデータを $\tilde{x}$ としたとき

$\tilde{x} \perp D$ ならば、 $\tilde{x}$ の分布は、以下で表せる。 $p(\tilde{x}|D) = \int_{\theta \in \Theta} p(\tilde{x}|\theta)p(\theta|D) d\theta$

解釈：「 $p(\tilde{x} \mid \theta)$ を、事後分布 $p(\theta \mid D)$ で重みをつけて平均したもの」

$\theta$ : 仮定した sampling model のパラメータ
$\Theta$ : $\theta$ の取りうる範囲の集合。

この形になると何が嬉しいのか

$p(\tilde{x}|D) = \int_{\theta \in \Theta} p(\tilde{x}|\theta)p(\theta|D) d\theta$

いま手元にあるデータ $D=(x_1, \dots, x_n)$ から、未観測のデータ $\tilde{x}$ の分布を予測したい。
そこで、 $\tilde{x}$ が何らかのパラメータ $\theta_{*}$ に依存する確率分布 $p( \tilde{x} \mid \theta_*)$ に従うと仮定する。
しかし、 $\theta_*$ は未知だから $p( \tilde{x} \mid \theta_*)$ は直接求められない。
一方多くの場合で、
- 確率（密度）関数
- 事後分布
は形がわかる(モデルで仮定している)
そこで、 $D$ を生かして $\theta_*$ を推定し、(← まさに前半でやってきたこと！！)
わからないパラメータに関しては積分で消えてもらう。

設定

$\begin{align*} \text{観測されたデータ: }&D = (x_1, \dots, x_n) \\ &\qquad \text{where }x_1, \dots, x_n \overset{\text{i.i.d.} }{\sim} \text{Bernoulli}(q) \\ \\ \text{将来観測されるであろうデータ: }&\tilde{x} (\perp D) \end{align*}$

条件付確率分布の性質より $p(\tilde{x}, D) = p(\tilde{x} \mid D) p(D) \ \Rightarrow \ p(\tilde{x} \mid D) = \frac{p(\tilde{x}, D)}{p(D)}$

(補足)Hoff の出し方との違い

Hoff の本では、

$\begin{align*} p(\tilde{x}|D) & = \int_{\theta \in \Theta} p(\tilde{x}, \theta \mid D) d\theta \\ & = \int_{\theta \in \Theta} p(\tilde{x}|\theta, D)p(\theta|D) d\theta \\ & = \int_{\theta \in \Theta} p(\tilde{x}|\theta)p(\theta|D) d\theta \quad (\because \tilde{x} \text{ and }D \text{ are conditionally i.i.d.}) \end{align*}$

(どちらもやってることは同じだが、中妻先生は周辺尤度を使いたかったらしい。)

ベルヌーイ分布の例

$\begin{align*} p(\tilde{x} \mid D) &= \int_0^1 p(\tilde{x} \mid q) p(q \mid D) dq \\ &= \int_0^1 q^{\tilde{x}} (1-q)^{1-\tilde{x}} \frac{q^{\alpha_* - 1} (1-q)^{\beta_* - 1}}{B(\alpha_*, \beta_*)} dq \\ &= \frac{1}{B(\alpha_*, \beta_*)} \int_0^1 q^{\tilde{x} + \alpha_* - 1} (1-q)^{\beta_* -\tilde{x} + 1- 1} dq \\ &= \frac{B(\alpha_* + \tilde{x}, \beta_* - \tilde{x} + 1)}{B(\alpha_*, \beta_*)} \qquad \text{where } \alpha_*, \beta_* \cdots \text{ 事後分布のパラメータ} \end{align*}$

よって、

$\begin{align*} & p(\tilde{x} = 1 \mid D) = \frac{B(\alpha_* + 1, \beta_*)}{B(\alpha_*, \beta_*)} & p(\tilde{x} = 0 \mid D) = \frac{B(\alpha_*, \beta_* + 1)}{B(\alpha_*, \beta_*)} \end{align*}$

補足：ガンマ関数

$\Gamma(\alpha)$ はガンマ関数(gamma function)と呼ばれる関数で、 $\Gamma(\alpha) = \int_0^\infty y^{\alpha - 1} e^{-y} dy$ で与えられる。

ガンマ関数の重要な性質

(1)
$\Gamma(\frac{1}{2}) = \sqrt{\pi}, \quad \Gamma(1)=1$

(2)
$\Gamma(\alpha + 1) = \alpha \Gamma(\alpha)$ , 特に $n$ が自然数のときには $\Gamma(n)= (n-1)!$

証明は、久保川(2017) p.44 を参照。

ベータ関数の性質

$\begin{align*} B(\alpha + 1, \beta) &= \frac{\Gamma(\alpha + 1) \Gamma(\beta)}{\Gamma(\alpha + 1 + \beta)} \qquad(\because \text{ ベータ関数の定義 }) \\ &= \frac{\alpha \Gamma(\alpha) \Gamma(\beta)}{(\alpha + \beta)\Gamma(\alpha + \beta)} \qquad(\because \text{ ガンマ関数の性質(2) }) \\ &= \frac{\alpha}{\alpha + \beta} B(\alpha, \beta) \qquad(\because \text{ ベータ関数の定義 }) \end{align*}$

同様にして、

$\begin{align*} B(\alpha, \beta + 1) = \frac{\beta}{\alpha + \beta} B(\alpha, \beta) \end{align*}$

ベータ分布の例

よって、

$\begin{align*} & p(\tilde{x} = 1 \mid D) = \frac{B(\alpha_* + 1, \beta_*)}{B(\alpha_*, \beta_*)} = \frac{\alpha_*}{\alpha_* + \beta_*} \\ & p(\tilde{x} = 0 \mid D) = \frac{B(\alpha_*, \beta_* + 1)}{B(\alpha_*, \beta_*)} = \frac{\beta_*}{\alpha_* + \beta_*} \end{align*}$

であるから、予測分布は、

$\begin{align*} p(\tilde{x} \mid D) &= \left( \frac{\alpha_*}{\alpha_* + \beta_*} \right)^{\tilde{x}} \left( \frac{\beta_*}{\alpha_* + \beta_*} \right)^{1 - \tilde{x}} \\ &= q_*^{\tilde{x}} (1-q_*)^{1-\tilde{x}} \end{align*}$

という成功確率 $q_* = \frac{\alpha_*}{\alpha_* + \beta_*}$ のベルヌーイ分布となる。

ガンマ分布とは

正の値のみをとる連続的な確率分布で、

$\begin{align*} p(x \mid \alpha, \theta) = \frac{x^{\alpha-1} e^{- \frac{x}{\theta} } }{\Gamma(\alpha) \theta^\alpha}, \quad x > 0, \alpha > 0, \theta > 0 \end{align*}$

という確率密度関数を持つ。

ベイズでは、 $\beta = \frac{1}{\theta}$ と変換して、

$\begin{align*} p(x \mid \alpha, \beta) = \frac{\beta^{\alpha} }{\Gamma(\alpha)} x^{\alpha-1} e^{- \beta x }, \quad x > 0, \alpha > 0, \beta > 0 \end{align*}$

を用いることが通例。

事後分布の導出

$\begin{align*} p(\lambda \mid D) &= \frac{p(D \mid \lambda) p(\lambda)}{p(D)} \\ &\propto p(D \mid \lambda) p(\lambda) \quad (\because p(D)\text{は}\lambda\text{に関して定数なので無視}) \\ &= \frac{\lambda^{\sum_{i=1}^n x_i} e^{-n \lambda}}{\prod_{i=1}^n x_i!} \frac{\beta_0^{\alpha_0} }{\Gamma(\alpha_0)} \lambda^{\alpha_0} e^{-\beta_0 \lambda} \\ &\propto \lambda^{\sum_{i=1}^n x_i} e^{-n \lambda} \times \lambda^{\alpha_0} e^{-\beta_0 \lambda} \quad (\because \lambda\text{に関して定数を無視}) \\ &= \lambda^{\sum_{i=0}^n x_i + \alpha_0 - 1} e^{-(n + \beta_0) \lambda} \end{align*}$

事後分布の導出

$\begin{align*} p(\lambda \mid D) &\propto \lambda^{\sum_{i=0}^n x_i + \alpha_0 - 1} e^{-(n + \beta_0) \lambda} \end{align*}$

の形から、

$\begin{align*} \lambda \mid D &\sim \mathcal{G}a(\alpha_*, \beta_{*}) \\ & \text{where } \alpha_{*} = \sum_{i=1}^n x_i + \alpha_0, \ \beta_{*} = n + \beta_0 \end{align*}$

であることが導出できる。

なぜか？(復習)

→　カーネルの形が同じだから

(補足) $\lambda$ の事後期待値

$\begin{align*} \text{E}[\lambda \mid x_1, \dots, x_n] &= \frac{\alpha_0+ \sum x_i}{\beta_0+n} \\ &= \frac{\beta_0}{\beta_0+n} \frac{\alpha_0}{\beta_0} + \frac{n}{\beta_0+n} \frac{\sum x_i }{n} \end{align*}$

$\beta_0$ : 事前のサンプルサイズ;
$\alpha_0$ : $\beta_0$ 個の事前の観測値の和.

十分大きな $n$ のとき、以下のように事前情報からの影響はほとんどなくなる。

$n >> \beta_0 \Rightarrow \text{E}[\lambda \mid x_1, \dots, x_n] \approx \bar{x}, \quad \text{Var}[\lambda \mid x_1, \dots, x_n] \approx \frac{\bar{x}}{n}.$

$\begin{align*} p(\tilde{x}|D) &= \int_{0}^{\infty} p(\tilde{x}|\lambda)p(\lambda|D) d\lambda \\ &= \int_{0}^{\infty} \frac{\lambda^{\tilde{x}} e^{-\lambda}}{\tilde{x}!} \frac{\beta_{*}^{\alpha_*} }{\Gamma(\alpha_*)} \lambda^{\alpha_* - 1} e^{-\beta_* \lambda} d\lambda \\ &= \frac{\beta_{*}^{\alpha_*} }{ \tilde{x}! \Gamma(\alpha_*)} \int_{0}^{\infty} \lambda^{\tilde{x} + \alpha_* - 1} e^{-(\beta_* + 1) \lambda} d\lambda \\ &= \frac{\beta_{*}^{\alpha_*} }{ \tilde{x}! \Gamma(\alpha_*)} \frac{\Gamma(\tilde{x} + \alpha_*)}{(\beta_* + 1)^{\tilde{x} + \alpha_*} } \quad (\because \text{ 菅澤先生おすすめの公式 })\\ &= \frac{\Gamma(\tilde{x} + \alpha_*)}{\tilde{x}! \Gamma(\alpha_*)} \frac{\beta_{*}^{\alpha_*} }{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \end{align*}$

ここで、ガンマ関数の性質(2) を使うと、 $\alpha_*$ が自然数である場合には

$\begin{align*} \frac{\Gamma(\tilde{x} + \alpha_*)}{\tilde{x}! \Gamma(\alpha_*)} = \frac{( \tilde{x} + \alpha_* - 1)!}{\tilde{x}! (\alpha_* - 1)!} = \begin{pmatrix} \tilde{x} + \alpha_* - 1 \\ \alpha_* - 1 \end{pmatrix} \end{align*}$

なので、

$\begin{align*} p(\tilde{x}|D) &= \begin{pmatrix} \tilde{x} + \alpha_* - 1 \\ \alpha_* - 1 \end{pmatrix} \frac{\beta_{*}^{\alpha_*} }{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \\ &= \begin{pmatrix} \tilde{x} + \alpha_* - 1 \\ \alpha_* - 1 \end{pmatrix} \left( \frac{\beta_*}{\beta_* + 1} \right)^{\alpha_*} \left( \frac{1}{\beta_* + 1} \right)^{\tilde{x} } \end{align*}$

よって、 $\tilde{x}|D \sim \text{NB}( \alpha_*, \frac{\beta_*}{\beta_* + 1})$ である。

(補足)負の 2 項分布

’成功’確率が $p$ のベルヌーイ試行について、 $r$ 回’成功’するまでに要した’失敗’の回数を $X$ とするとき、 $X$ の分布が負の 2 項分布となる。

’成功’を○、’失敗’を×で表示すると、最後は必ず○で終わるので、

××○×○× … ××○×××○

となり、合計 $r+k$ 回の試行のうち○が $r$ 回、×が $k$ 回である。

したがって確率分布は

$P(X=k \mid r,p) = \begin{pmatrix} r+k-1 \\ k \end{pmatrix} p^r q^k, \quad k=0,1,2,\ldots,$ で与えられる。この分布を負の 2 項分布(negative binomial)といい、 $NB(r,p)$ と書くことにする。

久保川(2017) p.36 を参照。

演習a 第１回発表 Python によるベイズ統計学入門2.1~3.1 2023-04-14 Kaoru Babasaki* and Keigo Takahashi* *Keio University, Graduate school of Economics