はじめに
統計学の勉強，研究を行っている大学院生です．
これまででは順序プロビットモデルの基本的なモデルと，Albert and Chib (1993) で提案された標準的なデータ拡大の方法 (AC93) ，再パラメータ化を用いて推定する方法について紹介してきました．
今回はデータ拡大の方法とブロックサンプリングを用いた MH 法について紹介します．
PX-DA
Liu and Wu (1999) の parameter expansion for data augmentation (PX-DA) を順序プロビットモデルに適用することを考える．Schliep and Hoeting (2015) では空間構造を入れた順序プロビットモデルに対するデータ拡大を扱っている．以下では Schliep and Hoeting (2015) における通常の順序プロビットモデルを推定する方法に基づくアルゴリズムを記述する．
この方法では，誤差項の分散を $\sigma^2=1$ に固定する識別制約を緩和し，分散 $\sigma^2$ を未知パラメータとして扱う Working Model を用いてサンプリングを行う．これにより，カットオフと潜在変数の間の強い相関が緩和され，アルゴリズムの収束が加速する．
モデル
元のモデルでは $\epsilon_i \sim N(0, 1)$ であるが，拡大モデルでは以下のように誤差分散を $\sigma^2$ とする．
z_i = x_i \tilde{\beta} + \epsilon_i, \quad \epsilon_i \stackrel{iid}{\sim} N(0, \sigma^2)
ここで $\tilde{\beta}$ と $\sigma^2$ は拡大モデルにおけるパラメータである．カットオフについても拡大パラメータ $\tilde{\alpha}$ を用いるが，位置の識別のために $\tilde{\alpha}_1 = 0$ は固定し，
\tilde{\alpha} = (\tilde{\alpha}_2, \dots, \tilde{\alpha}_{M-1})'
を推定する．
事前分布として，$\tilde{\beta}$ には平坦な事前分布（or 正規分布），$\sigma^2$ には逆ガンマ分布 $IG(n_0/2, S_0/2)$ を仮定する．
サンプリングの方法
拡大モデルのパラメータ $(\tilde{\beta}, \tilde{\alpha}, \sigma^2, z)$ をサンプリングし，各イテレーションの最後に元の識別モデルのパラメータ $(\beta, \alpha)$ へ変換する．
1. 潜在変数 $z$ のサンプリング
各 $i$ について，分散 $\sigma^2$ を持つ切断正規分布からサンプリングする．
z_i | \tilde{\beta}, \tilde{\alpha}, \sigma^2, y \sim TN_{(\tilde{\alpha}_{y_i-1}, \tilde{\alpha}_{y_i}]}(x_i\tilde{\beta}, \sigma^2)
2. 係数 $\tilde{\beta}$ のサンプリング
分散 $\sigma^2$ を考慮した線形回帰モデルとしてサンプリングする．
\tilde{\beta} | z, \sigma^2, y \sim N((X'X)^{-1}X'z, \sigma^2(X'X)^{-1})
3. 拡大パラメータ $\sigma^2$ のサンプリング
線形回帰モデルの分散として逆ガンマ分布からサンプリングする．残差平方和 $SSE = (z - X\tilde{\beta})'(z - X\tilde{\beta})$ として
\sigma^2 | z, \tilde{\beta}, y \sim IG \left( \frac{n + n_0}{2}, \frac{SSE + S_0}{2} \right)
4. カットオフ $\tilde{\alpha}$ のサンプリング
AC93 と同様に，順序制約を満たす一様分布からサンプリングする．
\tilde{\alpha}_k | \tilde{\alpha}_{-k}, z, y \sim U[\underline{\tilde{\alpha}}_k, \bar{\tilde{\alpha}}_k]
ここで下限・上限は拡大モデルの潜在変数 $z$ に基づいて計算される．
\begin{aligned}
\underline{\tilde{\alpha}}_k &= \max \left( \tilde{\alpha}_{k-1}, \max_{\{i: y_i = k\}} z_i \right) \\
\bar{\tilde{\alpha}}_k &= \min \left( \tilde{\alpha}_{k+1}, \min_{\{i: y_i = k+1\}} z_i \right)
\end{aligned}
5. 元のパラメータへの変換
サンプリングされた拡大パラメータを，識別制約 $\sigma^2=1$ を満たす元のパラメータ空間へ変換して保存する．
\beta = \frac{1}{\sqrt{\sigma^2}} \tilde{\beta}, \quad \alpha_k = \frac{1}{\sqrt{\sigma^2}} \tilde{\alpha}_k
なお，次回のイテレーションは，変換前の拡大パラメータ $(\tilde{\beta}, \tilde{\alpha}, \sigma^2)$ を初期値として継続する．
この手法は，潜在変数のスケール（分散）を動的に変化させることで，カットオフパラメータが局所解に留まるのを防ぎ，AC93 に比べて高い混合性能を示す．
Cowles96 の MH 法
Albert and Chib (1993) の標準的なデータ拡大法における収束の遅さは，カットオフ $\alpha$ を更新する際に潜在変数 $z$ で条件付けていることに起因している．Cowles (1996) は，この問題を解決するために，カットオフ $\tilde{\alpha}$ と潜在変数 $z$ をブロック化して同時に更新する MH アルゴリズムを提案した．
同時事後分布の分解と提案分布
カットオフ $\tilde{\alpha}$ と潜在変数 $z$ の同時完全条件付き分布は，以下のように分解できる．
p(\tilde{\alpha}, z | \beta, y) = p(\tilde{\alpha} | \beta, y) p(z | \tilde{\alpha}, \beta, y)
ここで，第1項 $p(\tilde{\alpha} | \beta, y)$ は，$z$ を積分消去した周辺事後分布であり，これは尤度関数と事前分布の積に比例する．
p(\tilde{\alpha} | \beta, y) \propto \pi(\tilde{\alpha}) \prod_{i=1}^n \left[ \Phi(\alpha_{y_i} - x_i\beta) - \Phi(\alpha_{y_i-1} - x_i\beta) \right]
第2項 $p(z | \tilde{\alpha}, \beta, y)$ は，AC93 で用いられる切断正規分布そのものである．Cowles の手法では，この分解を利用して MH ステップの提案分布 $q(\tilde{\alpha}, z | \tilde{\alpha}^{(t)}, z^{(t)})$ を以下のように構成する．
q(\tilde{\alpha}, z | \tilde{\alpha}^{(t)}, z^{(t)}) = q_1(\tilde{\alpha} | \tilde{\alpha}^{(t)}) p(z | \tilde{\alpha}, \beta, y)
ここでは，まずカットオフの候補 $\tilde{\alpha}^{\text{prop}}$ を何らかの提案分布 $q_1$ から生成し，次にその $\tilde{\alpha}^{\text{prop}}$ に基づいて潜在変数 $z^{\text{prop}}$ をモデルの条件付き分布から生成している．この提案分布を用いると，MH アルゴリズムの採択確率における $z$ の項が相殺され，採択確率はカットオフの周辺事後分布のみに依存する形となる．
サンプリングの方法
未知のカットオフ $\tilde{\alpha} = (\alpha_2, \dots, \alpha_{M-1})'$ と潜在変数 $z$ の更新ステップは以下の通りである．
1. カットオフ $\tilde{\alpha}$ の提案と採択判定
現在の値 $\tilde{\alpha}^{(t)}$ に基づき，提案分布 $q_1(\tilde{\alpha}^{\text{prop}} | \tilde{\alpha}^{(t)})$ から候補値 $\tilde{\alpha}^{\text{prop}}$ を生成する．提案分布には，例えば切断正規分布を用いたランダムウォークが使用される．採択確率 $r$ を計算する．
r = \min \left\{ 1, \frac{p(\tilde{\alpha}^{\text{prop}} | \beta, y) q_1(\tilde{\alpha}^{(t)} | \tilde{\alpha}^{\text{prop}})}{p(\tilde{\alpha}^{(t)} | \beta, y) q_1(\tilde{\alpha}^{\text{prop}} | \tilde{\alpha}^{(t)})} \right\}
ここで $p(\tilde{\alpha} | \beta, y)$ は前述の周辺尤度を用いる．確率 $r$ で候補を受け入れ $\tilde{\alpha}^{(t+1)} = \tilde{\alpha}^{\text{prop}}$ とし，棄却した場合は $\tilde{\alpha}^{(t+1)} = \tilde{\alpha}^{(t)}$ とする．
2. 潜在変数 $z$ の更新
新しいカットオフに基づいて $z$ を再サンプリングする．
z_i^{(t+1)} \sim TN_{(\alpha_{y_i-1}^{(t+1)}, \alpha_{y_i}^{(t+1)}]}(x_i\beta, 1)
3. 係数 $\beta$ のサンプリング
AC93 と同様に，線形回帰モデルとして更新する．
\beta | z^{(t+1)}, y \sim N(\hat{\beta}, \hat{B})
この手法は，カットオフの更新時に潜在変数 $z$ の制約を受けないため，AC93 に比べて Mixing が大幅に改善する．また，$\alpha$ が棄却されたステップでは $n$ 個の潜在変数を生成する必要がないため，計算効率の面でも有利である．
おわりに
今回の記事ではデータ拡大を使う方法と Cowles の MH アルゴリズムを用いる手法を紹介しました．
次回の記事ではでは EM アルゴリズムと漸近分布を用いて Cowles の MH アルゴリズムを改善した方法について最新の研究を紹介します．
その次の記事では数値実験を行い，これらの手法の計算効率について比較します．
執筆に関して一部 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.
- Cowles, M. K. (1996). Accelerating Monte Carlo Markov chain convergence for cumulative-link generalized linear models. Statistics and Computing, 6(2), 101-111.
- Liu, J. S., & Wu, Y. N. (1999). Parameter expansion for data augmentation. Journal of the American Statistical Association, 94(448), 1264-1274.
- Schliep, E. M., & Hoeting, J. A. (2015). Data augmentation and parameter expansion for independent or spatially correlated ordinal data. Computational statistics & data analysis, 90, 1-14.