演習a 第１回発表


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


2023-04-14


Kaoru Babasaki* and Keigo Takahashi*


*Keio University, Graduate school of Economics

はじめに

You can see this presentation from here!

sugasawa-lec.kaorubb.org/presentation1/

or Read me!

qrcode.png

自己紹介 1

Keigo Takahashi

  • 所属 : 慶應義塾大学大学院経済学研究科 修士 2 年
  • 出身 : 山梨県
  • 趣味 : キャンプ・サウナ

金融を軸に,学部では機械学習寄り,修士では計量経済・ベイズ統計よりで分析してます.

最近は数理寄りのファイナンスも触れるように,ちょくちょくやってます.

2.1 未知の比率に対する推論

2 章の概要

世論調査による内閣支持率 \(=q\) の推定を題材に,ベイズの定理等のベイズ統計学における基本概念を学ぶ.

  • 事前分布・尤度・事後分布 → ベイズの定理
  • 点推定 (損失関数)
  • 信用区間（→区間推定）
  • ベイズファクター・ SDDR（→ 仮説検定）
  • 予測分布

本発表では,以上の概念を数式を交えてできるだけ丁寧にさらっていく.

q の解釈

事前分布の話は次の節 2.2 でまとめて話す. ここでは推定対象である\(q\)について整理する.

確率\(q\)は,本質的には世論調査に応じた人が「支持する」と回答する確率である.

\begin{align*} Pr\{X_i = x_i \} = \begin{cases} q & (x_i = 1)\\ 1-q & (x_i = 0) \end{cases} \end{align*}

q の解釈-1

この \(q\) が、推測したい内閣支持率そのものになるのは,

  1. 抽出対象のサンプル数が極端に大きい
  2. 無作為抽出

を前提しているからである.

\(q\)は内閣を支持する人を抽出できたという意味を含んで「成功確率」と言える.

2.2 ベイズの定理による事後分布の導出

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

  • 真の値が未知であるパラメータ\(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(q)\)

※ 2.1 未知の比率に対する推論での内容

ベイズによる推論の出発点.

パラメータのとりうる値に対して分析者が抱いている確信の度合い・情報・バイアスを確率分布として表現したものである.

本章では一様分布とベータ分布の 2 つが事前分布が紹介されている.

一様分布

\begin{align*} p(q) = \begin{cases} \frac{1}{b-a} & (a \leq q \leq b) \\ 0 & (q < a , b < q) \end{cases} \end{align*}

uniform_prior.png

  • これは無情報事前分布である. パラメータに対して特に見立て等がない場合はこれを採用する.
  • 他の無情報事前分布には Jeffereys の事前分布がある.

ベータ分布

\begin{align*} p(q) & = \frac{q^{\alpha-1}(1-q)^{\beta-1}}{B(\alpha,\beta)} \end{align*}

※\( B(\alpha,\beta) = \int_{0}^{1} q^{\alpha-1}(1-q)^{\beta-1} dq \)

ベータ分布は \(\alpha > 1\),\(\beta > 1\) の範囲のベータ分布が使われることが多い.パラメータ\(q\) に対する分析者の予想を最頻値として反映できるため.

beta_prior_2.png

尤度\(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\)の値に依らない

補足:周辺尤度\(p(D)\)

事後分布の導出の際,周辺尤度はデータが観測された後固定化されるため,比例記号を使う際省かれている.

先程の説明にある「パラメータの値に寄らない」と同意だが,ここのイメージを尤度と共にもう少し具体的にする.

事後分布\(p(D|q)\)

事前分布\(p(q)\)と尤度\(p(D|q)\)を使って, パラメータ\(=q\) の事後分布を導く.

\(p(D)\)は,パラメータの値に依らず計算されるため, パラメータ の事後分布は 尤度 と 事前分布 に比例する.

