演習b 第1回発表


標準ベイズ統計学8.1~8.3
グループ比較と階層モデリング(前半)


2023-10-06


Kaoru Babasaki*


*Keio University, Graduate school of Economics

8.1 二つのグループを比較する

米国の二つの高校の数学の得点データ



  学校 1 学校 2
生徒数 約 600 人 約 600 人
標本サイズ 31 人 28 人
平均点 50.81 46.15


  • 二つの高校の得点に差はあると言えるか?

fig_box.png

頻度論的アプローチ

\(t\) 検定

  • 正規分布に従う 2 標本

    \begin{align*} y_{1,1}, \dots, y_{n_1, 1}, i.i.d. &\sim N(\theta_1, \sigma^2), \\ y_{1,2}, \dots, y_{n_2, 2}, i.i.d. &\sim N(\theta_2, \sigma^2) \end{align*}

    を考える。

  • 2 つの標本の母平均\(\theta_1, \theta_2\)が等しいか否かを検定
  • よって

    • 帰無仮説: \(H_0: \theta_1 = \theta_2\)
    • 対立仮説: \(H_1: \theta_1 \neq \theta_2\)

    とし、仮設検定を行う。

帰無仮説を仮定すると…

  • 母分散\(\sigma^2\)のプールされた推定量 \[ s_p^2 = \frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2} \] を用いると、 t 統計量 \[ t(\boldsymbol{y}_1, \boldsymbol{y}_2) = \frac{\bar{y}_1 - \bar{y}_2}{s_p \sqrt{\frac{1}{n_1} + \frac{1}{n_2} } } \] は、自由度\(n_1 + n_2 - 2\)の t 分布に従う。
  • 証明は久保川 (2017) を参照。
  • 有意水準を\(1-\alpha\)として検定を行う。
  • \(t_{n_1+n_2-2, \frac{\alpha}{2} }\)を \(P\left( T > t_{n_1+n_2-2, \frac{\alpha}{2} } \right) = \frac{\alpha}{2}\)となる点とすると、

    \begin{align*} &P_{H_0}(t(\boldsymbol{y}_1, \boldsymbol{y}_2) > t_{n_1+n_2-2, \frac{\alpha}{2} }) +P_{H_0}(-t(\boldsymbol{y}_1, \boldsymbol{y}_2) > t_{n_1+n_2-2, \frac{\alpha}{2} }) \\ = &P_{H_0}(|t(\boldsymbol{y}_1, \boldsymbol{y}_2)| > t_{n_1+n_2-2, \frac{\alpha}{2} }) \\ = &\alpha \end{align*}
  • よって\(|t(\boldsymbol{y}_1, \boldsymbol{y}_2)|\)が \(t_{n_1+n_2-2, \frac{\alpha}{2} }\)より大きい場合、すなわち p 値が\(\alpha\)より小さい場合、帰無仮説\(H_0\)を棄却する。

p 値に基づくモデル選択

  • 有意水準を 95% \(=1-0.05\)として検定を行う。
  • \(p < 0.05\)の場合
    • 2 つのグループが同じ分布をもつ(=母平均が等しい)というモデル(帰無仮説)を棄却する。
    • \(\theta_1 \neq \theta_2\)と結論づける。
    • 推定値として標本平均、すなわち \(\hat{\theta}_1=\bar{y}_1, \hat{\theta}_2=\bar{y}_2\) を用いる。
  • \(p > 0.05\)の場合
    • 帰無仮説は棄却できず、 2 つのグループが同じ分布をもつというモデルを採用。
    • \(\theta_1 = \theta_2\)と結論づける。
    • 推定値として、プールされた平均値、すなわち \(\hat{\theta}_1 = \hat{\theta}_2 = \frac{ \sum y_{i,1} + \sum y_{i,2} }{n_1 + n_2}\) を用いる。

二つの高校の数学の得点データの t 検定によるモデル選択

p 値に基づくモデル選択の問題点

  • 二項対立的な考え方のデメリット
    • 検定を行ったデータはランダムなサンプルのため、学校 1 のサンプルには成績の良い生徒がたまたま多く含まれていたり、学校 2 のサンプルには成績の悪い生徒がたまたま多く含まれていたりする可能性がある。
    • そのため、再びランダムサンプリングし同じプロセスを実行すると、真逆の結果になることも。
  • 各母集団を別々に推定する方法は?
    • 同じような母集団であることはある程度わかっているので、やや非効率。
    • いずれかのサンプルサイズが小さい場合、そのグループの推定値が不安定になる。

