ギブス・サンプリングの実装練習：(1-3)分散成分の事前分布を半コーシー分布にする

はじめに

これまで2回にわたってランダム効果モデルのギブス・サンプリングを説明してきました。またランダム効果 $\theta_i$ の分散 $\sigma^2_{\theta}$ について、事前分布として逆ガンマ分布や一様分布を設定した場合のベイズ推定結果をシミュレーションにより確認しました。

ところが、分散パラメータの推定に関して、推定に使用できるデータ数が少ない場合には逆ガンマ分布や一様分布ではうまくいかない可能性があることが知られています。
そこで、今回は最終回として、Gelman(2006)（参考文献[1]）で紹介されている「半コーシー分布(Half Cauchy distribution)」を分散の平方根 $\sigma_{\theta}$ （以下では「標準偏差」と言うことにします）の事前分布にすることを考えてみます*1。
この分布はコーシー分布の範囲を正の値に限定したもので、右裾の厚い形状の分布形を持ちます。
特に標準偏差に関してある程度の事前情報があるような場合（例えば、ばらつきの大きさはどれだけ大きくても10程度だろう）には、逆ガンマ分布や一様分布よりも推奨されています。

今回は、まず半コーシー分布について簡単に紹介した後、前回と同じようにギブス・サンプリングに必要な完全条件付き分布を確認します。
またRでのサンプリングコード例とともにシミュレーション結果を報告します。

なお、半コーシー分布を事前分布に設定するのは、例えばStanのようなソフトウェアであれば比較的簡単に行えるのですが、ギブス・サンプリングを自分で実装するにはちょっとした工夫が必要になります。
Gelman論文を参考に、そのあたりの解説も行っていきます。

半コーシー分布とは

半コーシー分布は、以下の確率密度関数をもつ分布です（参考文献[3]）。

$\displaystyle p(y | \mu, \sigma) = \frac{2}{\pi \sigma} \frac{1}{1 + \frac{(y - \mu)^2}{\sigma^2}} \hspace{20pt} (y \geq \mu)$

ただし $y < \mu$ における確率密度は $0$ です。
$\mu$ は分布の位置を決めるパラメータで、標準偏差の事前分布とする場合は通常 $\mu = 0$ のものが使われます。本記事でも以降は $\mu = 0$ としたものを考えることにします。
また $\sigma$ は分布の山の部分の広がりを決めるパラメータです（以下、スケールパラメータと呼ぶことにします）。このパラメータの値は定数として与えることも、さらに事前分布を考えることもあります。定数とする場合、例えば想定される標準偏差の大きさよりも少し大きめの値を設定します。

f:id:mstour:20211217204258j:plain — 半コーシー分布の確率密度関数

今回は、ランダム効果 $\theta_i$ が正規分布 $N(0, \sigma^2_{\theta})$ に従うとするモデルにて、標準偏差 $\sigma_{\theta}$ の事前分布を半コーシー分布とすることを考えます。

ところで、通常のコーシー分布については、次のような性質があります（参考文献[4]）。

$X$ を標準正規分布に従う確率変数、 $Z$ をパラメータ $(\frac{1}{2}, \frac{s^2}{2})$ の逆ガンマ分布に従う確率変数とする。この時、 $X \sqrt{Z}$ はスケールパラメータ $s$ のコーシー分布に従う*2。
スケールパラメータ $s$ のコーシー分布に従う変数に定数 $c$ を掛けたものは、スケールパラメータ $s|c|$ のコーシー分布に従う（ $|c|$ は $c$ の絶対値）。

この2つの性質から、「正規分布 $N(0, \sigma^2)$ に従う変数 $U$ 」と、「パラメータ $(\frac{1}{2}, \frac{1}{2})$ の逆ガンマ分布に従う変数 $V^2$ の平方根 $V$ 」とを掛けたものは以下の変形によりスケールパラメータ $\sigma$ のコーシー分布に従うことがわかります*3。

$\displaystyle U \sqrt{V^2} = \sigma \times \frac{U}{\sigma} V$

$U$ のとる値を正の値に制限した場合には、 $U \sqrt{V^2}$ の分布はスケールパラメータ $\sigma$ の半コーシー分布になります。
少しややこしい話になりましたが、この性質を利用することで、この後のギブス・サンプリングをうまく進めることができるようになります。

完全条件付き分布

半コーシー分布を設定した場合のギブス・サンプリングですが、サンプリングに必要となる完全条件付き分布を得るためには、元々のランダム効果モデルの表現を少し修正する必要があります。と言っても構造を変更して全く別のモデルにするわけではなく、用いられているパラメータの変換によって同じモデルを別の表現に直します。
ここで行う変換は「Parameter expansion」と呼ばれているものの一種で、サンプリングの高速化と真の分布への収束を良くするために古くから提案されているものです。

