Particle Markov chain Monte Carlo methods (PMCMC)について

Particle Markov chain Monte Carlo methods (PMCMC)

時系列の推定とモデル(のパラメータ)の推定においてParticle filter(SMC)とMCMCを組み合わせた手法があり、その分かりやすい解説としてParticle Markov chain Monte Carlo methods(pdf)というドキュメントを読んだのでその内容について記載します。

構成は
1. Introduction
2. アルゴリズム概要
3. 適応例(Lévy-driven stochastic volatility model,簡単な非線形モデル)
4. PMCMCの一般的定式化
5. Discussion
Appendix
Reference
色々な先生方による講評(半分以上を占める)

となっていますが、その２章の内容が以下になります。

普通のparticle filter(SMC)

SMCは

Resampling(前のステップで計算された重みヒストグラムを用いた粒子の取り直し)
Propagation(状態空間モデルの時間発展)
Weighting(粒子による重みヒストグラムの計算)

の手順の繰り返しで構成されています。
時点nの粒子のkの親となる粒子のindexを\( A_n^k \) k番目の粒子に至る"家系"(linage) \(X_{1:T}^k \)のindexを\( B_{1:T}^k \)とすると

Resampling

\( A^k_{n-1} \sim F(\cdot |W_n-1) \)

Propagation

\( X_n^k \sim q(\cdot |y_n, X_{n-1}^{A^k_{n-1}}) \)
\( X^k_{1:n}=(X_{1:n}^{A_{n-1}^k},X_n^k) \)

Weighting

\( w_n(X^k_{1:n})=\frac{f_\theta(X_n^k|X_{n-1}^{A_{n-1}^k}) g_\theta(y_n|X_n^k) }{q_\theta(X_n^k|y_n,X_{n-1}^{A_{n-1}^k})} \)
\( W_n^k=\frac{w_n(X^k_{1:n}) }{\sum_{i=1}^N w_n(X_{1:n}^i ) } \)
と書けます。
数値計算上はF,qからのサンプリングはヒストグラムからの抽出、状態空間モデルの式の粒子ごとの計算に相当します。
上記の手順では全時系列の経路をコピーすることになりますが、各ステップでの重みw_nの計算には直前のn-1の情報のみを使います。多次元、あるいは形状が複雑で積分で形状を計算することの難しい事後分布が
\( \hat{p}_\theta(x_{1:n}|y_{1:n})=\sum_{k=1}^N W_n^k \delta_{X^k_{1:n}} (x_{1:n}) \)
で近似されることになります。

SMCはオンラインでの時系列推定手法であるので後になるにつれて精度が高まっていく特徴があります。
経路全体の事後分布、全時系列データを用いたパラメータの事後分布を計算するには経路全体を対象としたサンプリング、推定が必要になります。そこで以下のような方法が挙げられています。

Particle Markov Chain Monte Carlo (PMCMC)

Particle independent Metropolis–Hastings sampler

パラメータ推定と時系列の推定の方法としてはSMCとEMアルゴリズムを用いて最尤値を求める方法*1あるいは
パラメータが時間的に変化する場合には自己組織型状態空間モデル*2が知られていますが、ここではMetropolis–Hastings(MH)法を用いて事後分布と時系列を推定する方法が説明されています。
時系列\(X_{1:T}\) のsamlingにはSMCで得られた事後分布の近似\( \hat{p}_\theta(\cdot|y_{1:T}) \)を用い、Metropolis–Hastings法で棄却サンプリングを行います。

iteration:
\( \hat{p}^*(\cdot|y_{1:T})=SMC() \\
X^*_{1:T} \sim \hat{p}(\cdot|y_{1:T}) \)
if( min(1, \( \frac{\hat{p}_\theta^*(\cdot|y_{1:T}) }{\hat{p}_\theta(\cdot|y_{1:T})} \) ) >rand() ):
\(X_{1:T}(i)=X^*_{1:T} \)
\(\hat{p}(\cdot|y_{1:T})(i)=\hat{p}^*(\cdot|y_{1:T}) \)
else:
\(X_{1:T}(i)=X_{1:T}(i-1) \)
\(\hat{p}(\cdot|y_{1:T})(i)=\hat{p}(\cdot|y_{1:T})(i-1) \)

