廿TT

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

Rcpp で PMMH(パーティクルマージナルメトロポリス・ヘイスティングス)

今日の川柳

R によるすごくかんたんなパーティクルフィルタの実装例 - 廿TT では分散パラメータを既知としていますが、分散パラメータも推定したい。

パーティクルフィルタの重みから周辺尤度を出してランダム・ウォーク・メトロポリスヘイスティングスでサンプリングします。

これらのサイトを参考にしました。

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;

// [[Rcpp::export]]
List locallevel(NumericVector y, int N_particles, double sigma, double tau) {
  int tmax = y.length();
  NumericMatrix xx(N_particles,tmax+1);
  NumericVector x(N_particles);
  NumericVector w(N_particles);
  IntegerVector ind(N_particles);
  IntegerVector int_seq=seq(0,N_particles-1);
  double ll=0;
  xx(_,0) = rep(y[0],N_particles);
    for(int t=0; t<tmax; t++){
      x = xx(_,t) + rnorm(N_particles,0,tau);
      w = exp(-pow(rep(y[t],N_particles)-x,2)/(2*pow(sigma,2)))/sigma;
      ll = ll+log(mean(w));
      x = RcppArmadillo::sample(x,N_particles,true,w/sum(w));
      xx(_,t+1) = x;
    }
  return List::create(Named("state")=xx,_["logLik"]=ll);
}

// [[Rcpp::export]]
NumericMatrix PMMH(NumericVector y, int iter, NumericVector ini, int thin){
  List Filt1 = locallevel(y, 100, exp(ini[0]), exp(ini[1]));
  double ll = Filt1["logLik"];
  double llprop = 0;
  NumericMatrix thetamat(iter,2);
  NumericVector theta0=ini;
  NumericVector thetaprop(2);
  bool accept;
  double u=0;
  for (int i=0; i<iter; i++) {
    for (int j=0; j<thin; j++) {
      thetaprop = rnorm(2, 0, 0.1);
      thetaprop += theta0;
      Filt1 = locallevel(y, 100, exp(thetaprop[0]), exp(thetaprop[1]));
      llprop = Filt1["logLik"];
      u = R::runif(0,1);
      accept = log(u) < (llprop - ll);
      if (accept) {
        theta0=thetaprop;
        ll=llprop;
      }
    }
    thetamat(i,_) = theta0;
  }
  return thetamat;
}

トレースプロットを見ると、収束は大丈夫そうです。

library(Rcpp)
library(tidyverse)
library(parallel)
sourceCpp("PF.cpp")
ini = rep(log(sd(Nile)),2)
system.time({
samples <-mclapply(1:4,function(i){PMMH(Nile,2000,ini,5)},mc.cores = detectCores())
})
sigmas <-sapply(samples,function(x)x[,1])
taus <-sapply(samples,function(x)x[,2])
col4<-c("royalblue","salmon","forestgreen","orange2")
matplot(sigmas,type="l",col = col4, lty=1)
matplot(taus,type="l",col = col4, lty=1)

f:id:abrahamcow:20180315040220p:plain

f:id:abrahamcow:20180315040229p:plain

得られた分散パラメータの点推定値を使ってフィルタリング。

Filt1 <-locallevel(y=Nile,N_particles = 1000,
                   sigma=exp(mean(sigmas[-c(1:500),])),
                   tau=exp(mean(taus[-c(1:500),])))

statemean <-apply(Filt1$state,2,mean)[-1]

plot(Nile)
lines(start(Nile)[1]:end(Nile)[1],statemean,col="royalblue",lwd=2)

f:id:abrahamcow:20180315040334p:plain

dlm でも同じモデルをやってみると、結果はほぼ一致することがわかります。

library(dlm)
N<-length(Nile)
build1 <- function(theta){
  dlmModPoly(order=1, dV=exp(theta[1]), dW=exp(theta[2]),m0=Nile[1],C0 = 0)
}
fit1 <- dlmMLE(Nile, parm=c(1, 1), build1)
fit1$par
modNile <- build1(fit1$par)
NileFilt <- dlmFilter(Nile, modNile)
lines(start(Nile)[1]:end(Nile)[1],dropFirst(NileFilt$m),col="red",lwd=2)

f:id:abrahamcow:20180315040445p:plain

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

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