廿TT

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

パーティクルフィルタとMCMCによる離散時間SIRモデルのパラメータ推定

今日の川柳

  1. SIR モデルと非定常ポアソン過程 - 廿TT
  2. Rcpp で PMMH(パーティクルマージナルメトロポリス・ヘイスティングス) - 廿TT
  3. ゼロの多いコンバージョンの状態空間モデル - 廿TT

これらのエントリの合わせ技です。

考えているモデルは
SIRモデルからはじめる微分方程式と離散時間確率過程(後編) - StatModeling Memorandum
とほぼ同じですが、観測モデルにはポアソン分布を仮定しています。

データは MERS2015.csv · GitHub に置いてあります。

f:id:abrahamcow:20180615033145p:plain

青がS、黄色がI、緑がRです。点が観測値です。

R のコードを貼ります。

library(tcltk)
dat <-read.csv(file = "MERS2015.csv",stringsAsFactors = FALSE)
dR <-diff(c(0,dat$Recovered))
Tim <-length(dR)
softmax <- function(x){
  tmp <- max(x,na.rm = TRUE)
  out <-exp(x-tmp)/sum(exp(x-tmp),na.rm = TRUE)
  out[is.na(out)] <- 0
  return(out)
}
logsumexp <- function(x,y){
  x + log(1+exp(y-x))
}
logmeanexp <- function(x){
  tmp <- max(x)
  tmp + log(sum(exp(x-tmp))) - log(length(x))
}
PF<- function(beta,gamma,smoothing=FALSE,np = 2000){
  dItoR <- dStoI <- S <- I <- R <-matrix(0,Tim,np)
  wt <-matrix(0,Tim,np)
  S[1,] <-rpois(np,rgamma(np,1,1/10))+sum(dR)
  I[1,] <-rpois(np,1)+1
  N <- S[1,] + I[1,]
  ll <-0 
  for(i in 2:Tim){
    dStoI[i,] <- rbinom(np,S[i-1,],beta*I[i-1,]/N)
    S[i,] <- S[i-1,]-dStoI[i,]
    dItoR[i,] <- rbinom(np,I[i-1,],gamma)
    I[i,] <- I[i-1,]+dStoI[i,] - dItoR[i,]
    R[i,] <- R[i-1,]+dItoR[i,]
    wt[i-1,] <-dpois(dR[i-1],dItoR[i,],log = TRUE)
    ind <-sample.int(np, replace = TRUE, prob = softmax(wt[i-1,]))
    S[i,] <- S[i,ind]
    I[i,] <- I[i,ind]
    R[i,] <- R[i,ind]
    dStoI[i,] <- dStoI[i,ind]
    dItoR[i,] <- dItoR[i,ind]
    N <- N[ind]
    ll <- ll + logmeanexp(wt[i-1,])
  }
  if(smoothing){
    for(i in (Tim-1):1){
      wt2 <- logsumexp(wt[i,],dbinom(dStoI[i+1,],S[i,],beta*I[i,]/N,log=TRUE)+
                         dbinom(dItoR[i+1,],I[i,],gamma,log=TRUE))
      ind <-sample.int(np, replace = TRUE, prob = softmax(wt2))
      S[i,] <- S[i,ind]
      I[i,] <- I[i,ind]
      R[i,] <- R[i,ind]
      dStoI[i,] <- dStoI[i,ind]
      dItoR[i,] <- dItoR[i,ind]
      N <- N[ind]
    } 
  }
  list(S=S,I=I,R=R,dStoI=dStoI,dItoR=dItoR,N=N,ll=ll)
}

iter <- 5000
thetamat <- matrix(0,iter,2)
theta0 <- c(0,-1.5)
PFout <-PF(plogis(theta0[1]),plogis(theta0[2]))
ll = PFout$ll
pb <- txtProgressBar(min = 1, max = iter, style = 3)
set.seed(1192)
system.time({
for (i in 1:iter) {
  thetaprop = theta0+rnorm(2, 0, 0.1)
  PFout <-PF(plogis(thetaprop[1]),plogis(thetaprop[2]))
  llprop = PFout$ll
  u = runif(1)
  accept = log(u) < (llprop - ll)
  if (accept) {
    theta0=thetaprop
    ll=llprop
  }
  setTxtProgressBar(pb, i)
  thetamat[i,] = theta0
}
})
# user  system elapsed 
# 360.037  32.151 394.774 

betahat <-mean(plogis(thetamat[,1]))
gammahat <-mean(plogis(thetamat[,2]))
PFout <-PF(betahat,gammahat,smoothing = TRUE)
Shat <- rowMeans(PFout$S)
Ihat <- rowMeans(PFout$I)
Rhat <- rowMeans(PFout$R)
plot(as.Date(dat$Date),Shat,ylim = c(0,max(Shat)),type = "l",col="royalblue",lwd=2,
     xlab="",ylab="population")
lines(as.Date(dat$Date),Ihat,col="orange2",lwd=2)
lines(as.Date(dat$Date),Rhat,col="forestgreen",lwd=2)
points(as.Date(dat$Date),dat$Recovered)

トレースプロットはこんな感じ。

f:id:abrahamcow:20180615033519p:plain

f:id:abrahamcow:20180615033501p:plain

plot(thetamat[,1],type="l", ylab="beta")
plot(thetamat[,2],type="l", ylab="gamma")
> print(betahat)
[1] 0.57834
> print(gammahat)
[1] 0.1755242

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

現象数理学入門

現象数理学入門

不明点あれば聞いてください。

雑感:推定結果はSの初期状態の分布の設定にけっこう依存する。