p 値に基づくモデル選択の別の解釈

  • 学校 1 の母平均\(\theta_1\)について

    • \(p < 0.05\)の場合、 \(\hat{\theta}_1=\bar{y}_1\)
    • \(p > 0.05\)の場合、 \(\hat{\theta}_1=\frac{ \sum y_{i,1} + \sum y_{i,2} }{n_1 + n_2}\)

    であった。

  • これらは以下のように解釈し直せる。

    \begin{align*} \hat{\theta}_1 &= w \bar{y}_1 + (1-w) \bar{y}_2 \quad \text{ where } w = \begin{cases} 1 & \text{ if } p < 0.05 \\ \frac{n_1}{n_1 + n_2} & \text{ otherwise } \end{cases} \end{align*}
  • このような見方だと、かなり極端な方法であることがわかる。
  • 重みを連続的に変化させる方法は? → ベイズ的アプローチ

ベイズ的アプローチ

サンプリングモデルの設定

\[\begin{aligned} Y_{i,1} &= \mu + \delta + \epsilon_{i,1} \\ Y_{i,2} &= \mu - \delta + \epsilon_{i,2} \\ \{\epsilon_{i,j}\} &\sim \text{i.i.d. normal}(0, \sigma^2) \end{aligned}\]

  • \(\theta_1 = \mu + \delta, \theta_2 = \mu - \delta\)とすると、

    • \(\delta = \frac{\theta_1 - \theta_2}{2}\) … 母集団の平均値の差の半分
    • \(\mu = \frac{\theta_1 + \theta_2}{2}\) … プールされた平均値

    と解釈できる。

事前分布の設定

\[\begin{aligned} p(\mu, \delta, \sigma^2) &= p(\mu) \times p(\delta) \times p(\sigma^2) \\ \mu &\sim \text{normal}(\mu_0, \gamma_0^2) \\ \delta &\sim \text{normal}(\delta_0, \tau_0^2) \\ \sigma^2 &\sim \text{inverse-gamma}( \frac{\nu_0}{2}, \frac{\nu_0 \sigma_0^2}{2} ) \end{aligned}\]

完全条件付き分布の計算

full conditional of \(\mu\)
\begin{align*} \{\mu | &\boldsymbol{y}_1, \boldsymbol{y}_2, \delta, \sigma^2\} \sim \text{normal}(\mu_n, \gamma_n^2), \text{ where} \\ & \mu_n = \gamma_n^2 \times \left[ \frac{\mu_0}{\gamma_0^2} + \frac{\sum_{i=1}^{n_1} (y_{i,1} - \delta)}{\sigma^2} + \frac{\sum_{i=1}^{n_2} (y_{i,2} + \delta)}{\sigma^2} \right] \\ & \gamma_n^2 = \left[ \frac{1}{\gamma_0^2} + \frac{n_1 + n_2}{\sigma^2} \right]^{-1} \end{align*}
  • 導出

    \[\begin{aligned} p(\mu | \boldsymbol{y}_1, \boldsymbol{y}_2, \delta, \sigma^2) &\propto p(\boldsymbol{y}_1, \boldsymbol{y}_2, \delta, \sigma^2 | \mu) \times p(\mu) \\ &= p(\boldsymbol{y}_1, \boldsymbol{y}_2| \delta, \sigma^2, \mu) \times p(\delta) p(\sigma^2) p(\mu) \\ &\propto p(\boldsymbol{y}_1| \delta, \sigma^2, \mu) \times p(\boldsymbol{y}_2| \delta, \sigma^2, \mu) \times p(\mu) \\ &= \prod_{i=1}^{n_1} p(y_{i,1}| \delta, \sigma^2, \mu) \times \prod_{i=1}^{n_2} p(y_{i,2}| \delta, \sigma^2, \mu) \times p(\mu) \\ &\propto \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_1} (y_{i,1} - \mu - \delta)^2 \right\} \times \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_2} (y_{i,2} - \mu + \delta)^2 \right\} \times \exp \left\{ -\frac{1}{2\gamma_0^2} (\mu - \mu_0)^2 \right\} \\ &= \exp \left\{ -\frac{1}{2\sigma^2} \left[ \sum_{i=1}^{n_1} (\mu^2 - 2(y_{i,1} - \delta)\mu + (y_{i,1} - \delta)^2) + \sum_{i=1}^{n_2} (\mu^2 - 2(y_{i,2} + \delta)\mu + (y_{i,2} + \delta)^2) \right] -\frac{1}{2\gamma_0^2} (\mu^2 - 2\mu_0\mu + \mu_0^2) \right\} \\ &\propto \exp \left\{ -\frac{1}{2\sigma^2} \left[ n_1 \mu^2 - 2\mu \sum_{i=1}^{n_1} (y_{i,1} - \delta) + n_2 \mu^2 - 2\mu \sum_{i=1}^{n_2} (y_{i,2} + \delta) \right] -\frac{1}{2\gamma_0^2} (\mu^2 - 2\mu_0\mu) \right\} \\ &= \exp \left\{ - \frac{1}{2} \left(\frac{n_1 + n_2}{\sigma^2} + \frac{1}{\gamma_0^2} \right) \mu^2 + \left(\frac{\sum_{i=1}^{n_1} (y_{i,1} - \delta) + \sum_{i=1}^{n_2} (y_{i,2} + \delta)}{\sigma^2} + \frac{\mu_0}{\gamma_0^2} \right) \mu \right\} \\ &= \exp \left\{ - \frac{1}{2 \gamma_n^2} \mu^2 + \frac{\mu_n}{\gamma_n^2} \mu \right\} \\ &\propto \exp \left\{ - \frac{1}{2 \gamma_n^2} \left(\mu^2 - 2\mu \mu_n + \mu_n^2 \right) \right\} \\ &= \exp \left\{ - \frac{1}{2 \gamma_n^2} (\mu - \mu_n)^2 \right\} \\ &\propto \text{dnorm}(\mu | \mu_n, \gamma_n) \end{aligned}\]

