こんにちは，株式会社Nospare・明治大学の小林です．本記事ではAnceshi et al. (2023)に基づき，二項プロビットモデルを中心に，事前分布と事後分布の分布族が同じになる共役性を持つようなモデルの定式化について紹介します．
背景
二項プロビットモデルの基本
プロビットモデルは，$0$か$1$の値を取る被説明変数に対する回帰モデルとして，ロジットモデルと並んで広く使われてきました．プロビットモデルの定式化は，観測されない潜在変数$z_i,\ i=1,\dots,n$を導入し，回帰モデル
z_i = x_i'\beta + \epsilon_i, \quad \epsilon\sim N(0,1)
を考え（$x_i$は$p$次元の説明変数のベクトル，$\beta$は係数パラメータ），二値変数$y_i$は潜在変数$z_i$の値が$0$を超える場合に$1$，超えない場合に$0$を取るようにします：
y_i = \left\{
\begin{array}{cc}
1 & \text{$z_i>0$の場合}\\
0 & \text{そうでない場合}
\end{array}
\right.
$y_i$の値に対応した$z_i$の領域で$z_i$を積分すると，$\Pr(y_i=1)=\Phi(x_i'\beta)$, $\Pr(y_i=0)=1-\Phi(x_i'\beta)$となり（$\Phi(\cdot)$は標準正規分布の分布関数），$y_i=1$となる確率を$x_i$で説明するモデルであることがわかります．似たようなモデルとしてロジットモデルもあります．以前ロジットモデルに関する記事も書いてあります．
プロビットモデルの尤度関数は
p(\mathbf{y}|\beta)=\prod_{i=1}^n\Phi(x_i'\beta)^{y_i}(1-\Phi(x_i'\beta))^{1-y_i}
と書けます（ここで$\mathbf{y}=(y_1,\dots,y_n)'$）．ベイズ推定を行う場合には$\beta$に対して事前分布を仮定することになりますが，尤度関数が複雑な形状をしているため，$\beta$の事後分布を解析的に調べることはできません．一方でマルコフ連鎖モンテカルロ法（MCMC）に基づく事後推定は簡単に行うことができます．潜在変数$z_i$を導入し，例えば$\beta$に正規事前分布を仮定すると，データ拡大法により$z_i$と$\beta$のサンプリングを交互に行うことで簡単に事後推定を行うことができます．プロビットモデルに対するデータ拡大法のオリジナルの論文はChib and Greenberg (1993)です．
共役性に対する動機
二項プロビットモデルと同様の定式化は，多変量プロビットモデル，多項プロビットモデル，トービットモデルなどに対しても行うことができ，データ拡大法によって簡単に推定を行うことができます．ところが，データ拡大法では各$i$ごとに潜在変数が導入されているため，数値計算のスケーラビリティがなかったり，高次元における効率性に問題があります．もしこれらのモデルに対して共役事前分布を特定することができれば，事後分布も事前分布と同じ分布族になるので解析的に事後分布を調べることが可能となり，これらの問題を解決することができるでしょう．
二項プロビットモデルに対する共役事前分布
二項プロビットモデルの尤度関数の再定式化
二項プロビットモデルに対する共役事前分布を構築する前に，尤度関数を次のように書き換えます：
\begin{split}
p(y|\beta)&=\prod_{i=1}^n\Phi(x_i'\beta)^{y_i}(1-\Phi(x_i'\beta))^{1-y_i}\\
&=\prod_{i=1}^n\Phi\left((2y_i-1)x_i'\beta\right) \\
&=\Phi_n(\text{diag}(2\mathbf{y}-\mathbf{1}_n)\mathbf{X}\beta;\mathbf{I}_n)\quad \quad (1)
\end{split}
ここで，$\mathbf{X}=(x_1,\dots,x_n)'$, $\Phi_n(\text{diag}(2\mathbf{y}-\mathbf{1}_n)\mathbf{X}\beta;\mathbf{I}_n)$は$n$変量正規分布$N_n(\mathbf{0},\mathbf{I}_n)$の分布関数を$\text{diag}(2\mathbf{y}-\mathbf{1}_n)\mathbf{X}\beta$で評価したもの，$\mathbf{1}_n$は要素が全て$1$の$n$次元ベクトル，$\mathbf{I}_n$は$n$次元単位行列です．
この形式はより次の一般的な形式
\begin{split}
p(\mathbf{y}|\beta)&=p(\bar{\mathbf{y}}_1|\beta)p(\bar{\mathbf{y}}_0|\beta)\\
&\propto \phi_{\bar{n}_1}(\bar{\mathbf{y}}_1-\bar{\mathbf{X}}_1\beta;\bar{\mathbf{\Sigma}}_1)
\Phi_{\bar{n}_0}(\bar{\mathbf{y}}_0+\bar{\mathbf{X}}_0\beta;\bar{\mathbf{\Sigma}}_0)\quad\quad (2)
\end{split}
の特殊な場合となります．
(2)式の第一項目は$N_{\bar{n}_1}(\mathbf{0},\bar{\mathbf{\Sigma}}_1)$ の密度関数を$\bar{\mathbf{y}}_1-\bar{\mathbf{X}}_1\beta$で評価したもの，$\bar{\mathbf{y}}_1$や$\bar{\mathbf{y}}_0$は$\mathbf{y}$によって決まるものとします．
(1)式の二項プロビットモデルの場合は，(2)式において$\bar{n}_1=0$（(2)式の第一項目が消える），$\bar{n}_0=n$，$\bar{\mathbf{y}}_0=\mathbf{0}$, $\bar{\mathbf{X}}_0=\text{diag}(2\mathbf{y}-\mathbf{1}_n)\mathbf{X}$, $\bar{\mathbf{\Sigma}}_0=\mathbf{I}_n$の場合に対応します．
Unified skew-normal（SUN）事前分布
SUN分布は，尤度関数が(2)式で与えられるようなモデルに対する共役事前分布となります．SUN分布は様々な表現形式がある非対称正規（skew-normal）分布を統一的に表現したもので，$SUN_{\bar{p},\bar{n}}(\xi,\Omega,\Delta,\gamma,\Gamma)$と表記します．$\beta$がこの分布に従うとき，密度関数は
p(\beta)=\phi_{\bar{p}}(\beta-\xi;\Omega)\frac{\Phi_{\bar{n}}(\gamma+\Delta'\bar{\Omega}\omega^{-1}(\beta-\xi);\Gamma-\Delta'\bar{\Omega}^{-1}\Delta)}{\Phi_{\bar{n}}(\gamma;\Gamma)}
で与えられます．ここで，$\bar{\Omega}$は$\bar{p}\times\bar{p}$の相関行列，$\Omega$は共分散行列で，$\omega=(\Omega\odot \mathbf{I}_{\bar{p}})^{1/2}$を使って$\Omega=\omega\bar{\Omega}\omega$と書けます．$\Delta$は非対称パラメータの$\bar{p}\times\bar{n}$行列で，すべての要素がゼロの場合に，この分布は密度関数の分子と分母がキャンセルされて対称な分布となります．$\xi$は位置パラメータ，$\gamma$と$\Gamma$はそれぞれ，SUNの生成過程における打ち切りパラメータと相関行列に当たります（詳細は論文を参照してください）．
SUN事後分布
特殊な場合（論文のLemma 1）
まずは簡単な場合として尤度関数が(2)式の第一項のみ
p(\mathbf{y}|\beta)=\phi_{\bar{n}_1}(\bar{\mathbf{y}}_1-\bar{\mathbf{X}}_1\beta;\bar{\mathbf{\Sigma}}_1)\quad \quad (3)
で与えられる特殊な場合を考えます．このとき，$\beta$の事前分布を$SUN_{\bar{p},\bar{n}}(\xi,\Omega,\Delta,\gamma,\Gamma)$と仮定すると，事後分布は
\beta|\mathbf{y}\sim SUN_{\bar{p},\bar{n}}(\xi_1,\Omega_1,\Delta_1,\gamma_1,\Gamma_1)
で与えられます（パラメータの詳細は論文を参照してください）．
一般の場合（論文のTheorem 1）
尤度関数が(2)式で与えられ，$\beta$の事前分布を$SUN_{\bar{p},\bar{n}}(\xi,\Omega,\Delta,\gamma,\Gamma)$と仮定すると，事後分布はやはりSUNになります
\beta|\mathbf{y}\sim SUN_{\bar{p},\bar{n}+\bar{n}_0}(\xi_\text{post},\Omega_\text{post},\Delta_\text{post},\gamma_\text{post},\Gamma_\text{post})
事後分布がSUNで与えられるので，事後期待値や事後分散はSUNの積率母関数から得ることができます：
\begin{split}
E[\beta|\mathbf{y}]&=\xi_\text{post}+\omega_\text{post}\Delta_\text{post}\psi\\
V(\beta|\mathbf{y})&=\Omega_\text{post}+\omega_\text{post}\Delta_\text{post}(\Psi-\psi\psi')\Delta_\text{post}'\omega_\text{post}
\end{split}
ここで，$\psi$は$\gamma_\text{post}$, $\Gamma_\text{post}$に依存する$(\bar{n}+\bar{n}_0)\times 1$次元ベクトル，$\Psi$はSUNの積率母関数の2階微分を含む$(\bar{n}+\bar{n}_0)\times (\bar{n}+\bar{n}_0)$行列です．
これらの数量は次の表現に基づいてモンテカルロ法によって推定することも可能です：
\begin{split}
\beta|\mathbf{y}& \buildrel d \over = \xi_\text{post}+ \omega_\text{post}(\mathbf{U}_0+\Delta_\text{post}\Gamma_\text{post}^{-1}\mathbf{U}_1)\\
\mathbf{U}_0&\sim N_{\bar{p}_0} (\mathbf{0},\bar{\Omega}_\text{post}-\Delta_\text{post}\Gamma_\text{post}^{-1}\Delta_\text{post}')\\
\mathbf{U}_1&\sim TN_{\bar{n}+\bar{n}_0}(-\gamma_\text{post};\mathbf{0},\Gamma_\text{post})
\end{split}
ここで最後の式の分布は$-\gamma_\text{post}$より下が打ち切られた$\bar{n}+\bar{n}_0$次元の切断正規分布を表します．
モデル選択などに用いられるデータの周辺分布（周辺尤度）や，事後予測分布も得られます．
二項プロビットモデルに対するSUN事前分布
最後に，二項プロビットモデルは$\bar{p}=p$, $\bar{n}_1=0$, $\bar{n}_0=n$，$\bar{\mathbf{y}}_0=\mathbf{0}$, $\bar{\mathbf{X}}_0=\text{diag}(2\mathbf{y}-\mathbf{1}_n)\mathbf{X}$, $\bar{\mathbf{\Sigma}}_0=\mathbf{I}_n$の場合です．$\bar{n}_1=0$より(2)式の第一項がなくなるので，$\bar{y}_1$, $\bar{X}_1$, $\bar{\Sigma}_1$は定義されず，事後分布の導出の際には無視されます（$\xi_1$や$\Omega_1$の定義で出現します）．
最後に
本記事ではSUN事前分布がプロビットモデルの係数パラメータに対して共役事前分布となるということを紹介しました．本記事では主に2項プロビットに焦点を当てましたが，多変量プロビットモデル，多項プロビットモデル，トービットモデルの尤度関数も(2)式の形式に書き換えることができるため，これらのモデルに対しても$\beta$の事後分布がSUNで与えられることになります．次の記事では，SUN事前分布を使った場合での計算手法に関して紹介を行いたいと思います．
一緒にお仕事しませんか！
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．データサイエンス研修・アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
