#はじめに
東京工業大学/株式会社Nospare の栗栖です．この記事では最新の研究成果(Kurisu, Ishihara and Sugasawa(2021))について紹介します．この研究では，正規分布に基づく階層モデルにおいて，事前分布に外れ値や他の分布が紛れ込んでいる場合にも対応できる信頼区間の構成方法を新たに提案し，その理論的妥当性を与えています．また数値実験と実データへの応用を通して提案手法のパフォーマンスを確認しています．

特に，本研究では経験ベイズ (empirical Bayes) の考え方とロバスト統計の方法 ($\gamma$-divergence) を組み合わせた手法を提案しています．経験ベイズの考え方を利用することで，

例えば小地域推定の問題への応用，即ち，地域別にサンプル数が少ないときに安定した推定値を与えることが可能になります．
より一般には，データセットを様々なサブグループに分けて解析を行う場合に各グループのサンプルサイズが小さくなってしまい，各グループの推定精度が低くなるという問題がありますが，他のグループの情報も適切に取り込んで各グループでの推定精度を向上させることが可能になります．

一方，伝統的な階層モデルに基づく経験ベイズ法ではデータが正規分布を通じて生成されていることが仮定されます．このため，実際に分析するデータが外れ値を含んでいたり，正規分布以外の分布との混合分布から生成されていたりする場合に，正規性の仮定の下で経験ベイズの考え方に基づく信頼区間では正確な分析ができない(理論上の被覆確率 (nominal coverage probability) よりも実際の被覆確率 (empirical coverage probability) が (大幅に) 小さくなってしまう) という問題が生じます．

そこで，本研究では経験ベイズの考え方にロバスト統計の考え方を組み合わせ，

データが正規分布から生成されていれば従来の方法と同様に信頼区間を構成し，
一方で正規分布と異なるデータ生成過程である場合には正規分布との乖離度をデータから自動的に学習して信頼区間の位置と幅を調整してくれる

方法を提案しています．
本研究で提案した方法のRコードは以下のリンクから利用可能です.
https://github.com/sshonosuke/GDCI

#モデル(汚染無し)

まず以下の階層正規モデルを考えます．$i=1,\dots, n$に対して，

\begin{align*}
y_i|\theta_i, \sigma_i &\sim N(\theta_i, \sigma_i^2),\ \theta_i \sim N(\mu_{\ast}, \tau_{\ast}^2).  \label{model} \tag{1}
\end{align*}

また$n$個のデータは独立であるとします．

経験ベイズ信頼区間

上記のモデルにおいて，$\psi_{\ast} = (\mu_{\ast}, \tau_{\ast}^2)$ とすると，$y_i$の値で条件付けた下での $\theta_i$ の事後分布は平均 $\tilde{\theta}_i(\psi_{\ast})$，分散 $\tilde{s}_i^2 (\psi_{\ast})$の正規分布となります．即ち，

\theta_i|y_i \sim N(\tilde{\theta}_i(\psi_{\ast}), \tilde{s}_i^2 (\psi_{\ast}))

となります．ただし，

\begin{align}
\tilde{\theta}_i(\psi_{\ast}) &= y_i - {\sigma_i^2 \over \tau_{\ast}^{2} + \sigma_i^2}(y_i - \mu_{\ast}),\ \tilde{s}_i^2 (\psi_{\ast}) = {\tau_{\ast}^2 \sigma_i^2 \over \tau_{\ast}^2 + \sigma_i^2} \label{1} \tag{2}
\end{align}

です．この結果を用いると，$\theta_i$ の$100(1-\alpha)%$ (経験ベイズ)信頼区間は

\widehat{C}_{i,1-\alpha} = (\widehat{\theta}_i - z_{\alpha/2}\widehat{s}_i, \widehat{\theta}_i + z_{\alpha/2}\widehat{s}_i)

という形で与えられます．ここで，$z_{\alpha/2}$ は標準正規分布の上側$100(\alpha/2)%$分位点であり，

\begin{align*}
\widehat{\theta_i} &= \tilde{\theta}_i(\widehat{\psi}),\ \widehat{s}_i^2 = \tilde{s}_i^2 (\widehat{\psi})
\end{align*}

です．また $\widehat{\psi}$ は $y_1,\dots,y_n$ の周辺尤度から計算される $\psi_{\ast}$の最尤推定量です．