\[ p(q|D) \propto p(D|q)p(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\]

自然共役事前分布

  • ベルヌーイ分布に,ベータ分布を想定したら,事後分布もベータ分布になった.
  • ある特定の尤度に,ある特定の事前分布を想定したら,事後分布も同じ分布になる.

この様な事前分布は,自然共役事前分布と呼ばれる.

ベルヌーイ分布に対してはベータ分布,ポアソン分布に対してはガンマ分布.

メリット

  1. 事後分布の形状がわかっているので,パラメータのみの演算だけとなり計算が楽.
  2. 新しい事後分布の導出をするとき,古い事後分布を事前分布に使用できる.

1 が画期的であり,積分計算とか不要. 2 については,カルマンフィルタでその醍醐味を知ることができる(勉強中)

2.3 未知のパラメータに関する議論

概要

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

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

  • 点推定
  • 区間推定
  • 仮説検定

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

点推定

分布の代表値は平均値・中央値・最頻値の 3 つ.

ベイズ統計学における点推定の値 \(\delta\) も,以上 3 つの代表値を用いる.

どの代表値を用いるかは,損失関数 \((L(q,\delta))\) による.

パラメータの真の値との損失・誤差を小さくする 代表値を点推定の値 \((=\delta)\) とする.

代表的な損失関数\(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)区間

信用区間

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

このとき, \(100(1-c)\) 信用区間 \([a_c,b_c]\) は

\begin{align*} Pr \{ q < a_c \mid D \} = \frac{c}{2} , Pr \{ q > b_c \mid D \} = \frac{c}{2} \end{align*}

を満たす \([a_c,b_c]\)

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 に近い値なのであればその仮説は成り立たないと判断.

複数の仮説\(H_i\)の比較

確率が 1 or 0 に近いといっても,どれくらい近いと仮説が支持・棄却されるかはまだ少し曖昧.

複数の仮説\(H_i\)を 事後オッズ比 もしくは ベイズファクター で比較・検証する.

※ \(H_i\)の範囲\(S_i\)は,通常\(\cap S_i = \phi\) かつ \(\cup S_i=\)パラメータの定義域 の必要がある.

事後オッズ比・ベイズファクター\(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)\)を使用.

以上から,それぞれの仮説の事後確率

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

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} \)

自己紹介 2

Kaoru Babasaki

  • 出身：川崎
  • 所属: 慶應義塾大学経済学研究科修士課程 2 年
  • 趣味:
    • サーフィン, スノーボード, スケートボード
    • Linux, emacs, キーボード
  • Kevin(8 歳) kevin.png

2.4 将来の確率変数の値の予測

目的

いま手元にあるデータから、将来観測されるであろう予測することを考える。

例

  • 内閣支持率
  • 製品の売上
  • 病気の発症率

この節のポイント

データを\(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)} \]

導出

\begin{align*} p(\tilde{x} \mid D) &= \frac{p(\tilde{x}, D)}{p(D)} \\ &= \frac{ \int_0^1 p(\tilde{x}, D \mid q) p(q) dq }{\int_0^1 p(D \mid q) p(q) dq} \\ &= \int_0^1 p(\tilde{x} \mid q) \frac{p(D \mid q) p(q)}{\int_0^1 p(D \mid q) p(q) dq} dq \quad (\because \tilde{x} \text{と} D \text{は(qの条件付き)独立} ) \\ &= \int_0^1 p(\tilde{x} \mid q) p(q \mid D) dq \end{align*}

(補足)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*}

ベータ関数の性質

ベータ関数の定義
\begin{align*} B(\alpha, \beta) = \int_0^1 x^{\alpha - 1} (1-x)^{\beta - 1} dx = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\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_*}\) のベルヌーイ分布となる。

Julia で実演

3.1 ポアソン分布のベイズ分析

問題設定

分析したいデータが、まれに起こるイベントの発生回数

  • 小さな美容院への１日の来店数
  • 特定の地区で１日に発生した犯罪や交通事故の件数

のような、非負の整数データをモデリングする際、ポアソン分布がよく用いられる。

ポアソン分布とは

確率関数が \[ p(x \mid \lambda) = \frac{\lambda^x e^{-\lambda}}{x!}, \quad x = 0, 1, 2, \cdots, \lambda > 0 \] で与えられる確率分布のことを、パラメータ\(\lambda\)を持つポアソン分布という。

mean-variance relationship

\(x \sim \) Poisson(\(\lambda\)) ならば、