Full conditional of \(\delta\)

\[\begin{aligned} \{\delta | &\boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \sigma^2\} \sim \text{normal}(\delta_n, \tau_n^2), \text{ where} \\ & \delta_n = \tau_n^2 \times \left[ \frac{\delta_0}{\tau_0^2} + \frac{\sum_{i=1}^{n_1} (y_{i,1} - \mu)}{\sigma^2} -\frac{\sum_{i=1}^{n_2} (y_{i,2} - \mu)}{\sigma^2} \right] \\ & \tau_n^2 = \left[ \frac{1}{\tau_0^2} + \frac{n_1 + n_2}{\sigma^2} \right]^{-1} \end{aligned}\]

  • 導出

    \[\begin{aligned} p(\delta | \boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \sigma^2) &\propto p(\boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \sigma^2 | \delta) \times p(\delta) \\ &\propto p(\boldsymbol{y}_1, \boldsymbol{y}_2| \mu, \sigma^2, \delta) \times p(\mu) p(\sigma^2) p(\delta) \\ &\propto p(\boldsymbol{y}_1| \mu, \sigma^2, \delta) \times p(\boldsymbol{y}_2| \mu, \sigma^2, \delta) \times p(\delta) \\ &= \prod_{i=1}^{n_1} p(y_{i,1}| \mu, \sigma^2, \delta) \times \prod_{i=1}^{n_2} p(y_{i,2}| \mu, \sigma^2, \delta) \times p(\delta) \\ &\propto \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_1} (y_{i,1} - \mu - \delta)^2 \right\} \times \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_2} (y_{i,2} - \mu + \delta)^2 \right\} \times \exp \left\{ -\frac{1}{2\tau_0^2} (\delta - \delta_0)^2 \right\} \\ &= \exp \left\{ -\frac{1}{2\sigma^2} \left[ \sum_{i=1}^{n_1} (\delta^2 - 2(y_{i,1} - \mu)\delta + (y_{i,1} - \mu)^2) + \sum_{i=1}^{n_2} (\delta^2 + 2(y_{i,2} - \mu)\delta + (y_{i,2} - \mu)^2) \right] -\frac{1}{2\tau_0^2} (\delta^2 - 2\delta_0\delta + \delta_0^2) \right\} \\ &\propto \exp \left\{ -\frac{1}{2\sigma^2} \left[ n_1 \delta^2 - 2\delta \sum_{i=1}^{n_1} (y_{i,1} - \mu) + n_2 \delta^2 + 2\delta \sum_{i=1}^{n_2} (y_{i,2} - \mu) \right] -\frac{1}{2\tau_0^2} (\delta^2 - 2\delta_0\delta) \right\} \\ &= \exp \left\{ - \frac{1}{2} \left(\frac{n_1 + n_2}{\sigma^2} + \frac{1}{\tau_0^2} \right) \delta^2 + \left(\frac{\sum_{i=1}^{n_1} (y_{i,1} - \mu) - \sum_{i=1}^{n_2} (y_{i,2} - \mu)}{\sigma^2} + \frac{\delta_0}{\tau_0^2} \right) \delta \right\} \\ &= \exp \left\{ - \frac{1}{2 \tau_n^2} \delta^2 + \frac{\delta_n}{\tau_n^2} \delta \right\} \\ &\propto \exp \left\{ - \frac{1}{2 \tau_n^2} \left(\delta^2 - 2\delta \delta_n + \delta_n^2 \right) \right\} \\ &= \exp \left\{ - \frac{1}{2 \tau_n^2} (\delta - \delta_n)^2 \right\} \\ &\propto \text{dnorm}(\delta | \delta_n, \tau_n) \end{aligned}\]