元のモデルは、観測値が以下の式で表されるとするものでした。

$\displaystyle Y_{ij} = \theta_i + e_{ij} \hspace{10pt} (i = 1, \cdots, K; \hspace{5pt} j = 1, \cdots, J)$

ただし $\theta_i \sim N(\mu, \sigma^2_{\theta}), e_{ij} \sim N(0, \sigma^2_e)$ です。

ここで、ランダム効果 $\theta_i$ を以下のように変換します。

$\displaystyle \theta_i = \mu + \xi \eta_i, \\ \eta_i \sim N(0, \sigma^2_{\eta}) \\ \hspace{10pt} (i = 1, \cdots, K)$

この時 $\theta_i$ の（ $\xi, \mu$ を与えた条件付きでの）分散は

$\displaystyle V(\theta_i) = \xi^2 V(\eta_i) = \xi^2 \sigma^2_{\eta}$

となるので、 $\sigma_{\theta} = |\xi| \sigma_{\eta}$ という対応関係が成り立ちます。

さて、 $\xi$ の事前分布を正規分布 $N(0, \sigma^2_{\xi})$ 、 $\sigma_{\eta}$ の事前分布をパラメータ $(\frac{1}{2}, \frac{1}{2})$ の逆ガンマ分布とすると、半コーシー分布の紹介のところで述べたように $|\xi|$ と $\sigma_{\eta}$ を掛けた $|\xi| \sigma_{\eta}$ はスケールパラメータ $\sigma_{\xi}$ の半コーシー分布に従います。
さらに $\xi, \sigma_{\eta}$ はそれぞれ他の全てのパラメータと観測値に対して「条件付き共役(conditionally conjugate)」、つまり完全条件付き分布も同じく $\xi$ は正規分布、 $\sigma_{\eta}$ は逆ガンマ分布になることが知られています。

したがって、正規分布と逆ガンマ分布だけを用いたギブス・サンプリングにより $\xi, \sigma_{\eta}$ のサンプリングを（他のパラメータとともに）行い、各反復で得られるサンプルを用いて $\sigma_{\theta} = |\xi| \sigma_{\eta}$ を計算することで、事前分布を半コーシー分布とする $\sigma_{\theta}$ の（事後分布からの）サンプルが得られるというわけです*4。

サンプリングに必要な完全条件付き分布の導出は、前の2記事と同様に地道に数式を展開していけばできますので、結果だけ紹介することにします。
ただし、各パラメータの事前分布は以下のように定義します。

$\displaystyle \xi \sim N(0, \sigma^2_{\xi}),$

$\displaystyle \eta_i \sim N(\mu, \sigma^2_{\eta}) \hspace{10pt} (i = 1, \cdots, K) \hspace{20pt} \eta = (\eta_1, \cdots, \eta_K)^T,$

$\displaystyle \sigma^2_{\eta} \sim IG(a_{\eta}, b_{\eta}),$

$\displaystyle \sigma^2_e \sim IG(a_e, b_e),$

$\displaystyle \mu \sim N(\mu_0, \sigma^2_0)$

なお $a_{\eta}, b_{\eta}, a_e, b_e, \mu_0, \sigma^2_0, \sigma^2_{\xi}$ は定数とします。 $\sigma_{\xi}$ が半コーシー分布の広がりを決めるスケールパラメータになります。

$\xi$ の完全条件付き分布

以下の平均・分散パラメータをもつ正規分布になります。
平均： $\frac{\sigma^2_{\xi} \sum_{i=1}^{K} \sum_{j=1}^{J} (Y_{ij} - \mu) \eta_i}{J \sigma^2_{\xi} \sum_{i=1}^{K} \eta^2_i + \sigma^2_e}$
分散： $\frac{\sigma^2_e \sigma^2_{\xi}}{J \sigma^2_{\xi} \sum_{i=1}^{K} \eta^2_i + \sigma^2_e}$

$\sigma^2_{\eta}$ の完全条件付き分布

パラメータ $a_{\eta}+\frac{K}{2}, b_{\eta}+\frac{\sum_{i=1}^K \eta_i^2}{2}$ の逆ガンマ分布になります。

$\sigma^2_e$ の完全条件付き分布

パラメータ $a_e+\frac{KJ}{2}, b_e+\frac{\sum_{i=1}^K \sum_{j=1}^J (Y_{ij} - \mu - \xi \eta_i)^2}{2}$ の逆ガンマ分布になります。

$\mu$ の完全条件付き分布

以下の平均・分散パラメータをもつ正規分布になります。
平均： $\frac{\sigma^2_0 \sum_{i=1}^K \sum_{j=1}^J (Y_{ij} - \xi \eta_i) + \sigma^2_e \mu_0}{K J \sigma^2_0 + \sigma^2_e}$
分散： $\frac{\sigma^2_0 \sigma^2_e}{K J \sigma^2_0 + \sigma^2_e}$

