東京大学・株式会社Nospareの菅澤です．
前回の外れ値とロバスト推定(2)では線形回帰モデルのロバスト推定について紹介しました．今回は，より一般の(パラメトリックな)統計モデルに対するロバスト推定の方法について紹介します．
線形回帰モデルにおけるHuberの方法とその一般化
尤度方程式
まず具体的な統計モデルとして線形回帰モデル
$$
y_i\sim N(x_i^\top \beta, \sigma^2), \ \ \ i=1,\ldots,n
$$
を考えます．まずは簡単のため，分散パラメータ$\sigma^2$は既知として議論します．
このモデルに対する対数尤度関数は以下で与えられます．
\ell(\beta, \sigma^2)=-\frac{n}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-x_i^\top\beta)^2
この関数を$\beta$で微分することで尤度方程式
\frac{\partial}{\partial\beta}\ell(\beta, \sigma^2)=\sum_{i=1}^n s_i(\beta,\sigma^2)=0, \ \ \ s_i(\beta,\sigma^2)=\frac{x_i}{\sigma^2}(y_i-x_i^\top\beta)
が得られます．ここで，$s_i(\beta,\sigma^2)$はスコア関数と呼ばれています．上記の方程式を$\beta$について解くことで$\beta$の最尤推定量が得られます．
Huberの方法
前回の記事で紹介したように，外れ値に対してロバストな回帰係数$\beta$の推定方法としてHuberの方法がありました．スコア関数を使うと，Huberの推定量は以下のような推定方程式の解として与えられます．
\sum_{i=1}^n w_is_i(\beta,\sigma^2)=0, \ \ \ \ \
w_i=\begin{cases}
1 & (\sigma^{-1}|y_i-x_i^\top \beta|<K) \\
\sigma K/|y_i-x_i^\top \beta| & (\sigma^{-1}|y_i-x_i^\top \beta|\geq K)
\end{cases}
これは尤度方程式に$w_i$という重みを導入することで調整した重み付き尤度方程式になっていることがわかります．重み$w_i$の形から，標準化残差$\sigma^{-1}|y_i-x_i^\top \beta|$が一定以上大きくなると$w_i$が小さくなっていくことがわかります．これは，残差が大きい外れ値に対応するスコア関数の影響を小さくしていることに相当します．極端なケースとして$|y_i-x_i^\top \beta|\to \infty$ (回帰線から無限遠離れた外れ値)を仮想的に考えると$w_i\to 0$となり，重み付き尤度方程式において外れ値が完全に無視されます．
密度関数に基づいた重み
Huberの方法は残差を用いた直感的な重みの設計になっていますが，これは線形回帰モデル特有のもので一般の統計モデルに拡張できる考えではありません．また，実際は$\sigma^2$のロバスト推定も同時に行う必要がありますが，Huberの方法は必ずしも自然な同時推定の方法を与えていません．($\sigma^2$のロバスト推定は別で考える必要があります．)
そこで，汎用的な重みとしてモデル$N(x_i^\top \beta, \sigma^2)$の密度関数に基づいた重みを考えます．パラメータ$(\beta,\sigma^2)$が所与のもと，密度関数の値が極端に小さいデータは外れ値とみなすことができるので，以下の形の重み$w_i$を用いることが考えられます．
w_i=\phi(y_i;x_i^\top \beta, \sigma^2)^{\alpha}, \ \ \ \ \alpha\geq 0
ここで，$\phi(x;a,b)$は平均$a$，分散$b$の正規分布の密度関数であり，$\alpha$はチューニングぱパラメータです．$\alpha>0$のもとで$w_i$は残差$|y_i-x_i^\top \beta|$の減少関数になり，$\alpha$の値が大きいと$w_i$はより速く減少します．すなわち，大きい$\alpha$の値を用いると外れ値に対してより小さい重みを与える(外れ値の情報をより積極的に取り除く)ことになります．
この$w_i$を用いたときの重み付き尤度関数は以下のようになります．
\sum_{i=1}^n \phi(y_i;x_i^\top \beta, \sigma^2)^{\alpha}s_i(\beta,\sigma^2)=0 \\
\Leftrightarrow
\sum_{i=1}^n \phi(y_i;x_i^\top \beta, \sigma^2)^{\alpha}\frac{\partial}{\partial\beta}\log \phi(y_i;x_i^\top \beta, \sigma^2)=0. \ \ \ \ \ \ \ (1)
$E[\phi(y_i;x_i^\top \beta, \sigma^2)s_i(\beta,\sigma^2)]=0$となるので，この方程式は不偏な推定方程式(推定関数の期待値が0)になっています．また，この推定方程式は$\alpha=0$で通常の尤度方程式に戻ります．
この推定方程式の面白い点としては，この重み付き推定方程式を導出する目的関数が定まる点です．具体的には
\frac{1}{\alpha}\sum_{i=1}^n \phi(y_i;x_i^\top \beta, \sigma^2)^{\alpha}
を$\beta$で微分することで重み付き推定方程式(1)が得られます．
次節では，このようなアイデアを一般の統計モデルに対して議論してみます．
Density Power Divergence
重み付き尤度方程式と目的関数
データ$y_1,\ldots,y_n$に対して，統計モデル$f(\cdot;\theta)$のロバストな当てはめを考えます．前節で導入した重み付き推定方程式を一般化すると以下のようになります．
\sum_{i=1}^n f(y_i;\theta)^{\alpha}\frac{\partial}{\partial\theta}\log f(y_i;\theta)=0.
$\alpha=0$として得られる通常の尤度方程式は不偏な推定方程式ですが，$\alpha>0$の場合は一般的に不偏な推定方程式になりません．(前節の例のように$\alpha>0$のときも不偏な推定方程式になるケースもあります．) そのため，バイアス調整をした以下のような推定方程式を扱う必要があります．
\sum_{i=1}^n \bigg\{f(y_i;\theta)^{\alpha}\frac{\partial}{\partial\theta}\log f(y_i;\theta) - \underbrace{E\left[f(y_i;\theta)^{\alpha}\frac{\partial}{\partial\theta}\log f(y_i;\theta)\right]}_{\text{バイアス調整項}} \bigg\}=0 \\
\Leftrightarrow
\sum_{i=1}^n \left\{f(y_i;\theta)^{\alpha}\frac{\partial}{\partial\theta}\log f(y_i;\theta) -\int f(x;\theta)^{1+\alpha}\frac{\partial}{\partial\theta}\log f(x;\theta)dx \right\}=0. \ \ \ \ \ \ \ \ (2)
$f(\cdot;\theta)$が確率関数の場合，上式で現れる積分はサポート上の和として定義されます．重み付き推定方程式(2)を解くことでロバストな$\theta$の推定量を得ることができます．
また，方程式(2)に対応する目的関数は以下のようになります．
D_{\alpha}(\theta)=\frac{1}{\alpha}\sum_{i=1}^n f(y_i;\theta)^{\alpha} - \frac{n}{1+\alpha}\int f^{1+\alpha}(x;\theta)dx. \ \ \ \ \ (3)
実際に(3)を$\theta$で微分することで推定方程式(2)を導出することができます．
この目的関数に対して
\lim_{\alpha\to 0}\left[D_{\alpha}(\theta)+n\left(1-\frac{1}{\alpha}\right) \right]
=\sum_{i=1}^n\log f(y_i;\theta)
が成り立つため，小さい$\alpha$を用いた$D_{\alpha}(\theta)$による推定は，通常の(非ロバストな)最尤法とほぼ同等なものになります．
ダイバージェンス
分布間の距離を測るツールとしてダイバージェンス(divergence)と呼ばれるものがあります．有名なものとしてKLD(Kullback-Leibler divergence)があり，統計モデルとデータ生成分布の間のKLDの経験推定量(データでダイバージェンスを近似した量)を最小化したものが最尤推定量として知られています．このように最尤推定量とKLDには密接な関係があります．
ロバストな目的関数(3)に対応するダイバージェンスはDPD (density power divergece)として知られています．これはBasu et al. (1998)によって導入されました．DPDに関する詳細な議論については論文をご参照ください．Basu et al. (1998)ではDPDに基づく(目的関数(3)に基づく)ロバスト推定量の漸近的な性質についても触れています．
線形回帰モデル (再)
再び線形回帰モデルの例を考えてみましょう．線形回帰モデル$y_i\sim N(x_i^\top\beta, \sigma^2)$は独立であるが同一分布ではありませんので，目的関数(3)を直接的に用いることはできません．同一分布でない場合，目的関数(3)は以下のように再定義することができます(Ghosh and Basu, 2013)．
D_{\alpha}(\theta)=\frac{1}{\alpha}\sum_{i=1}^n f_i(y_i;\theta)^{\alpha} - \frac{1}{1+\alpha}\sum_{i=1}^n \int f_i^{1+\alpha}(x;\theta)dx.
ここで$f_i(\cdot;\theta)$は各$i$ごとに異なる密度関数を表します．
線形回帰モデルに対して上記の目的関数は
D_{\alpha}(\beta,\sigma^2)=\frac{1}{\alpha}\sum_{i=1}^n \phi(y_i;x_i^\top\beta,\sigma^2)^{\alpha} - n(2\pi\sigma^2)^{-\alpha/2}(1+\alpha)^{-3/2}.
となります．この関数を$(\beta,\sigma^2)$について最小化することで回帰係数と誤差分散を同時にロバスト推定することができます．ちなみに$\beta$について上記の目的関数を微分すると(1)の推定方程式が得られます．
補足
他のダイバージェンス
今回は一般的なロバスト推定の方法としてDPDのみを紹介しましたが，他にも効果的なダイバージェンスがいくつか提案されています．代表的なものとしては$\gamma$-divergence (Jones et al., 2001)です．ロバスト性の観点でDPDよりも$\gamma$-divergenceの方が好ましい性質をもつことがFujisawa and Eguchi (2008)の論文で明らかにされています．
チューニングパラメータの選択
またDPDを使う際にはチューニングパラメータ$\alpha$を指定する必要があります．このパラメータの選択の際には効率性をロバスト性のトレードオフに注意する必要があります．最尤推定量は効率的な推定量として知られていますが，その効率性をある程度犠牲にしてロバスト性を担保するのがDPDによる推定になっています．
- 大きい$\alpha$: 外れ値に対するロバスト性が強くなるが，外れ値がないときの統計的な効率性が落ちる
- 小さい$\alpha$: 外れ値がないときの効率性は最尤法と同等になるが，外れ値に対するロバスト性が弱くなる
基本的にデータにどの程度外れ値が混入するかは未知なので，データに適合的に$\alpha$を選択する必要があります．これについてはいくつかの方法が提案されており，漸近リスクに基づく方法(Basak et al., 2021)や選択規準に基づく方法(Sugasawa and Yonekura, 2021)があります．
おわりに
今回はDPDを用いた一般的なロバスト推定法について紹介しました．その他のダイバージェンスを使ったロバスト推定や前回の記事で紹介したロバスト推定法との数値的な比較などは，今後の記事で紹介していこうと思います．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