Full conditional of\(\sigma^2\)

\[\begin{aligned} \{\sigma^2 | &\boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \delta\} \sim \text{inverse-gamma}( \frac{\nu_n}{2}, \frac{\nu_n \sigma_n^2}{2}), \text{ where} \\ & \nu_n = \nu_0 + n_1 + n_2 \\ & \nu_n \sigma_n^2 = \nu_0 \sigma_0^2 + \sum_{i=1}^{n_1} (y_{i,1} - [\mu + \delta])^2 + \sum_{i=1}^{n_2} (y_{i,2} - [\mu - \delta])^2 \end{aligned}\]

  • 導出

    \[\begin{aligned} & \quad p(\sigma^2 | \boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \delta) \\ &\propto p(\boldsymbol{y}_1, \boldsymbol{y}_2, \mu, \delta | \sigma^2) \times p(\sigma^2) \\ &= p(\boldsymbol{y}_1, \boldsymbol{y}_2| \mu, \delta, \sigma^2) \times p(\mu) p(\delta) p(\sigma^2) \\ &\propto p(\boldsymbol{y}_1| \mu, \delta, \sigma^2) \times p(\boldsymbol{y}_2| \mu, \delta, \sigma^2) \times p(\sigma^2) \\ &= \prod_{i=1}^{n_1} p(y_{i,1}| \mu, \delta, \sigma^2) \times \prod_{i=1}^{n_2} p(y_{i,2}| \mu, \delta, \sigma^2) \times p(\sigma^2) \\ &\propto (\sigma^2)^{- \frac{n_1}{2} } \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_1} (y_{i,1} - [\mu + \delta])^2 \right\} \times (\sigma^2)^{- \frac{n_2}{2} } \exp \left\{ -\frac{1}{2\sigma^2} \sum_{i=1}^{n_2} (y_{i,2} - [\mu - \delta])^2 \right\} \times (\sigma^2)^{-\frac{\nu_0}{2} - 1} \exp \left\{ -\frac{\nu_0 \sigma_0^2}{2\sigma^2} \right\} \\ &= (\sigma^2)^{- \frac{n_1 + n_2 + \nu_0}{2} - 1} \exp \left\{ -\frac{1}{2\sigma^2} \left[ \sum_{i=1}^{n_1} (y_{i,1} - [\mu + \delta])^2 + \sum_{i=1}^{n_2} (y_{i,2} - [\mu - \delta])^2 + \nu_0 \sigma_0^2 \right] \right\} \\ &= (\sigma^2)^{- \frac{\nu_n}{2} - 1} \exp \left\{ -\frac{\nu_n \sigma_n^2 }{2\sigma^2} \right\} \\ &\propto \text{dinverse-gamma}(\sigma^2 | \frac{\nu_n}{2}, \frac{\nu_n \sigma_n^2}{2}) \end{aligned}\]

数学試験データの分析

(おまけ) 数値実験

8.2 複数のグループを比較する

