廿TT

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

ガンマ・ポアソンの状態空間モデル

PMMH(パーティクルマージナルメトロポリスヘイスティングス; Rcpp で PMMH(パーティクルマージナルメトロポリス・ヘイスティングス) - 廿TT)でもう少し遊ぶ。

パーティクルフィルタは正規分布以外の分布でも使えるということなので、未観測の状態変数がガンマ分布、観測モデルがポアソン分布のモデルを試してみる。

システムモデル:
\lambda_i \sim {\rm Gamma}(\alpha, \alpha/\lambda_{i-1})

観測モデル:
y_i \sim {\rm Poisson}(\lambda_i)

なんの役に立つのかはわかりませんけれども。

\alpha は未知パラメータなので PMMH で推定する。

解析対称はイギリスの炭鉱事故の発生件数のデータで、3. Tutorial — PyMC 2.3.6 documentation から取得した。

library(parallel)
y <- c(4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
         3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
         2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
         1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
         0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
         3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
         0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1)

GammaFilter <-function(y,alpha,N_particles){
Tmax <- length(y)
lambdaFilt <-matrix(NA,Tmax,N_particles)
lambdatmp <-rgamma(N_particles,y[1]+0.001,1+0.001)
ll <- 0
for(t in 1:Tmax){
  lambdatmp <-rgamma(N_particles,alpha,alpha/lambdatmp)
  w <-dpois(y[t],lambdatmp)
  ll <- ll+log(mean(w))
  w <- w/sum(w)
  lambdatmp <-sample(lambdatmp,N_particles,replace = TRUE,prob = w)
  lambdaFilt[t,] <- lambdatmp
}
list(lambdaFilt=lambdaFilt,loglik=ll)
}
llfunc <- function(alpha,N_particles){
  Filt1 <-GammaFilter(y,alpha,N_particles)
  Filt1$loglik
}

PMMH <- function(iter=1000){
alpha <- runif(1) + rnorm(1,0,1)
ll <- llfunc(exp(alpha),500)
alphavec <- numeric(iter)
alphavec[1] <- alpha
for(i in 2:iter){
  alpha <- alphavec[i-1] + rnorm(1,0,1)
  propll <- llfunc(exp(alpha),500)
  if(log(runif(1))<propll-ll){
    alphavec[i] <- alpha
    ll <- propll
  }else{
    alphavec[i] <- alphavec[i-1]
  }
  }
return(alphavec)
}

samples <-mclapply(1:3,function(i)PMMH(),mc.cores = 3)
samples <-simplify2array(samples)
matplot(samples,type="l",lty=1,col=c("royalblue","orange","forestgreen"))

こんな感じで  \log(\alpha) が求まった。

f:id:abrahamcow:20180317181608p:plain

求まったパラメータの点推定値でフィルタリングするとこんな感じ。

Filt <- GammaFilter(y,exp(mean(samples[-c(1:200),])),1000)
plot(y,type="b")
lines(apply(Filt$lambdaFilt , 1 , mean),col="royalblue", lwd=2)

f:id:abrahamcow:20180317181735p:plain

以上です。

Rによるベイジアン動的線形モデル (統計ライブラリー)

Rによるベイジアン動的線形モデル (統計ライブラリー)