はじめに
(1) の記事では順序プロビットモデルの基本的なモデルと，Albert and Chib (1993) で提案された標準的なデータ拡大の方法 (AC93) について紹介しました．
今回は(1) に引き続き，再パラメータ化を用いる手法に関して Nandram and Chen (1996) の再パラメータ化の方法と Albert and Chib (2001) の再パラメータ化と最適化を用いる方法を紹介します．
カットオフを固定する再パラメータ化の方法 (NC96)
Nandram and Chen (1996) は，収束の遅さを改善するための再パラメータ化を提案した．通常，順序プロビットモデルでは誤差項の分散を $\sigma^2=1$ に固定して識別するが，ここでは分散を未知パラメータ $\delta^2$ として扱い，代わりにカットオフを固定している．
モデルの変換
$M=3$ のケースを考える．元のパラメータ $(\beta, \alpha_2, z)$ から，以下の変換を行う．
\delta = 1/\alpha_2, \quad \alpha^* = \delta\alpha, \quad \beta^* = \delta\beta, \quad z^* = \delta z
この変換により，変換後のカットオフ $\alpha^*$ は以下のように固定の値となる．
\alpha^*_0 = -\infty, \quad \alpha^*_1 = 0, \quad \alpha^*_2 = 1, \quad \alpha^*_3 = \infty
$r$ を $\beta$ の次元としたとき，ヤコビアンは $J=\delta^{-(n+r+2)}$ となり，変数変換後の変数
(\delta, \beta^*, z^*)
に対する同時事後分布は以下のように導出される．
p(\delta, \beta^*, z^*|y) \propto \delta^{-(n+r+2)} \exp\left( -\frac{1}{2\delta^2} \sum_{i=1}^n (z_i^* - x_i\beta^*)^2 \right) \prod_{i=1}^n I(\alpha^*_{y_i-1} < z_i^* \le \alpha^*_{y_i})
サンプリングの方法
$\alpha^*$ は固定値であることに注意して，以下のステップでギブスサンプラーを実行する．(1) と同様に定数に比例する事前分布を考える．
1. $\beta^*$ のサンプリング
分散 $\delta^2$ を持つ線形回帰モデルの係数としてサンプリングする．
\beta^* | z^*, \delta, y \sim N((X'X)^{-1}X'z^*, \delta^2 (X'X)^{-1})
2. 変換後の潜在変数 $z^*$ のサンプリング
固定された区間 $(\alpha^_{y_i-1}, \alpha^_{y_i}]$ における切断正規分布からサンプリングする．
z_i^* | \beta^*, \delta, y \sim TN_{(\alpha^*_{y_i-1}, \alpha^*_{y_i}]}(x_i\beta^*, \delta^2)
3. スケールパラメータ $\delta^2$ のサンプリング
逆ガンマ分布からサンプリングする．
\delta^2 | \beta^*, z^*, y \sim IG \left( \frac{n+r+1}{2}, \frac{1}{2} (z^* - X\beta^*)'(z^* - X\beta^*) \right)
4. 元のパラメータへの復元
各イテレーションで以下の計算を行い，元のパラメータ $(\beta, \alpha_2)$ を保存する．
\beta = \frac{1}{\delta}\beta^*, \quad \alpha_2 = \frac{1}{\delta}
この手法は，カットオフと係数の間の相関が高い場合に特に有効であり，AC93 の方法よりも高い Mixing 性能を示すことが知られている．一方，$M>3$ の場合，すべてのカットオフを固定することはできないため，一部を固定し残りの $\alpha^*$ をサンプリングする手順が必要となる．
$M > 3$ の場合への一般化
$M > 3$ の場合，識別のために最初の2つのカットオフを固定し，残りのカットオフを推定対象とする．具体的には，変換後のカットオフについて
\alpha^*_1 = 0, \alpha^*_2 = 1
と固定する．$M > 3$ の場合，残りの
\alpha^*_3, \dots, \alpha^*_{M-1}
は未知パラメータであるため，ギブスサンプラーの中にこれらをサンプリングするステップを追加する必要がある．
修正されたアルゴリズムは以下の通りである．
1-3. のパラメータのサンプリング
前述のステップ 1〜3 と同様に，
\beta^*, z^*, \delta^2
をサンプリングする．ただし，潜在変数のサンプリングにおいては，現在のカットオフの値を用いる．
4. 残りのカットオフ $\alpha^*_k$ のサンプリング
$k = 3, \dots, M-1$ について，以下の条件付き分布（一様分布）からサンプリングを行う．
\alpha^*_k | \alpha^*_{-k}, z^*, y \sim U[\underline{\alpha}^*_k, \bar{\alpha}^*_k]
ここで，$\alpha^*_{M} = \infty$ である．
\begin{aligned}
\underline{\alpha}^*_k &= \max \left( \alpha^*_{k-1}, \max_{\{i: y_i = k\}} z^*_i \right) \\
\bar{\alpha}^*_k &= \min \left( \alpha^*_{k+1}, \min_{\{i: y_i = k+1\}} z^*_i \right)
\end{aligned}
元のスケールのパラメータは以下で計算される．
\beta = \frac{1}{\delta}\beta^*, \quad \alpha_k = \frac{1}{\delta}\alpha^*_k \quad (k=2, \dots, M-1)
この拡張により，カテゴリ数が多い場合でもスケールパラメータ $\delta$ を用いたサンプリングが可能となる．ただし，追加されたカットオフのサンプリングは AC93 の方法と同じ構造を持つため，この拡張では問題点が完全には解消されない．
再パラメータ化と最適化を用いる方法 (AC01)
Albert and Chib (2001) は，順序制約のあるカットオフパラメータ $\alpha$ に対して対数変換を用いた再パラメータ化を行い，順序制約のない実数空間上のパラメータとして推定する方法を提案している．
これにより，AC93 において潜在変数 $z$ とカットオフ $\alpha$ の間で生じる依存関係を断ち切り，MH アルゴリズムを用いた効率的なサンプリングが可能となる．
再パラメータ化
順序プロビットモデルのカットオフパラメータ $\alpha$ は，$\alpha_0 = -\infty, \alpha_1 = 0, \alpha_M = \infty$ が固定されており，推定対象は $\tilde{\alpha} = (\alpha_2, \dots, \alpha_{M-1})'$ である．これらは $\alpha_1 < \alpha_2 < \dots < \alpha_{M-1}$ という順序制約を満たす必要がある．
AC01 では，この制約を解消するために以下の変換を行う．新たなパラメータ $\varphi = (\varphi_2, \dots, \varphi_{M-1})'$ を次のように定義する．
\varphi_k = \log(\alpha_k - \alpha_{k-1}), \quad k = 2, \dots, M-1
この変換により，$\varphi_k$ は実数値をとることができる．逆変換は $\alpha_1 = 0$ であることから以下で与えられる．
\alpha_k = \sum_{j=2}^k \exp(\varphi_j)
この変換に伴うヤコビアンは以下となる．
J = \left| \frac{\partial \alpha}{\partial \varphi} \right| = \prod_{k=2}^{M-1} \frac{\partial \alpha_k}{\partial \varphi_k} = \prod_{k=2}^{M-1} \exp(\varphi_k)
事後分布と提案分布の構築
AC01 のアプローチでは，カットオフ $\alpha$（変換後は $\varphi$）の更新において，潜在変数 $z$ を積分消去した周辺尤度を用いる．これにより $z$ に依存しないサンプリングが可能となる．パラメータ $\varphi$ に対する条件付き事後分布は，以下のように記述される．
\begin{aligned}
\pi(\varphi | \beta, y) &\propto L(\alpha(\varphi), \beta) \pi(\varphi) |J| \\
&= \left[ \prod_{i=1}^n \left\{ \Phi(\alpha_{y_i}(\varphi) - x_i\beta) - \Phi(\alpha_{y_i-1}(\varphi) - x_i\beta) \right\} \right] \pi(\varphi) \prod_{k=2}^{M-1} \exp(\varphi_k)
\end{aligned}
ここで $\pi(\varphi)$ は変換後のパラメータに対する事前分布であり，通常は多変量正規分布 $N(\varphi_0, V_0)$ などが用いられる．
提案分布を構築するために，Albert and Chib (2001) は事後密度関数の対数モードとその近傍での曲率を利用する．具体的には，現在の $\beta$ の下で $\log \pi(\varphi | \beta, y)$ を最大化するモード $\hat{\varphi}$ と，その点でのヘッセ行列の逆行列 $H_{\hat{\varphi}}^{-1}$ を数値的に求める．これらを用いて，MH アルゴリズムの提案分布 $q(\varphi)$ を自由度 $\nu$ の多変量 $t$ 分布などで設定する．
q(\varphi) \propto \left[ 1 + \frac{1}{\nu} (\varphi - \hat{\varphi})' H_{\hat{\varphi}} (\varphi - \hat{\varphi}) \right]^{-(\nu + M - 2)/2}
サンプリングの方法
AC01 に基づく MCMC アルゴリズムは，Gibbs サンプラーと MH アルゴリズムを組み合わせた手法となる．
1. カットオフ $\alpha$（$\varphi$）のサンプリング
現在の値 $\beta$ を所与として，$\varphi$ を更新する．提案分布 $q(\varphi | \hat{\varphi}, H_{\hat{\varphi}}^{-1})$ から候補値 $\varphi^{\text{prop}}$ を生成する．採択確率 $r$ は以下で計算される．
r = \min \left\{ 1, \frac{\pi(\varphi^{\text{prop}} | \beta, y) q(\varphi^{\text{curr}})}{\pi(\varphi^{\text{curr}} | \beta, y) q(\varphi^{\text{prop}})} \right\}
確率 $r$ で候補値を採択し $\alpha$ を更新する．不採択の場合は前回の値を維持する．
2. 潜在変数 $z$ のサンプリング
更新された $\alpha$ を用いて，AC93 と同様に切断正規分布からサンプリングする．
z_i | \beta, \alpha, y \sim TN_{(\alpha_{y_i-1}, \alpha_{y_i}]}(x_i\beta, 1)
3. 係数 $\beta$ のサンプリング
更新された $z$ を用いて，線形回帰モデルとしてサンプリングする．
この手法は，カットオフの更新に $z$ を条件付けないため，AC93 のような狭い区間でカットポイントが動かなくなる問題は発生せず，混合が大幅に改善する．一方で，ステップごとに数値最適化を行う計算コストが発生する．
おわりに
今回の記事では再パラメータ化を行う手法に関連して二つの方法を紹介しました．
NC96 の手法は $M=3$ の場合には有効ですが，$M>3$ の場合には AC93 と同様の問題が発生してしまいます．
AC01 の手法はサンプリングの効率という点では従来の手法よりも良くなりそうですが，ステップごとに数値最適化が必要となるので，計算時間は長くなってしまいます．
この数値計算が計算時間及ぼす影響については (5) で行う数値実験の記事で確認します．
以降の記事の予定です．
(3) では MH アルゴリズムを用いる方法と prameter expantion を用いる方法を紹介します．
(4) では EM アルゴリズムと漸近分布を用いた MH アルゴリズムの方法を紹介します．
(5) の記事では数値実験を行い，これらの手法の計算効率について比較します．
執筆に関して一部 Gemini を使用しました．
References
- Albert, J. H., & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association, 88(422), 669-679.
- Albert, J. H., & Chib, S. (2001). Sequential ordinal modeling with applications to survival data. Biometrics, 57(3), 829-836.
- Nandram, B., & Chen, M. H. (1996). Reparameterizing the generalized linear model to accelerate Gibbs sampler convergence. Journal of Statistical Computation and simulation, 54(1-3), 129-144.