階層データ

  • 前節で見たような、入れ子になった集団の構造をもつデータは、階層データ (hierarchical data) やマルチレベルデータ (multilevel data) と呼ばれる。
    • 複数の病院に入院している患者のデータ
    • 動物のあるグループ内の遺伝子
    • ある国の地域の群に属する人々

8.2.1 交換可能性と階層モデル

(復習) 交換可能性 (exchangeability)

  • \(p(y_1, \dots, y_n)\)を\(Y_1, \dots, Y_n\)の同時密度とする。もし、\(\{1, \dots, n\}\)の任意の置換\(\pi\)に対して \[ p(y_1, \dots, y_n) = p(y_{\pi(1)}, \dots, y_{\pi(n)}) \] が成り立つならば、\(Y_1, \dots, Y_n\)は交換可能であるという。
  • 添字の付け方が出力についての情報を持たないということ。
  • 例:
    • 交換可能: コインを 5 回投げて表が 3 回出る確率 \[p(1,0,0,1,1) = p(0,1,1,0,1) = p(1,1,0,0,1) = \dots\]
    • 交換可能でない: 大学 1,2,3,4 年次それぞれの累積 GPA

      \begin{align*} p(GPA_1, GPA_2, GPA_3, GPA_4) &\neq p(GPA_4, GPA_1, GPA_3, GPA_2) \\ p(3.5, 3.0, 2.5, 2.0) &\neq p(2.0, 3.5, 2.5, 3.0) \end{align*}

(復習) デ・フィネッティの定理

任意の\(i \in \{1,2, \dots\}\)について\(Y_i \in \mathcal{Y}\)とする。任意の\(n\)に対して、 \(Y_1, \dots, Y_n\)に対する信念のモデルは \(\{1, \dots, n\}\)の任意の置換\(\pi\)に対して、交換可能であるとする。つまり、 \[ p(y_1, \dots, y_n) = p(y_{\pi(1)}, \dots, y_{\pi(n)}) \] が成り立っているとする。このとき、パラメータ\(\phi\)と、\(\phi\)に関する事前分布\(p(\phi)\)と、標本モデル\(p(y|\phi)\)が存在して、モデルは次のように書ける。 \[ p(y_1, \dots, y_n) = \int \left\{ \prod_1^n p(y_i|\phi) \right\} p(\phi) d\phi \]

  • 要するに、
  • データ\(Y_1, \dots, Y_n\)に交換可能性を仮定すると、

    • サンプリングモデルの関数系\(p(y|\phi)\)と、
    • 事前分布\(p(\phi)\)

    が存在し、以下のような定式化が可能である。

\begin{align*} \phi &\sim p(\phi) \\ \{Y_1, \dots, Y_n\} &\sim \text{ i.i.d. } p(y|\phi) \end{align*}

階層データに対するデ・フィネッティの定理の適用

グループ\(j\)のデータに対する適用
  • グループ\(j\)のデータを\(\boldsymbol{Y}_j = \{Y_{1,j}, \dots, Y_{n_j, j}\}\)とし、データ全体を\(\{\boldsymbol{Y}_1, \dots, \boldsymbol{Y}_m\}\)とする。
  • グループ\(j\)のデータは交換可能であるとすると(妥当な仮定)、あるパラメータ\(\phi_j\)が存在して、

    \begin{equation} \tag{1} \begin{aligned} \phi_j &\sim p(\phi_j) \\ \{Y_{1,j}, \dots, Y_{n_j, j}\} &\sim \text{ i.i.d. } p(y|\phi_j) \end{aligned} \end{equation}

    と定式化できる。

グループ固有のパラメータ\(\phi_j\)に対する適用
  • (1)の定式化のみで推論を行ってよいか?
    • \(\phi_1, \dots, \phi_m\)が独立を暗に仮定したことになる。
    • では、\(\phi_1, \dots, \phi_m\)が独立であるか?
  • もし独立であれば、 \(\phi_1, \dots, \phi_{m-1}\)の値を知っていても \(\phi_m\)に関する情報は変わらない。→ 不自然!
  • よって、独立でない、つまり、 \(p(\phi_1, \dots, \phi_m) \neq \prod_1^m p(\phi_j)\)とする。
  • 一方、グループ全体がより大きなグループの母集団からのサンプルである場合、グループ固有のパラメータ\(\phi_j, j=1, \dots, m\)の交換可能性は妥当。
  • したがって、

    \begin{align*} \{ \phi_1, \dots, \phi_m \} &\sim \text{ i.i.d. } p(\phi|\psi) \\ \psi &\sim p(\psi) \end{align*}

    と定式化できる。