\begin{align*} E[x \mid \lambda] = Var[x \mid \lambda] = \lambda \end{align*}

ポアソン分布のプロット

poisson.gif

事前分布の設定

ここでは、sampling model である ポアソン分布のパラメータ\(\lambda\)に対して、事前分布として、ガンマ分布を仮定する。

\begin{align*} \lambda \sim \text{Gamma}(\alpha_0, \beta_0) \end{align*}

ガンマ分布とは

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

\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*}

を用いることが通例。

ガンマ分布の性質

\[ \lambda \sim \text{Gamma}(a,b) \] のとき、

\[\begin{align*} \text{E}(\lambda) &= a/b \\ \text{Var}(\lambda) &= a/b^2. \\ \text{mode}[\lambda] &= \begin{cases} \frac{a-1}{b} & \text{if } a > 1 \\ 0 & \text{if } a \le 1 \end{cases} \end{align*}\]

覚えたほうがいい公式(by 菅澤先生)

\begin{align*} & \quad\int_0^{\infty} f(x) dx = 1 \\ &\Leftrightarrow \int_0^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha - 1} e^{- \beta x} dx = 1 \\ &\Leftrightarrow \int_0^{\infty} x^{\alpha - 1} e^{- \beta x} dx = \frac{\Gamma(\alpha)}{\beta^{\alpha} } \\ \end{align*}

ガンマ分布のプロット

gamma.gif

いざ、ベイズ

\begin{align*} x_1, \cdots, x_n &\overset{\text{i.i.d.} }{\sim} \text{Poisson}(\lambda) \\ \lambda &\sim \text{Gamma}(\alpha_0, \beta_0) \end{align*}

これでモデルの仮定が整ったので、あとは、ベイズの定理に従って事後分布を導出するだけ。

尤度

\begin{align*} p(D \mid \lambda) &= p(x_1, \cdots, x_n \mid \lambda) \\ &= \prod_{i=1}^n p(x_i \mid \lambda) \quad \because \text{i.i.d.の仮定} \\ &= \prod_{i=1}^n \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} \\ &= \frac{\lambda^{\sum_{i=1}^n x_i} e^{-n \lambda}}{\prod_{i=1}^n x_i!} \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*} 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}. \]

事後予測分布の導出

\[ p(\tilde{x}|D) = \int_{\theta \in \Theta} p(\tilde{x}|\theta)p(\theta|D) d\theta \] を思い出すと、

\begin{align*} p(\tilde{x}|D) &= \int_{0}^{\infty} p(\tilde{x}|\lambda)p(\lambda|D) d\lambda \\ \end{align*}

なので、

(続き) botsu

\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{1}{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \int_{0}^{\infty} (\beta_* + 1) \cdot [(\beta_* + 1) \lambda]^{\tilde{x} + \alpha_* - 1} e^{-(\beta_* + 1) \lambda} d\lambda \\ &= \frac{\beta_{*}^{\alpha_*} }{ \tilde{x}! \Gamma(\alpha_*)} \frac{1}{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \int_{0}^{\infty} [(\beta_* + 1) \lambda]^{\tilde{x} + \alpha_* - 1} e^{-(\beta_* + 1) \lambda} \left(\frac{d}{d \lambda}(\beta_* + 1) \lambda \right) d \lambda \\ &= \frac{\beta_{*}^{\alpha_*} }{ \tilde{x}! \Gamma(\alpha_*)} \frac{1}{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \int_{0}^{\infty} [(\beta_* + 1) \lambda]^{\tilde{x} + \alpha_* - 1} e^{-(\beta_* + 1) \lambda} d((\beta_* + 1) \lambda ) \quad (\because \text{ 置換積分}) \\ &= \frac{\Gamma(\tilde{x} + \alpha_*)}{\tilde{x}! \Gamma(\alpha_*)} \frac{\beta_{*}^{\alpha_*} }{(\beta_* + 1)^{\tilde{x} + \alpha_*}} \end{align*}
\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 を参照。

例

Julia で実演

おわり

ご清聴ありがとうございました。

質問などありましたらどうぞ。