ポリア・ガンマ分布を用いた多項ロジスティック回帰のギブスサンプラー

ポリア・ガンマ分布

まず Polson et al. (2013, [1205.0310] Bayesian inference for logistic models using Polya-Gamma latent variables) の主な結果を紹介する.

ポリア・ガンマ分布の密度関数を $p(\omega) = \mathrm{PG}(\omega | b, 0)$ とすると,

$\displaystyle \frac{(e^{\eta})^a}{(1+e^{\eta})^b} = 2^{-b} \exp( (a-b/2) \eta) \int \exp(-\omega \eta^2/2) \, p(\omega) \, d\omega$

が成り立つ. ポリア・ガンマ分布というのは, あまり聞いたことがない（と思う）けど, この性質が成り立つように定義された分布（そういうものがある）だと思えばいいだろう.

これより $\eta$ を所与にしたときの $\omega$ の分布は,

$\displaystyle p(\omega | \eta)=\frac{\exp(-(\eta^2/2) \, \omega) \, p(\omega)}{\int \exp(-(\eta^2/2) \, \omega) \, p(\omega) \, d\omega}\ \displaystyle =\mathrm{PG}(\omega | b, \eta).$

Polson et al. (2013) ではこれを使って多項ロジスティック回帰について, ギブスサンプリングのためのfull conditional distribution（自分自身以外の変数を全部所与としたときの事後分布）を導出した.

多項ロジスティック回帰

多項ロジスティック回帰, つまり次のモデルを考える：

$\displaystyle y_i \sim \mathrm{Multi}(n_i, f(x_i' \beta))$

ここで, $f(x)$ はベクトル値関数

$\displaystyle f(x) = \left(\frac{\exp(x_1)}{\sum_j \exp(x_j)}, \ldots, \frac{\exp(x_k)}{\sum_j \exp(x_j)}\right)$

である（ソフトマックス関数と呼ばれる）.

モデルの識別可能性のため, 係数 $\beta$ のどこか1列 ( $\beta_j$ ) を0に固定するのが定石である. ここでは $\beta_1 = (0,\ldots,0)'$ とする.

$\displaystyle c_{ij} = \log \sum_{j'\neq j} \exp(x_i' \beta_{j'})$

とし, $\eta_{ij} = x_i' \beta_j - c_{ij}$ とすると, 多項ロジスティック回帰の $\beta_j$ 以外の変数を所与としたときのサンプル一つあたりの尤度 $L_i$ は,

$\displaystyle L_i \propto \frac{\exp(\eta_{ij})^{y_{ij}}}{\{1+\exp(\eta_{ij})\}^{m_i}}$

先のポリア・ガンマ分布の性質を使うと,

$\displaystyle L_i \propto 2^{-m_i} \exp({y_{ij}-m_i/2} \, \eta_{ij})\int \exp(-\omega \eta_{ij}^2/2) \, p(\omega) \, d\omega$

と書け, 指数関数の中身は

$\displaystyle (y_i-n_i/2)\eta_{ij} - \omega_i\eta_{ij}^2 = -\frac{\omega_i}{2} (\eta_{ij}^2 - 2\eta_{ij}(y_{ij}-m_i/2)/\omega_i )$

$\displaystyle =-\frac{\omega_i}{2} (\eta_{ij} - (y_{ij}-m_i/2)/\omega_i )^2+(y_{ij}-m_i/2)^2/\omega_i$

と変形できる. ベクトルと行列を使ってまとめたいので,

$\eta_j=(\eta_{1j}, \ldots, \eta_{nj})'$ , $c_j=(c_{1j},\ldots, c_{nj})'$ , $z_j=(z_{1j}, \ldots , z_{nj})'$ , $z_{ij} = (y_{ij}-m_i/2)/\omega_{ij}$ と置くと,