まとめると
\begin{align*} \{y_{1,j}, \dots, y_{n_j, j} | \phi_j \} &\sim \text{ i.i.d. } p(y|\phi_j)& \quad & \text{(グループ内の変動)} \\ \{ \phi_1, \dots, \phi_m | \psi \} &\sim \text{ i.i.d. } p(\phi|\psi) & \quad & \text{(グループ間の変動)} \\ \psi &\sim p(\psi) & \quad & \text{(事前分布)} \end{align*}
  • このように定式化することで、階層構造を考慮したモデルを構築できる。
  • 階層モデルのメリット:
    • 効率的: \(\theta_j\)の推定において、グループ全てに共通するハイパーパラメータ\(\psi\)を導入することで、間接的に別のグループの情報も利用できている。
    • overfitting の抑制: データ数が少ないグループにおいて、そのグループのデータのみを用いて推定すると、オーバーフィッティングが起こりやすい。しかし、階層モデルは、他グループのデータも利用することで、オーバーフィッティングを抑制する。

8.3 階層正規モデル

階層正規モデル: 定式化

\begin{align*} \{y_{1,j}, \dots, y_{n_j, j}\} &\sim \text{i.i.d.} \mathcal{N}(\theta_j, \sigma^2) \text{ for } j=1, \dots, m \\ \theta_j &\sim \mathcal{N}(\mu, \tau^2) \end{align*}

fig8_3.png

  • ここではまだグループごとの分散の不均一性は考慮せず、均一分散を仮定していることに注意。

事前分布の設定

\begin{align*} 1/\sigma^2 &\sim \text{gamma}(\nu_0/2, \nu_0\sigma_0^2/2) \\ 1/\tau^2 &\sim \text{gamma}(\nu_0/2, \nu_0\tau_0^2/2) \\ \mu &\sim \text{normal}(\mu_0, \gamma_0^2). \end{align*}

8.3.1 事後推論

鍵となる分解

\[\begin{aligned} &p(\theta_1, \dots, \theta_m, \mu, \tau^2, \sigma^2 | \boldsymbol{y}_1, \dots, \boldsymbol{y}_m) \\ \propto &p(\theta_1, \dots, \theta_m, \mu, \tau^2, \sigma^2) \times p(\boldsymbol{y}_1, \dots, \boldsymbol{y}_m | \theta_1, \dots, \theta_m, \mu, \tau^2, \sigma^2) \\ \propto &p(\mu, \tau^2, \sigma^2) \times p(\theta_1, \dots, \theta_m | \mu, \tau^2, \sigma^2) \times p(\boldsymbol{y}_1, \dots, \boldsymbol{y}_m | \theta_1, \dots, \theta_m, \mu, \tau^2, \sigma^2) \\ \propto &p(\mu) p(\tau^2) p(\sigma^2) \left\{ \prod_{j=1}^m p(\theta_j | \mu, \tau^2) \right\} \left\{ \prod_{j=1}^m \prod_{i=1}^{n_j} p(y_{i,j} | \theta_j, \sigma^2) \right\}. \\ \end{aligned}\]

Full conditional distributions of \(\mu\) and \(\tau^2\)

\[\begin{aligned} \{\mu | \theta_1, \dots, \theta_m, \tau^2\} &\sim \text{normal}(\frac{ m\bar{\theta}/\tau^2 + \mu_0/\gamma_0^2 }{ m/\tau^2 + 1/\gamma_0^2 }, [m/\tau^2 + 1/\gamma_0^2]^{-1}) \\ \{1/\tau^2 | \theta_1, \dots, \theta_m, \mu\} &\sim \text{gamma}(\frac{\eta_0 + m}{2}, \frac{\eta_0 \tau_0^2 + \sum_{j=1}^m (\theta_j - \mu)^2}{2}). \\ \end{aligned}\]

導出: \(\mu\)

\(\mu\)に関しては、section 5.2 と全く同じ方法で導出できる。 \[\begin{aligned} y \sim p(\theta, \sigma^2) &\rightarrow \theta \sim p(\mu, \tau^2) \\ \theta \sim p(\mu_0, \tau_0^2) &\rightarrow \mu \sim p(\mu_0, \gamma_0^2) \\ \sigma^2 &\rightarrow \tau^2 \\ \end{aligned}\] とおきかえるだけ。