ロバスト経験ベイズ信頼区間

正規分布に基づく経験ベイズ信頼区間 $\widehat{C}_{i,1-\alpha}$ は $\theta_i$ が実際に正規分布に従う場合は正しい被覆確率( coverage probability) を与えます．一方で $\theta_i$ の分布に外れ値が含まれる場合や $\theta_i$ が正規分布でない場合は coverage probability が低くなってしまい，そのパフォーマンスが悪くなることが指摘されています (Armstrong et al.(2020))．そこで，この研究では $\theta_i$ が正規分布ではない場合，特に $\theta_i$ の分布が外れ値を含んでいたり，正規分布と他の分布の混合分布で表されていたりする場合にも適切な coverage probability を与える $\theta_i$ の信頼区間を与えることを目標としています．特に事後分布の平均と分散についての関係式である Tweedies' formula (Efron (2011)) と外れ値に対して頑健な手法である$\gamma$-divergence (Jones et al.(2001)) を組み合わせた方法を提案しています．

Tweedie's formula によれば，階層正規モデルにおける $\theta_i$ の事後分布は以下で与えられます．

\begin{align*}
\tilde{\theta}_i(\psi_{\ast}) &= y_i + \sigma_i^2 {d \over dy_i}\log m(y_i),\ \tilde{s}_i^2 (\psi_{\ast}) = \sigma_i^2 + \sigma_i^4 {d^2 \over dy_i^2}\log m(y_i).
\end{align*}

ここで $m(y_i)$ は階層正規モデルにおける $y_i$ の周辺分布 $N(\mu_{\ast}, \tau_{\ast}^2 + \sigma_i^2)$ の密度関数です．実際，$m(y_i) = \phi(y_i;\mu_{\ast},\tau_{\ast}^2 + \sigma_i^2)$ (平均 $\mu_{\ast}$，分散 $\tau_{\ast}^2 + \sigma_i^2$ の正規分布の密度関数) とすると ($\ref{1}$) が得られます．

さらに事後平均，分散を外れ値や他の分布の混入に対して頑健にするため，$\log m(y_i)$ を $\gamma$-divergence を用いてロバスト化した

\mathcal{D}_{\gamma}(y_i; \psi) = {1 \over \gamma}\phi(y_i;\mu, \tau^2 + \sigma_i^2)^{\gamma}c_{\gamma}(\tau^2)

(ただし $c_{\gamma}(\tau^2) = \{2\pi(\tau^2 + \sigma_i^2)\}^{\gamma^2/2(1+\gamma)}$) に置き換えたものを考えます．このとき，$w_{\gamma}(y_i;\psi) = \phi(y_i;\mu, \tau^2 + \sigma_i^2)^{\gamma}c_{\gamma}(\tau^2)$として，

\begin{align*}
\tilde{\theta}_i^{(\gamma)}(\psi_{\ast}) &= y_i + \sigma_i^2 {d \over dy_i}\mathcal{D}_{\gamma}(y_i; \psi_{\ast}),\\
&= y_i - w_{\gamma}(y_i;\psi_{\ast}){\sigma_i^2 \over \tau_{\ast}^{2} + \sigma_i^2}(y_i - \mu_{\ast})\\
\tilde{s}_i^{2(\gamma)} (\psi_{\ast}) &= \sigma_i^2 + \sigma_i^4 {d^2 \over dy_i^2}\mathcal{D}_{\gamma}(y_i; \psi_{\ast})\\
&= \sigma_i^2 + w_{\gamma}(y_i;\psi_{\ast}){\sigma_i^4 \over (\tau_{\ast}^2 + \sigma_i^2)^2}\left\{\gamma(y_i - \mu_{\ast})^2 - (\tau_{\ast}^2 + \sigma_i^2)\right\}.
\end{align*}

となります．これらを用いることでロバスト化した経験ベイズ信頼区間

\widehat{C}_{i,1-\alpha}^{(\gamma)} = (\widehat{\theta}_i^{(\gamma)} - z_{\alpha/2}\widehat{s}_i^{(\gamma)}, \widehat{\theta}_i^{(\gamma)} + z_{\alpha/2}\widehat{s}_i^{(\gamma)})