$\eta$ の完全条件付き分布

以下の平均ベクトル・共分散行列をもつ多変量正規分布になります。
平均ベクトル： $\left( \frac{\sigma^2_e \sigma^2_{\eta}}{J \sigma^2_{\eta} \xi^2 + \sigma^2_e} \right) \frac{J \xi}{\sigma^2_e} (\overline{Y} - \mu \boldsymbol{1})$
共分散行列： $\left( \frac{\sigma^2_e \sigma^2_{\eta}}{J \sigma^2_{\eta} \xi^2 + \sigma^2_e} \right) I$
ただし $\overline{Y} = (\overline{Y}_1, \cdots, \overline{Y}_K)^T$ 、 $\boldsymbol{1}$ は1を並べたベクトルです。

ギブス・サンプリングではこれらの分布から順番にサンプリングを行いますが、各反復の最後に現時点のサンプルを用いて以下の変換式で $\sigma_{\theta}$ の事後サンプルを得ます。

$\displaystyle \sigma_{\theta} = |\xi| \sigma_{\eta}$

では、Rで実装してみましょう。

Rコードとシミュレーション結果

今回の設定でのギブス・サンプリングは、例えば以下のようなコードで実行できます。
基本的な内容は逆ガンマ分布や一様分布の場合と同様ですが、サンプリングステップの最後（forループの最後のあたり）で $\sigma^2_{\theta}$ と $\theta$ を他のパラメータサンプルを用いて計算しています（ただし、今回は $\theta$ の結果は確認していません）。

library(tidyverse)
library(MASS)
library(MCMCpack)

set.seed(819)

#----- sample data -----
K = 6 # number of subjects
J = 5 # obs per subject
mu = 5
s2_theta = 4
s2_e = 16

theta = rnorm(K, mean = mu, sd = sqrt(s2_theta))
e = rnorm(K*J, mean = 0, sd = sqrt(s2_e))
Y = rep(theta, each = J) + e
Ymat = matrix(Y, nrow = K, ncol = J, byrow = T)
Ybar = apply(Ymat, 1, mean) # Yの各個人の平均

#----- Gibbs sampling(Half-Cauchy) -----
# MCMC設定
nsamp = 11000
burnin = 1000
thin = seq(burnin+1, nsamp, by = 1)

# hyperparameter
# mu
mu_0 = 0
s2_0 = 10000
# xi
s2_xi = 20^2 # sqrtがHalf-Cauchyのscaleになる
# sigma_e
a_e = 0.001
b_e = 0.001
# sigma_eta
a_et = 1/2
b_et = 1/2

# 結果の入れ物
xi_s = array(NA, length(thin))
s2_eta_s = array(NA, length(thin))
s2_e_s = array(NA, length(thin))
mu_s = array(NA, length(thin))
eta_s = array(NA, c(K, length(thin)))
s2_theta_s = array(NA, length(thin)) # sampleから導出
theta_s = array(NA, c(K, length(thin))) # sampleから導出

# Gibbs sampling初期値
xi_g = runif(1, min = -2, max = 2)
s2_eta_g = runif(1, min = 0, max = 2)
s2_e_g = runif(1, min = 0, max = 2)
mu_g = runif(1, min = -2, max = 2)
eta_g = runif(K, min = -2, max = 2)

### sampling
for(s in 1:nsamp){
  if((round(s/500) - (s/500)) == 0){print(s)}
  # xi
  a = (s2_xi*sum((Y-rep(mu_g, K*J))*rep(eta_g, each = J)))/(J*s2_xi*sum(eta_g^2)+s2_e_g)
  b = (s2_e_g*s2_xi)/(J*s2_xi*sum(eta_g^2)+s2_e_g)
  xi_g = rnorm(1, mean = a, sd = sqrt(b))
  xi_s[s==thin] = xi_g

  # s2_eta
  a = a_et + K/2
  b = b_et + sum(eta_g^2)/2
  s2_eta_g = rinvgamma(1, a, b) # MCMCpack
  s2_eta_s[s==thin] = s2_eta_g
  
  # s2_e
  a = a_e + (K*J)/2
  b = b_e + sum((Y-rep(mu_g, K*J)-xi_g*rep(eta_g, each = J))^2)/2
  s2_e_g = rinvgamma(1, a, b)
  s2_e_s[s==thin] = s2_e_g
  
  # mu
  a = (s2_0*sum(Y-xi_g*rep(eta_g, each = J))+s2_e_g*mu_0)/(K*J*s2_0+s2_e_g)
  b = (s2_0*s2_e_g)/(K*J*s2_0+s2_e_g)
  mu_g = rnorm(1, mean = a, sd = sqrt(b))
  mu_s[s==thin] = mu_g
  
  # eta
  coef1 = (s2_e_g*s2_eta_g)/(J*s2_eta_g*(xi_g^2)+s2_e_g)
  coef2 = (J*xi_g)/s2_e_g
  a = coef1 * coef2 * (Ybar - rep(mu_g, K))
  b = diag(rep(coef1, K))
  eta_g = mvrnorm(1, mu = a, Sigma = b)
  eta_s[, s==thin] = eta_g
  
  # s2_theta
  s2_theta_g = (abs(xi_g) * sqrt(s2_eta_g))^2
  s2_theta_s[s==thin] = s2_theta_g
  
  # theta
  theta_g = mu_g + xi_g * eta_g
  theta_s[s==thin] = theta_g
}

