読了：Vats & Gupta (2012) Rのmcmcseパッケージ

Vats, D., Gupta, K. (2021) An Introduction to Estimating Monte Carlo Standard Errors with R Package mcmcse.
　MCMC標準誤差を出力するRパッケージ mcmcse のvignette. 実戦投入前の儀式として読んだ。

　たとえばですね。目標分布が $$ \left( \begin{array}{c} X_1 \\ X_2 \end{array} \right) \sim N \left( \left( \begin{array}{c} \mu_1 \\ \mu_2 \end{array} \right) , \left( \begin{array}{cc} \omega_1 & \rho \\ \rho & \omega_2 \end{array} \right) \right) $$ ただし\(\omega_1 > 0, \omega_2 > 0, \rho^2 < \omega_1 \omega_2\)だったとするじゃないですか。Gibbsサンプラーでマルコフ連鎖を作ったとしますよね、こんな風に。$$ X_1 | X_2 \sim N \left( \mu_1 + \frac{\rho}{\omega_2} (X_2 – \mu_2), \omega_1 – \frac{\rho^2}{\omega_2} \right) $$ $$ X_2 | X_1 \sim N \left( \mu_2 + \frac{\rho}{\omega_1} (X_1 – \mu_1), \omega_2 – \frac{\rho^2}{\omega_1} \right) $$ これ、mcmcse::BVN_Gibbs()で簡単に生成できますけどね。
　さて。mcmcseパッケージは\(n\)行\(p\)列のMCMC出力を扱います。行が標本、列が要素です。たとえば行\(i\)を\(y_i = (y^{(1)}_i, y^{(2)}_i)\)だとしましょう。
　もしあなたが分布\(F\)の平均\(\mu = E_F y\)に関心があるのなら、その推定量は標本平均\(\mu_n = \frac{1}{n} \sum_{t=1}^n y_t\)ですよね。colMeans()してください。
　もしあなたが二乗和の期待値\(E_F (y^{(1)2}+y^{(2)2})\)に関心があるんなら、その推定量は、\(g(x_1, x_2) = x^2_1 + x^2_2\)として$$ \mu_{g,n} = \frac{1}{n} \sum_{t=1}^n g(y_t)$$ ですよね。apply()してmean()してください。
　このように、モンテカルロ推定値を得るにはbaseパッケージで十分です。

　問題は、モンテカルロ推定値のSEを出したいときです。中心極限定理が成り立てば、$$ \sqrt{n}(\mu_n – E_F y) \rightarrow N_p(0, \Sigma)$$ $$ \sqrt{n} (\mu_{g,n} – E_F(||y||^2) \rightarrow N_p(0, \Sigma_g)$$ です[矢印は分布収束]。この\(\Sigma, \Sigma_g\)をどう出すか。
　mcmcパッケージは、マルコフ連鎖についての中心極限定理が成り立つと仮定して、モンテカルロ推定値のSEを出すパッケージです。

　主な関数は以下の通り。

  • mcse(): \(\Sigma\)が\(1 \times 1\)のとき、\(\mu_n\) SE \( \sqrt{\Sigma/n} \) の一致推定値を出力する。
  • mcse.mat(): \( \Sigma/n \) 対角の平方根の一致推定値を出力する。
  • mcse.multi(): \( \Sigma \) の一致推定値を出力する。
  • mcse.initseq(): initial sequence estimatorを使い、\(\Sigma\)の漸近的に保守的な推定値を出す[←???]。Dai & Jones (2017 J.Multivariate Analysis)をみよ。

　以下の引数があります。

  • method: bm(バッチ平均法), obm(スペクトル分散法), bartlett(修正Bartlett法), tukey(Tukey-Hanning windows)がある。大きな行列ならbmがお勧め。[この辺全然ワカラン]
  • r: lugsailパラメータ。ラグ・ウィンドをどれだけリフトするかを示す。r=1がふつうの場合で、それより大きくすると、\(\Sigma\)の過小評価が取り除かれる代わりに推定量の変動が大きくなる。1,2,3がお勧め。
  • […他にもいろいろ書いてあるけどメモ省略]

　2つのパラメータの同時信頼区間をつくるconfRegion()関数もご用意している。

　ESSについて。
　サンプリングする前に、minESS()関数で必要な最小ESSを調べておくと良い。変数数、alpha, eps (tolerance level)をいれると最小ESSが出る。変数数、alpha, ESSをいれてepsを出すこともできる。[へー]
　ess()で要素ごとのESSが出せる。multiESS()で多変量ESSが出せる。

　いちおう描画機能もある。estvssample()で、横軸に標本サイズ、縦軸に推定値をとったチャートを描く[??? running meanプロットのこと?]
　云々。