読者です 読者をやめる 読者になる 読者になる

廿TT

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

Rcpp で独立メトロポリスヘイスティングス

独立メトロポリス・ヘイスティングス法を用いたベイズ推測の簡単な例題 - 廿TT でやったのと同じことを Rcpp で書いてみた。
ハローワールド。

C++ のコードはこう。

#include <Rcpp.h>
using namespace Rcpp;
double lik(double lambda, NumericVector dat) {
  double sumdat;
  sumdat = sum(dat);
  return(dat.length()*log(lambda)-sumdat*lambda);
}
// [[Rcpp::export]]
NumericVector MHexp(int N, NumericVector dat) {
  NumericVector chain(N);
  double a;
  double r;
  RNGScope scope;
  chain[0] =R::runif(0, 10);
  for(int i = 0; i < N; i++) {
    a = R::runif(0, 10);
    r = lik(a,dat)-lik(chain[i-1],dat);
    if(R::runif(0,1) < exp(r)){
      chain[i] = a;
    }else{
      chain[i] = chain[i-1];
    }
  }
  return(chain);
}

R からはこうやって使う。

library(Rcpp)
sourceCpp('MHexp.cpp')
dat <- rexp(100,2)
system.time({out <-MHexp(100000,dat)})
#   ユーザ   システム       経過  
#      0.03       0.00       0.03 
plot(out,type="l")

f:id:abrahamcow:20170329065503p:plain