\begin{align*} p(\mu | \theta_1, \dots, \theta_m, \tau^2) &\propto p(\mu) \left\{ \prod_{j=1}^m p(\theta_j | \mu, \tau^2) \right\} \\ &\propto \exp \left\{ -\frac{1}{2\gamma_0^2} (\mu - \mu_0)^2 \right\} \exp \left\{ -\frac{1}{2\tau^2} \sum_{j=1}^m (\theta_j - \mu)^2 \right\} \\ \end{align*}

ここで、指数部分は\(-\frac{1}{2}\)を一旦無視して、

\begin{align*} \frac{1}{\gamma_0^2}( \mu^2 - 2\mu\mu_0 + \mu_0^2 ) +\frac{1}{\tau^2}( \sum_{j=1}^m \theta_j^2 - 2\mu\sum_{j=1}^m \theta_j + m\mu^2 ) = a\mu^2 - 2b\mu + c \\ a = \frac{1}{\gamma_0^2} + \frac{m}{\tau^2}, \quad b = \frac{\mu_0}{\gamma_0^2} + \frac{\sum_{j=1}^m \theta_j}{\tau^2}, \quad c = c(\theta_1, \dots, \theta_m, \mu_0, \gamma_0^2, \tau^2) \\ \end{align*}

と置けるので、

\begin{align*} p(\mu | \theta_1, \dots, \theta_m, \tau^2) &\propto \exp \left\{ -\frac{1}{2} (a\mu^2 - 2b\mu ) \right\} \\ &= \exp \left\{ -\frac{1}{2} a(\mu^2 - 2\frac{b \mu}{a} + \frac{b^2}{a^2}) + \frac{b^2}{2a} \right\} \\ &\propto \exp \left\{ -\frac{1}{2} a(\mu - \frac{b}{a})^2 \right\} \\ &= \exp \left\{ -\frac{1}{2} \left( \frac{\mu - b/a}{1/\sqrt{a} } \right)^2 \right\} \\ &\propto \text{dnorm}(\mu | \frac{b}{a}, \frac{1}{\sqrt{a} }) \\ \end{align*}
導出: \(\tau^2\)

\(1/\tau^2\)に関しては、section 6.3 と全く同じ方法で導出できる。

\[\begin{aligned} \theta &\rightarrow y \\ \mu &\rightarrow \theta \\ 1/\tau^2 &\rightarrow \tilde{\sigma}^2 \\ \end{aligned}\]

とおきかえるだけ。

\begin{align*} p(1/\tau^2 | \theta_1, \dots, \theta_m, \mu) &\propto p(\tau^2) \left\{ \prod_{j=1}^m p(\theta_j | \mu, \tau^2) \right\} \\ &\propto \left( \left(\frac{1}{\tau^2}\right)^{\eta_0/2 - 1} \exp \left\{ -\frac{\eta_0 \tau_0^2}{2\tau^2} \right\} \right) \times \left( \left(\frac{1}{\tau^2}\right)^{m/2} \exp \left\{ -\frac{ \sum_{j=1}^m (\theta_j - \mu)^2 }{2\tau^2} \right\} \right) \\ &= \left( \frac{1}{\tau^2} \right)^{(\eta_0 + m)/2 - 1} \times \exp \left\{ -\frac{1}{2\tau^2} \left( \eta_0 \tau_0^2 + \sum_{j=1}^m (\theta_j - \mu)^2 \right) \right\} \\ &\propto \text{dgamma} \left( \frac{1}{\tau^2} | \frac{\eta_0 + m}{2}, \frac{\eta_0 \tau_0^2 + \sum_{j=1}^m (\theta_j - \mu)^2}{2} \right) \\ \end{align*}

Full conditional distributions of \(\theta_j\)

\(\begin{aligned} p(\theta_j | \mu, \tau^2, \sigma^2, \boldsymbol{y}_1, \dots, \boldsymbol{y}_m) \propto p(\theta_j | \mu, \tau^2) \prod_{i=0}^{n_j} p(y_{i,j} | \theta_j, \sigma^2) \end{aligned}\) より、 \[ \{\theta_j | y_{1,j}, \dots, y_{n_j,j}, \sigma^2 \} \sim \text{normal}(\frac{ n_j\bar{y}_j/\sigma^2 + \mu/\tau^2 }{ n_j/\sigma^2 + 1/\tau^2 }, [n_j/\sigma^2 + 1/\tau^2]^{-1}) \]