Particle marginal Metropolis–Hastings sample(PMMH)

パラメータθの推定を行うのに時系列\(X_{1:T}\)とθを同時にupdateする方法です。
周辺尤度(marginal likehood) \(p_\theta(y_{1:T})\)を用いることからこの名前があるようです。

iteration:
\( \theta^* \sim q(\cdot|\theta(i-1) ) \)
\( \hat{p}_{\theta^*} (\cdot|y_{1:T}) = SMC(\theta^*, y_{1:T}) \)
\( X^*_{1:T} \sim \hat{p}_{\theta^*} (\cdot|y_{1:T}) \)
if( min(1, \( \frac{ \hat{p}_{\theta^*} (\cdot|y_{1:T}) p(\theta^*) }{ \hat{p}_{\theta(i-1)} (\cdot|y_{1:T}) p(\theta(i-1)) } \frac{ q(\theta(i-1)|\theta^*) }{ q(\theta^*|\theta (i-1) )}\) >rand() ):
\(\theta(i)=\theta^* \)
\(X_{1:T}(i)=X_{1:T}^* \)
\(\hat{p}_{\theta(i)}(y_{1:T})=\hat{p}_{\theta^*}(y_{1:T}) \)
else:
\(\theta(i)=\theta(i-1) \)
\(X_{1:T}(i)=X_{1:T}(i-1) \)
\(\hat{p}_{\theta(i)}(y_{1:T})=\hat{p}_{\theta(i-1)}(y_{1:T}) \)

Conditional Particle Filter

前の処理で推定された時系列
\( X_{1:T}=(X_1^{B_1}, X_2^{B_2}, \cdots , X_T^{B_T} ) \)
を1つだけとっておいてそれを種のようにして用いてSMCを行う手法です。これによって以下のParticle Gibbs samplerで用いる近似事後分布\(\hat{p}_\theta(\cdot|y_{1:T}) \)を得ます。

n=1:
for \(m \neq B_1^{k}: X_1^m \sim q_1() \)
compute \(w_1(X_1^k), W_n^k \propto w_1(X_1^k) \)
for n in 2 … T:
\(A^k_{n-1} \sim F(\cdot |W_{n-1}) (k \neq B_n) \)
\(X^k_n \sim q(\cdot |y_n,X_{n-1}^{A^k_{n-1}}) (k \neq B_n) \)
compute \(w_n(X_{1:n}^k) , W_n^k \propto w_n(X^k_{1:n}) \)

Aはindexなのでそのサンプリングは単に乱数を選ぶだけになります。

Particle Gibbs sampler

Conditional Particle Filterで得られた近似事後分布を用いてGibbs samplingを行います。前のiterationで作られた\(X_{1:T}\)を繰り返し用いるのででMetropolis–Hastings samplingは行う必要がありません。
\(X_{1:T}, y_{1:T}\)が得られればθのサンプリングは難しくないということですが、一般的な場合を考えるとそう簡単ではなさそうです。ノイズの分散などがθの場合は単純な分布からサンプリングするだけなので簡単そうです。

iteration:
\( \theta(i) \sim p{\cdot|y_{1:T},X_{1:T}(i-1) } \)
\( \hat{p}_{\theta(i)} (\cdot|y_{1:T}) = SMC_{conditional}(\theta(i), X_{1:T}(i-1), y_{1:T}) \)
\( X_{1:T}(i) \sim \hat{p}_{\theta(i)} (\cdot|y_{1:T}) \)

家系\(B_{1:T}(i)\)は\(X_{1:T}(i)\)から自ずと定まります。
粒子を再利用して効率的に経路を得る手法(論文4.6)やauxiliary particle filterで退化を防ぐ方法などがあるそうです。