ブートストラップを用いたバイアス補正（1）

はじめに

千葉大学・株式会社Nospareの川久保です．今回と次回で，ブートストラップを用いたバイアス補正の手法について紹介します．今回は手法の一般論について説明して，次回で具体的な応用例を紹介したいと思います．

プラグイン推定量のバイアス

説明を容易にするため，パラメトリックブートストラップと呼ばれる手法に沿って説明します．本来ブートストラップ法はノンパラメトリックな手法ですが，パラメトリックブートストラップは，ブートストラップサンプルをパラメトリックな確率分布から乱数生成します．多変量解析などでは，パラメトリックモデルにおいても解析的な期待値計算が容易でないことがしばしばあり，そのような場合に用いられます．ブートストラップ法に関しては，栗栖先生の記事で最新の研究成果の紹介がなされています．

データ$y = (y_1,\dots,y_n)^\top$が，パラメータ$\theta$の確率分布$F_\theta$から生成されていると仮定します．ここで，$g(\theta)$という$\theta$の関数を推定する問題に関心があるとします．パラメータ$\theta$の推定量として$\hat{\theta}$が得られたときに，$g(\theta)$のナイーブな推定量としては，プラグイン推定量$g(\hat{\theta})$が考えられます．ところが，一般に（$\hat{\theta}$が$\theta$の不偏推定量であったとしても）$g(\hat{\theta})$は$g(\theta)$の不偏推定量になるとは限りません．

$g(\theta) = O(1)$であり，プラグイン推定量$g(\hat{\theta})$のバイアスが，
$$
b_1(\theta) = E_\theta[g(\hat{\theta})] - g(\theta) = O(n^{-1}), \tag{1}
$$
であるとします．ただし，$E_\theta[\cdot]$は，確率分布$F_\theta$についての期待値を表すとします．ここで，バイアス$b_1(\theta)$を推定し，プラグイン推定量$g(\hat{\theta})$から引いてあげることで，$g(\theta)$のバイアス補正された推定量を求めます．バイアス$b_1(\theta)$の推定量として，そのプラグイン推定量$b_1(\hat{\theta})$を考えたときにもやはりバイアスが生じてしまいますが，そのバイアスが,
$$
b_2(\theta) = E_\theta[b_1(\hat{\theta})] - b_1(\theta) = o(n^{-1}), \tag{2}
$$
というオーダーである（$n^{-1}$より速くバイアスが0に収束する）とします．このとき，バイアス補正された$g(\theta)$の推定量$g(\hat{\theta}) - b_1(\hat{\theta})$のバイアスは，

\begin{align}
E_\theta[g(\hat{\theta}) - b_1(\hat{\theta})] - g(\theta) &= E_\theta[g(\hat{\theta})] - g(\theta) - E_\theta[b_1(\hat{\theta})] \\
&= b_1(\theta) - E_\theta[b_1(\hat{\theta})] \quad (\because \ (1)式) \\
&= o(n^{-1}) \quad (\because \ (2)式)
\end{align}

となり，バイアスのオーダーが小さくなります．このとき，$g(\hat{\theta}) - b_1(\hat{\theta})$は$g(\theta)$の二次漸近不偏推定量であると言います．

バイアスのブートストラップ推定量

バイアス$b_1(\theta)$が，$\theta$の関数として陽に表現されている場合には，バイアス補正推定量$g(\hat{\theta}) - b_1(\hat{\theta})$を求めることができますが，それが難しいケースがしばしば存在します．そこで，バイアス$b_1(\theta)$をブートストラップ法で推定することを考えます．$b_1(\hat{\theta})$は，形式的には(1)式の$\theta$を$\hat{\theta}$に置き換えればいいわけですが，$E_\theta[\cdot]$を$E_\hat{\theta}[\cdot]$に置き換えるとは，どういうことを意味するでしょうか．それは，確率分布$F_\hat{\theta}$について期待値計算を行うことを意味します．そして，期待値演算の中身のプラグイン推定量$g(\hat{\theta})$は，$F_\hat{\theta}$から生成されたデータ（これをブートストラップサンプルと言います）にもとづいた推定量に置き換わります．つまり，$y^* = (y_1^\ast,\dots,y_n^\ast)^\top \sim F_\hat{\theta}$としたとき，$y^\ast = (y_1^\ast,\dots,y_n^\ast)^\top$から構成する$\theta$の推定量を$\hat{\theta}^\ast$と表すと，$g(\hat{\theta}^\ast)$に置き換わるわけです．つまり，$b_1(\hat{\theta})$は，
$$
b_1(\hat{\theta})=　E_\hat{\theta}[g(\hat{\theta}^\ast)] - g(\hat{\theta}), \tag{3}
$$
と表現できます．これは，$B$個のブートストラップサンプルを$y^\ast(1),\dots,y^\ast(B) \overset{\mathrm{iid}}{\sim} F_\hat{\theta}$とし，$y^\ast(b)$にもとづいた$\theta$の推定量を$\hat{\theta}^\ast(b)$とすると，
$$
b_1(\hat{\theta}) \approx {1 \over B}\sum_{b=1}^B g(\hat{\theta}^*(b)) - g(\hat{\theta})
$$
とモンテカルロ近似できます．このバイアス推定量を用いて，$g(\hat{\theta})$のバイアス補正を行うと，
$$
g(\hat{\theta}) - \left\{ E_\hat{\theta}[g(\hat{\theta}^\ast)] - g(\hat{\theta}) \right\} = 2g(\hat{\theta}) - E_\hat{\theta}[g(\hat{\theta}^\ast)], \tag{4}
$$
という$g(\theta)$の二次漸近不偏推定量を導出できます．