\(\mu\)の完全条件付き分布と同様に、 section 5.2 と全く同じ方法で導出できる。 \[\begin{aligned} y_i \sim p(\theta, \sigma^2) &\rightarrow y_{i,j} \sim p(\theta_j, \sigma^2) \\ \theta \sim p(\mu, \tau_0^2) &\rightarrow \theta_j \sim p(\mu, \tau^2) \\ \end{aligned}\] とおきかえるだけ。

Full conditional of \(\sigma^2\)

\[\begin{aligned} p(\sigma^2 | \theta_1, \dots, \theta_m, \boldsymbol{y}_1, \dots, \boldsymbol{y}_m) &\propto p(\sigma^2) \prod_{j=1}^m \prod_{i=1}^{n_j} p(y_{i,j} | \theta_j, \sigma^2) \\ &\propto (\sigma^2)^{-(\nu_0/2 + 1)} \exp \left\{ -\frac{\nu_0 \sigma_0^2}{2 \sigma^2} \right\} (\sigma^2)^{ \frac{-\sum_{j=1}^m n_j}{2} } \exp \left\{ -\frac{1}{2 \sigma^2} \sum_{j=1}^m \sum_{i=1}^{n_j} (y_{i,j} - \theta_j)^2 \right\} \\ &= (\sigma^2)^{- \left( \frac{1}{2} \left[\nu_0 \sum_{j=1}^m n_j \right] + 1 \right)} \exp \left\{ -\frac{1}{2 \sigma^2} \left[ \nu_0 \sigma_0^2 + \sum_{j=1}^m \sum_{i=1}^{n_j} (y_{i,j} - \theta_j)^2 \right] \right\} \\ \end{aligned}\]

より、

\[ \left\{ \frac{1}{\sigma^2} | \boldsymbol{\theta}, \boldsymbol{y}_1, \dots, \boldsymbol{y}_m \right\} \sim \text{gamma} \left( \frac{1}{2} \left[\nu_0 \sum_{j=1}^m n_j \right], \frac{1}{2} \left[ \nu_0 \sigma_0^2 + \sum_{j=1}^m \sum_{i=1}^{n_j} (y_{i,j} - \theta_j)^2 \right] \right) \]

完全条件付き分布のまとめ

\[\begin{aligned} \{\mu | \theta_1, \dots, \theta_m, \tau^2\} &\sim \text{normal}(\frac{ m\bar{\theta}/\tau^2 + \mu_0/\gamma_0^2 }{ m/\tau^2 + 1/\gamma_0^2 }, [m/\tau^2 + 1/\gamma_0^2]^{-1}) \\ \{1/\tau^2 | \theta_1, \dots, \theta_m, \mu\} &\sim \text{gamma}(\frac{\eta_0 + m}{2}, \frac{\eta_0 \tau_0^2 + \sum_{j=1}^m (\theta_j - \mu)^2}{2}). \\ \{\theta_j | y_{1,j}, \dots, y_{n_j,j}, \sigma^2 \} &\sim \text{normal}(\frac{ n_j\bar{y}_j/\sigma^2 + \mu/\tau^2 }{ n_j/\sigma^2 + 1/\tau^2 }, [n_j/\sigma^2 + 1/\tau^2]^{-1}) \\ \left\{ 1/\sigma^2 | \boldsymbol{\theta}, \boldsymbol{y}_1, \dots, \boldsymbol{y}_m \right\} &\sim \text{gamma} \left( \frac{1}{2} \left[\nu_0 \sum_{j=1}^m n_j \right], \frac{1}{2} \left[ \nu_0 \sigma_0^2 + \sum_{j=1}^m \sum_{i=1}^{n_j} (y_{i,j} - \theta_j)^2 \right] \right) \end{aligned}\]

参考文献

Hoff, Peter D., 入江 薫, 菅澤 翔之助, and 橋本 真太 郎. 2022. 標準ベイズ統計学. 東京: 朝倉書店.
久保川 達也. 2017. 現代数理統計学の基礎. 共立講座数学の魅力 ; 11. 東京: 共立出版.