さて、今回は

半コーシー分布
逆ガンマ分布
一様分布

のそれぞれを事前分布にした場合の $\sigma^2_{\theta}$ のサンプリング結果を比較してみます。
なお、半コーシー分布は上記コード例の通りスケールパラメータを20（真の値である4の5倍）としました。また逆ガンマ分布と一様分布は無情報分布として用いることを想定し、逆ガンマ分布のパラメータは $(0.001, 0.001)$ 、一様分布の範囲は $(0, 100)$ に設定しました*5。
半コーシー分布の場合のプロットの作図コード例は以下の通りです。今回も稀にとても大きな値がサンプリングされるので、表示範囲を0〜20に制限しています。

val = s2_theta_s
s2_theta_df = as.data.frame(val) %>%
  mutate(Iter = 1:length(thin))
ggplot(s2_theta_df, aes(x = Iter, y = val)) + geom_line() +
  geom_hline(yintercept = s2_theta, color = "red") +
  ylim(c(0, 20)) +
  labs(y = "sigma^2_theta", title = "Half-Cauchy prior") + theme_bw()

トレースプロットを(1)半コーシー分布、(2)逆ガンマ分布、(3)一様分布、の順に示します。

f:id:mstour:20211217194127p:plain — (1)半コーシー分布の場合のトレースプロット

f:id:mstour:20211217194150p:plain — (2)逆ガンマ分布の場合のトレースプロット

f:id:mstour:20211217200145p:plain — (3)一様分布の場合のトレースプロット

他の場合と比べて、半コーシー分布を事前分布とした場合は真の値（赤色の直線で表示）に近いところにサンプルが集まっているように見えます。

サンプルから推定される事後分布の要約統計量は次のようになりました。

事前分布	平均値	中央値	95%信用区間
半コーシー分布	4.27	1.45	0.00 - 25.57
逆ガンマ分布	0.94	0.10	0.00 - 6.93
一様分布	3.74	1.26	0.00 - 23.25

今回のシミュレーション設定においては、平均値と中央値は半コーシー分布と一様分布とで近い結果になりました。一方、逆ガンマ分布に関しては分布の中心が非常に小さな値になっていることもわかります。

まとめ

ランダム効果モデルにおける分散パラメータに関して、Gelman(2006)にて事前分布として提案された「半コーシー分布」を用いる場合のギブス・サンプリングについて説明し、Rでの実装例とシミュレーション実験結果を紹介しました。
今回行った実験では、半コーシー分布と一様分布での推定結果は少なくとも要約統計量に関して言えばあまり違いはなく、半コーシー分布を用いることによる劇的な効果を示すことはできませんでした。しかしながら、半コーシー分布を用いるメリットはいろいろな所で言われており、選択肢の一つとして理解しておくと良いように思いました。

参考文献

[1] Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper). Bayesian analysis, 1(3), 515-534.
[2] 松浦健太郎(2016). StanとRでベイズ統計モデリング. 共立出版.
[3] Half-Cauchy distribution — Probability Distribution Explorer documentation
[4] Cauchy distribution - Wikipedia

*1:Gelman論文ではより一般の「半t分布」について取り上げられていますが、本記事では話を簡単にするために半コーシー分布のみを扱います。なお、t分布の自由度を1にしたものがコーシー分布です。

*2:前述の通り、位置を表すパラメータは0としています。

*3:Gelman論文の本文では、正規分布の分散を1にして、逆ガンマ分布のパラメータによって半コーシー分布のスケールパラメータを調節するような説明がされていますが、個人的には論文付録のコード例に記載の方法（逆ガンマ分布のパラメータを $(\frac{1}{2}, \frac{1}{2})$ に固定し、正規分布の標準偏差でスケールパラメータを定義する）がわかりやすかったので、以降の導出もそちらで進めています。

*4:厳密な説明ではないかもしれませんが、先はまだ長いのでこの程度でご容赦ください。

*5:サンプリングの初期値設定なども半コーシー分布の場合に揃えています。したがって、前回の記事とは結果が異なっています。