が得られます．ここで，$\widehat{\theta}_i^{(\gamma)} = \tilde{\theta}_i^{(\gamma)}(\widehat{\psi}^{(\gamma)})$，$\widehat{s}_i^{(\gamma)} = \tilde{s}_i^{(\gamma)}(\widehat{\psi}^{(\gamma)})$ であり，$\widehat{\psi}^{(\gamma)} = (\widehat{\mu}_{\gamma}, \widehat{\tau}_{\gamma}^{2})$ は

\hat{\psi}^{(\gamma)} = \text{argmax} _{\psi}\sum_{i=1}^{n}\mathcal{D}_{\gamma}(y_i; \psi)

で定義される $\psi_{\ast}$ の推定量です．

ロバスト化パラメータの選択

信頼区間 $\widehat{C}_{i,1-\alpha}^{(\gamma)}$ を実際にデータから計算するにあたり，ロバスト性をコントロールするパラメータ $\gamma$ を選択する必要があります．$\gamma$ の選択にあたり，次のような方法を考えます．

最終的な $\gamma$ の候補となる値のリスト $\Gamma = \{\gamma_1, \dots, \gamma_M\}$ ($0=\gamma_1 < \cdots < \gamma_M=1$) を用意する．
$\gamma$ の値として以下を満たす $\gamma_{op}$ を採用する．

\gamma_{op} = \text{argmin} _{\gamma \in \Gamma}\sum_{i=1}^{n}{\widehat{s}_i^{2(\gamma)} \over \sigma_i^2}.

ロバスト経験ベイズ信頼区間の性質(1)

階層正規モデル ($\ref{model}$) が正しい場合，適当な条件の下で

P(\gamma_{op} = 0) \to 1

となることがわかります．この性質は，モデルが正しく特定化されている場合にはロバスト経験ベイズ信頼区間 $\widehat{C}_{i,1-\alpha}^{(\gamma_{op})}$ は通常の経験ベイズ信頼区間 $\widehat{C}_{i,1-\alpha}$ と漸近的に等しくなることを意味しています．さらに

P\left(\theta_i \in \widehat{C}_{i,1-\alpha}^{(\gamma_{op})}\right) = 1-\alpha + o(1)

となり，漸近的に正しい coverage probability を与えることがわかります．

#モデル(汚染有り)

次に $\theta_i$ の分布が (正規分布)+(外れ値 or 他の分布) として表現される場合，即ち以下のような場合を考えます．

\begin{align}
y_i|\theta_i, \sigma_i &\sim N(\theta_i, \sigma_i^2),\ \theta_i \sim f_{\psi_{\ast}} \label{model2} \tag{3}
\end{align}

ここで$f_{\psi_{\ast}}(\theta_i)$は以下で与えられる確率密度関数をもつ分布に従うとします．

f_{\psi_{\ast}}(\theta_i) = (1-\omega)\phi(\theta_i; \mu_{\ast}, \tau_{\ast}^2) + \omega \delta,\ \omega \in (0,1).

特に $\delta$ が一点分布や離散分布の場合は $\theta_i$ が外れ値を含む場合に相当します．

ロバスト経験ベイズ信頼区間の性質(2)

モデル ($\ref{model2}$) の場合，適当な条件の下で

P(\gamma_{op}>0) \to 1.

が成り立ちます．これは$\theta_i$ が正規分以外の成分をもつ場合，信頼区間 $\widehat{C}_{i,1-\alpha}^{(\gamma_{op})}$ は適切にロバスト化されることを意味しています．また，$\gamma>0$に対して

P\left(\theta_i \in \widehat{C}_{i,1-\alpha}^{(\gamma)}\right) = 1-\alpha + O\left( (\omega \rho)^{1/4} + \gamma\right)

となることがわかります．ここで，$\rho$ は正規分布の密度関数 $\phi(\theta_i; \mu_{\ast}, \tau_{\ast}^2)$ とその他の密度関数 $\delta$ がどれほど"分離可能"であるかを測るパラメータです．外れ値が大きな値であったり，密度関数 $\delta$ と $\phi(\theta_i; \mu_{\ast}, \tau_{\ast}^2)$ の位置が大きくずれていたりすると $\rho$ の値は小さくなります．この場合，$\delta$ による"汚染割合" $\omega$ が小さくなくても $(\omega \rho)^{1/4}$ の項は小さくなります．

また $\rho$ がそれほど小さくなくとも，$\delta$ による汚染割合 $\omega$ が小さければ，この場合も $(\omega \rho)^{1/4}$ の項は小さくなります．

