東京大学・株式会社Nospareの菅澤です．
今回は 重み付き尤度ブートストラップ（WLB; weighted likelihood bootstrap） によって事後分布を近似計算する方法について紹介します．
はじめに
ブートストラップ法はデータからのリサンプリングを繰り返すことで統計量のバラつきを数値的に見積もることが可能な方法として知られています．
今，$X=(X_1,\ldots,X_n)$を観測データとし，統計量（平均や分散など）を$T(X)$としておきます．Efron (1979)によって提案されたノンパラメトリックブートストラップ法は，$X$からの復元抽出によって得られたブートストラップ標本$X_{\ast}$を用いて統計量$T(X_{\ast})$を計算するステップを繰り返すことで，$T(X)$の標本分布を${T(X_{\ast}^{(1)}),\ldots,T(X_{\ast}^{(B)})}$で近似する方法です．
また，この方法をベイズ的に解釈したベイジアンブートストラップと呼ばれる考えがRubin (1981)によって提案されています．端的に説明すると，ノンパラメトリックブートストラップが標本を多項分布的にリサンプリングするのに対し，ベイジアンブートストラップは各標本の重みを事後分布として与えてリサンプリングする方法です．
今回紹介するWLBはNewton and Raftery (1994)によって提案された方法で，ベイジアンブートストラップのある種の拡張となる方法です．
WLB (weighted likelihood bootstrap)とは
データ$X=(x_1,\ldots,x_n)$の個々の$x_i$が独立に確率分布$f(x_i;\theta)$に従うとします．$\theta$に対する事前分布を$\pi(\theta)$とすると，$\theta$の事後分布は
\pi(\theta \mid X)\propto \pi(\theta)L(\theta), \ \ \ \ \ L(\theta)\equiv \prod_{i=1}^n f(x_i;\theta)
で与えられます．この事後分布を近似する代表的な方法はマルコフ連鎖モンテカルロ(MCMC)法で，$\pi(\theta|X)$から間接的な乱数生成を行うことで$\theta$のベイズ推測を行うことができます．一方で，モデルの形によっては事後分布からの乱数生成が（たとえMCMCを使ったとしても）簡単じゃないケースも多くあります．WLBは「事後分布からの乱数生成は難しいが，重み付き尤度の最大化は簡単」な状況で事後分布を近似するアルゴリズムです．
具体的には，以下の関数の最大化を考えます．
\tilde{L}(\theta)=\prod_{i=1}^n f(x_i;\theta)^{w_i}.
ここで，$w=(w_1,\ldots,w_n)$はランダムな重みです．$\tilde{L}(\theta)$の最大化は重み付き対数尤度$\sum_{i=1}^n w_i \log f(x_i;\theta)$の最大化と等価です．
重み付き尤度を最大化する$\theta$を$\tilde{\theta}(w)$と表すと，WLBは$B$回生成したランダムな重み（$w^{(1)},\ldots,w^{(B)}$）に対して，$\{\tilde{\theta}(w^{(1)}),\ldots,\tilde{\theta}(w^{(B)})\}$によって事後分布を近似します．この際に，$B$回重み付き関数を最大化する必要があることに注意が必要です．したがって，このアルゴリズムが計算速度の意味で効率的か否かはモデルの形に依って大きく異なります．
余談ですが，近年は大規模なデータや高次元のデータに対してベイズ推測を行う際に「MCMCでのサンプリングは計算コストがかかるが，尤度関数の最大化（または事後モードの計算）は容易」というケースがあり，WLBを発展させたアルゴリズムによる近似ベイズ推測のアルゴリズム開発も行われています．この話は一般化ベイズ法とも深く関わってくる話題でして，その点については今後別の記事で紹介しようと思います．
WLBを実行するためには重み$w$を生成する確率分布を決める必要がありますが，Newton and Raftery (1994)ではディリクレ分布${\rm Dir}(1,\ldots,1)$を考えています．これは各$w_i$を指数分布の乱数をスケーリングした形$w_i\propto Y_i\sim {\rm Exp}(1)$と取ることと同値になります．（これはディリクレ分布が複数の独立なガンマ分布をスケーリングすることで導出できることに由来します．）
重みを$w\sim {\rm Dir}(1,\ldots,1)$とすることの一般化として$w_i\propto Y_i^{\alpha}$とした重みの設定もNewton and Raftery (1994)で考えられています．$\alpha>1$のとき，$w$の分布がディリクレ分布よりもover-dispersedな分布になります．
重点ウェイトによる調整
上記のWLBのアルゴリズムでは事前分布$\pi(\theta)$の情報が考慮されていません．$n\to\infty$では基本的に事前分布の影響が無視されるため，尤度パートだけ考えたWLBでも事後分布の良い近似になります．一方で，そこまで大きくない$n$のもとでの近似結果はあまり良くない可能性があります．事前分布の影響も考慮したより精度の高い近似を得るためには，各$\tilde{\theta}_b\equiv \tilde{\theta}(w^{(b)})\ (b=1,\ldots,B)$に対して重点ウェイト
u_b \propto \frac{\pi(\tilde{\theta}_b)L(\tilde{\theta}_b)}{\hat{g}(\tilde{\theta}_b)}
を与えれば解決できます．ここで$\hat{g}$は$\tilde{\theta}_b$の生成分布であり，生成されたサンプルを用いて推定する必要があります．
近年では，そもそものWLBのアルゴリズム内で事前分布に重みを与えて近似事後分布を導出する方法なども提案されています．（こちらについても今後の記事で紹介していきたいと思います．）
理論的な性質
$\hat{\theta}$を尤度関数$L(\theta)$を最大にする最尤推定量，$I_n(\theta)$をフィッシャー情報行列とします．
このとき，モデル$f(x;\theta)$に関する適当な正則条件のもと，WLBによる事後分布$\sqrt{nI_n(\hat{\theta})}(\tilde{\theta}-\hat{\theta})$が漸近的に標準正規分布に従います．また，通常の事後分布に対して$\sqrt{nI_n(\hat{\theta})}(\theta-\hat{\theta})$が漸近的に標準正規分布に従うことを考慮すると，WLBの事後分布と通常の事後分布がこのスケールで同一の挙動を持つことがわかります．（あくまで$n^{-1/2}$のスケールでの挙動が一致するだけで，より高次のスケールでの挙動は必ずしも一致しません．）
この結果より，$n$がある程度大きければWLBによって生成した$\theta$のサンプルによって近似的な事後推論が実行できることがわかります．
ベイジアンブートストラップとの関係
データ$x_1,\ldots,x_n$が$\theta=(\theta_1,\ldots,\theta_K)$の確率で$K$個の異なる値（一般性を失うことなく$\{1,\ldots,K\}$とします）をとる多項分布に従う状況を考えます．また，$y_j\ (j=1,\ldots,K)$を$j$番目のカテゴリに属するデータの個数とします．すなわち$y_j=\sum_{i=1}^n I(x_i=j)$です．
このモデルはデータの生成過程をノンパラメトリックに表現したモデルであると捉えることができます．ブートストラップ標本の生成をこの多項モデルにおける$\theta$の事後シミュレーションに置き換えたものがベイジアンブートストラップ(Rubin, 1981)です．
尤度関数および重み付き尤度関数は
L(\theta)=\prod_{j=1}^K\theta_j^{y_j}, \ \ \ \
\tilde{L}(\theta)=\prod_{j=1}^K\theta_j^{n\gamma_j}
となります．ここで，$\gamma_j=\sum_{i=1}^n w_iI(x_i=j)$はWLBの重みを各カテゴリごとに和をとったものです．ディリクレ分布に従う確率ベクトルの部分和のベクトルはまたディリクレ分布になる性質を用いると，$w\sim {\rm Dir}(1,\ldots,1)$のもとで$(\gamma_1,\ldots,\gamma_K)\sim {\rm Dir}(y_1,\ldots,y_K)$が成り立ちます．また，重み付き尤度$\tilde{L}(\theta)$を$\theta$について最大化すると$\tilde{\theta}_j=\gamma_j$となるため，WLBで与えられる$\theta$の分布は${\rm Dir}(y_1,\ldots,y_K)$であることがわかります．この分布は，$\theta$に対して非正則な事前分布${\rm Dir}(0,\ldots,0)$を用いた際の事後分布に相当します．したがって，多項モデルではWLBとベイジアンブートストラップが一致します．WLBは一般のモデル$f(x;\theta)$に対して適用できるため，その意味でWLBがベイジアンブートストラップの拡張であるとみなすことができます．
おわりに
今回はweighted likelihood boostrap (WLB) の方法について紹介しました．WLBを提案したNewton and Raftery (1994)では仮定したモデルが「正しい」ことを想定していますが，近年は一般化ベイズ法との関連で誤特定された状況でのWLB法(Lyddon et al., 2019)も開発されています．また，spike-and-slab事前分布などMCMCでの計算が非常に困難な事後分布をWLBのアイデアで近似する研究(Nie and Rovkova, 2022)などもあります．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
