廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Rcpp を用いたギブスサンプリングのかんたんな例題

参考文献

Rによるモンテカルロ法入門

Rによるモンテカルロ法入門

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

ギブスサンプリングのなんたるかについてはこれらの本をご覧ください.

設定

X_i が独立に平均 \mu, 分散 \sigma^2正規分布に従うとする.

\mu の事前分布に平均 \mu_0, 分散 \sigma^2_0正規分布, \sigma^2 の事前分布にシェイプパラメータ \alpha_0/2, レートパラメータ \beta_0/2 の逆ガンマ分布を仮定する.

\mu, \sigma の事後分布は以下に比例する.


 \displaystyle \frac{1}{(\sqrt{2\pi \sigma^2})^n} \exp\left(-\frac{\sum(x_i-\mu)^2}{2\sigma^2}\right) \times\\ 
\displaystyle \quad \frac{1}{\sqrt{2\pi \sigma^2_0}} \exp\left(-\frac{(\mu-\mu_0)^2}{2\sigma^2_0}\right) \times \\
\displaystyle \quad \frac{(\beta_0/2)^{\alpha_0/2}}{\Gamma(\alpha_0/2)} (\sigma^2)^{-\frac{\alpha_0}{2}+1}\exp\left( -\frac{\beta_0}{2\sigma^2}\right)\\
\displaystyle \propto (\sigma^2)^{-\frac{\alpha_0+n}{2}+1}\exp\left(-\frac{\sum(x_i-\mu)^2}{2\sigma^2}-\frac{(\mu-\mu_0)^2}{2\sigma^2_0}-\frac{\beta_0}{2\sigma^2} \right)

\mu の条件付き事後密度は, 以下に比例する.


 \displaystyle \exp\left(-\frac{(\mu-\mu_0)^2}{2\sigma^2_0}\right)\exp\left(-\frac{n(\bar{x} -\mu)^2}{2\sigma^2}\right) \\
\displaystyle \propto_{\mu} \exp\left( (\sigma^{-2}_{0}+n\sigma^{-2})(\mu-\frac{\sigma^{-2}_{0}\mu_0+n\sigma^2\bar{x}}{\sigma^{-2}_{0}+n\sigma^{-2}}) \right)

この密度は正規分布に対応する.

\sigma^2 の条件付き事後密度は, 以下に比例する.


 \displaystyle (\sigma^2)^{-\frac{\alpha_0+n}{2}+1}\exp\left(-\frac{\sum(x_i-\mu)^2 + \beta_0}{2\sigma^2}\right)

この密度は逆ガンマ分布に対応する.

コード

Rcpp のコード:NormalGibbs.cpp

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
List GibbsNormal(NumericVector x, double init,
                 double mu0, double sigma20, double a, double b,
                 int iter) {
  NumericVector mu(iter);
  NumericVector sigma2(iter);
  int n = x.length();
  double xbar = mean(x);
  sigma2[0] = init;
  double B;
  B = 1/sigma20+n/sigma2[0];
  mu[0] = R::rnorm(((1/sigma20)*mu0+n*(1/sigma2[0])*xbar)/B,1/sqrt(B));
  double ra1;
  double sh1 = a+n;
  for(int i=1; i<iter; i++){
    B = 1/sigma20+n/sigma2[i-1];
    mu[i] = R::rnorm(((1/sigma20)*mu0+n*(1/sigma2[i-1])*xbar)/B,1/sqrt(B));
    ra1 = 0;
    for(int j=0; j<n ;j++){
      ra1 = ra1+pow(x[j]-mu[i],2);
    }
    ra1 = ra1 + b;
    sigma2[i] = 1/R::rgamma(sh1,1/ra1);
  }
  return List::create(Named("mean") = mu , _["variance"]=sigma2);
}

これを R から動かしてみる.

初期値 init は \sigma^2 の初期値のみを要求する.

事前分布を無情報的にするため, \sigma^2_0,  \alpha_0, \beta_0 はそれぞれ 10000, 0.02, 0,02 と分散が大きくなるように選んだ.

library(Rcpp)
sourceCpp("NormalGibbs.cpp")
x <- c(91,504,557,609,693,727,764,803,857,929,970,1043,1089,1195,1384,1713)
logx <- log(x)

fit1 <-GibbsNormal(x = logx, init = 1,
                   mu0 = 0, sigma20 = 10000,
                   a = 0.02, b = 0.02,iter = 2000)
plot(fit1$mean,type="l")
plot(fit1$variance,type="l")

\mu, \sigma^2 それぞれのトレースプロットはこんな感じ.

f:id:abrahamcow:20170916033726p:plain

f:id:abrahamcow:20170916033811p:plain

信用区間ベイズ信頼区間)と信頼区間を並べてみたくなったので, 正規分布の分散の信頼区間を出す関数を以下に定義する.

varCI <- function(x){
  df <- length(x) - 1
  varx <- var(x)
  lower = varx * df / qchisq(0.05/2, df, lower.tail = FALSE)
  upper = varx * df / qchisq(1 - 0.05/2, df, lower.tail = FALSE)
  c(lower = lower, upper = upper)
}


事後期待値は最尤推定値とだいぶ近くなる.
分散の信用区間は信頼区間よりもだいぶ狭めに出るようだ. コード間違ってるかな? 不安…

> mean(logx)
[1] 6.631066
> mean(fit1$mean[-c(1:200)])
[1] 6.632218
> var(logx)
[1] 0.4267377
> mean(fit1$variance[-c(1:200)])
[1] 0.4586105
> quantile(fit1$mean[-c(1:200)],c(0.025,0.975))
    2.5%    97.5% 
6.282746 6.958671 
> t.test(logx)$conf.int[1]
[1] 6.282972
> quantile(fit1$variance[-c(1:200)],c(0.025,0.975))
     2.5%     97.5% 
0.2722280 0.7618296 
> varCI(logx)
    lower     upper 
0.2328643 1.0221854 

質問

for(int j=0; j<n ;j++){
  ra1 = ra1+pow(x[j]-mu[i],2);
}

の部分をもっとシンプルに書けないでしょうか?