なお，似たような手法として，ジャックナイフと呼ばれる手法によってもバイアスの推定を行うことができます．

その他のバイアス補正推定量

(4)式以外の$g(\theta)$のバイアス補正推定量としては，
$$
{ \{ g(\hat{\theta}) \}^2 \over E_\hat{\theta}[ g(\hat{\theta}^\ast) ] }, \tag{5}
$$
というものも考えられます．この推定量も二次漸近不偏ですが，それは以下のようにして示すことができます．

\begin{align}
{ \{ g(\hat{\theta}) \}^2 \over E_\hat{\theta}[ g(\hat{\theta}^\ast) ] } &= { \{ g(\hat{\theta}) \}^2 \over g(\hat{\theta}) + b_1(\hat{\theta}) } \quad (\because (3)式)\\
&= {g(\hat{\theta}) \over 1 + b_1(\hat{\theta}) / g(\hat{\theta})} \\
&= g(\hat{\theta}) \left\{ 1 - {b_1(\hat{\theta}) \over g(\hat{\theta}) } + o_p(n^{-1}) \right\} \quad (\because b_1(\theta) = O(n^{-1}), g(\theta) = O(1) ) \\
&= g(\hat{\theta}) - b_1(\hat{\theta}) + o_p(n^{-1}).
\end{align}

(5)式の推定量の利点は，推定対象$g(\theta)$が正の値であった場合，ブートストラップ推定量も正の値となることが保証されている点です．一方，(4)式の推定量は，それが保証されていません．

ダブルブートストラップを用いた高次のバイアス補正

(4)式は二次漸近不偏推定量でしたが，さらに高次のバイアス補正を行うことも可能です．(4)式の推定量のバイアスを$b_3(\theta)$とすると，
$$
b_3(\theta) = E_\theta \left[ 2g(\hat{\theta}) - E_\hat{\theta}[g(\hat{\theta}^\ast)] \right] - g(\theta), \tag{6}
$$
となります．$2g(\hat{\theta}) - E_\hat{\theta}[g(\hat{\theta}^\ast)]$から$b_3(\theta)$の推定量を引いてあげれば，高次のバイアス補正が行えるわけです．以前の議論と同様に考えてあげると，$b_3(\theta)$のプラグイン推定量$b_3(\hat{\theta})$は以下のように表現できます．

\begin{align}
b_3(\hat{\theta}) &= E_\hat{\theta} \left[ 2g(\hat{\theta}^\ast) - E_{\hat{\theta}^\ast}[ g(\hat{\theta}^{\ast\ast}) ] \right] - g(\hat{\theta}) \\
&= E_\hat{\theta}[ 2g(\hat{\theta}^\ast) ] - E_\hat{\theta}\left[ E_{\hat{\theta}^\ast}[ g(\hat{\theta}^{\ast\ast}) ] \right] - g(\hat{\theta})
\end{align}

形式的には，(6)式の$\theta$は$\hat{\theta}$に，$\hat{\theta}$は$\hat{\theta}^\ast$に，$\hat{\theta}^{\ast}$は$\hat{\theta}^{\ast\ast}$に置き換わっているわけですが，$E_{\hat{\theta}^\ast}[\cdot]$はブートストラップサンプル$y^\ast \sim F_\hat{\theta}$から構成した推定量$\hat{\theta}^\ast$をパラメータに持つ確率分布$F_{\hat{\theta}^\ast}$についての期待値を表します．そして$\hat{\theta}^{\ast\ast}$は，$F_{\hat{\theta}^\ast}$からのブートストラップサンプル$y^{\ast\ast} \sim F_{\hat{\theta}^\ast}$にもとづいた$\theta$の推定量です．$y^{\ast\ast}$はブートストラップサンプル$y^\ast$から構成した$\hat{\theta}^\ast$にもとづいてさらに発生させたブートストラップサンプルなので，ダブルブートストラップサンプルと呼びます．$b$番目のブートストラップサンプルによる$\hat{\theta}^\ast(b)$にもとづいた$C$個のダブルブートストラップサンプルを$y^{\ast\ast}(b,1),\dots,y^{\ast\ast}(b,C) \overset{\mathrm{iid}}{\sim} F_{\hat{\theta}^\ast(b)}$とし，$y^{\ast\ast}(b,c)$による$\theta$の推定量を$\hat{\theta}^{\ast\ast}(b,c)$とすると，
$$
E_\hat{\theta}\left[ E_{\hat{\theta}^\ast}[ g(\hat{\theta}^{\ast\ast}) ] \right] \approx {1 \over BC} \sum_{b=1}^B \sum_{c=1}^C g(\hat{\theta}^{\ast\ast}(b,c))
$$
とモンテカルロ近似できます．これらの結果を用いると，高次のバイアス補正推定量は，
$$
2g(\hat{\theta}) - E_\hat{\theta}[g(\hat{\theta}^\ast)] - b_3(\hat{\theta}) = 3g(\hat{\theta}) - 3E_\hat{\theta}[g(\hat{\theta}^\ast)] + E_\hat{\theta} \left[ E_{\hat{\theta}^\ast}[g(\hat{\theta}^{\ast\ast})] \right]
$$
と導出することができます．この推定量は，高次のバイアス補正はされていますが，計算負荷が高いという欠点があります．また，バイアスは小さくなっても，推定量の効率性を考慮に入れると，望ましい推定量と言えるかはわかりません．

ダブルブートストラップによるバイアス補正は，この他にも，推定対象$g(\theta)$が$\theta$の関数として陽に書けないときの二次漸近不偏推定量の構成でも登場します．こちらについては，次回の記事で，具体例を交えながら紹介したいと思います．

おわりに

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