$O(\gamma)$ の項は通常の経験ベイズ信頼区間をロバスト化したことによる影響で，信頼区間 $\widehat{C}_{i,1-\alpha}^{(\gamma)}$ の coverage probability を増加させる方向に働きます．

#数値実験

以下では3つのシナリオに関して提案手法の有限標本でのパフォーマンスを確認してみます．

$y_i \sim N(\theta_i, 1)$, $\theta_i \sim N(0,V)$,
$y_i \sim N(\theta_i, 1)$, $\theta_i \sim 0.95N(0,V) + 0.05N(7,V)$,
$y_i \sim N(\theta_i, 1)$, $\theta_i \sim 0.9N(0,0.04) + 0.1N(10V, 0.04)$.

ただし，$V \in \{1,2\}$, $i=1,\dots, n$, $n=200$ とします．

下の表はそれぞれのシナリオにおける，

提案手法 (GD)，
経験ベイズ信頼区間 (EB)，
Efron(2016) (EF),
Armstrong et al.(2020) (AKM)

のそれぞれの方法での coverage probability (CP) と 95% 信頼区間の average length (AL) を与えています．それぞれのシナリオについて 1000回繰り返し計算を行っています．論文中では他の設定や先行研究の手法 (Armstrong et al.(2020), Efron (2016))との比較も行っています．さらに応用例として，東京都の都心部の犯罪件数のデータの分析を行い，先行研究との比較も含めた分析結果を論文中で紹介しています．

下の表を見ると，

事前分布に汚染がない場合 (シナリオ1) はロバストネスを調整するパラメータ $\gamma$ の値が $0$ として推定されている (経験ベイズ信頼区間と一致する) こと，
信頼区間の coverage probaiblity に関して，提案手法(GD)は経験ベイズ信頼区間(EB)と比べて同じ or 高くなっていること，
信頼区間の average length に関して，提案手法(GD)の方が経験ベイズ信頼区間(EB)よりも狭いこと

などがわかります．何れの結果も理論的結果と整合的であり，特に $\gamma$ が正の値で推定された場合の結果はロバスト経験ベイズ信頼区間の性質(2)における $O(\gamma)$ の項の影響がプラスに働く (coverage probability を向上させる方向に働く) ことを示唆しています．

シナリオ	V	$\gamma$	GD (CP)	EB (CP)	EF (CP)	AKM (CP)	GD (AL)	EB (AL)	EF (AL)	AKM (AL)
(1)	0.1	0.00	94.6	94.6	95.5	72.0	1.36	1.36	1.39	1.03
	1	0.00	94.5	94.5	93.4	97.6	2.76	2.76	2.70	3.23
(2)	0.1	0.17	99.4	94.6	98.8	96.7	2.76	3.30	1.80	3.86
	1	0.09	95.5	94.7	94.1	97.2	3.40	3.44	2.80	3.99
(3)	0.1	0.00	93.0	93.0	93.6	79.8	1.43	1.43	1.44	1.23
	1	0.14	99.4	95.0	99.7	97.1	2.56	3.72	1.73	4.18

#まとめ
この記事では Kurisu, Ishihara and Sugasawa(2021) の研究紹介として，階層正規モデルの事前分布に外れ値や他の分布が紛れ込んでいる場合にも頑健な信頼区間の構成方法と数値実験の結果について解説しました．現在関連するテーマの研究にも取り組んでいるのでまた機会があれば紹介したいと思います．

株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．

参考文献
[1] Armstrong, T.B., Kolesar, M., and Plagborg-Moller, M. (2020). Robust empirical Bayes confidence intervals. Working paper version. R&R at Econometrica.
[2] Efron, B. (2011). Tweedies's formula and section bias. Journal of the American Statistical Association 106, 1602-1614.
[3] Efron, B. (2016). Empirical Bayes deconvolution estimates. Biometrika 103, 1-20.
[4] Jones, M., Hjort, N.L., Harris, I.R., and Basu, A. (2001). A comparison of related density based minimum divergence estimators. 88, 865-873, Biometrika.
[5] Kurisu, D., Ishihara, T., and Sugasawa, S. (2021). Robust and efficient empirical Bayes confidence intervals using $\gamma$-divergence. arXiv:2108.11551.