R: BayesianToolsをためす(2)

前回のつづきです。今度はZero-Infrated Poissonモデル。

Rコード。

library(BayesianTools)

set.seed(123)

data <- read.csv("https://raw.githubusercontent.com/ito4303/naro_toukei/master/example4.csv")

N <- nrow(data)
Y <- data$Num
X <- data$Light

loglik <- function(param) {
  # param[1]: inclusion probability
  # param[2]: intercept of Poisson regression
  # param[3]: slope of Poisson regression
  ll <- 0
  lambda <- exp(param[2] + param[3] * X)
  for (i in 1:N) {
    if (Y[i] > 0) {
      ll <- ll +
        dbinom(1, 1, param[1], log = TRUE) +
        dpois(Y[i], lambda[i], log = TRUE)
    } else {
      p1 <- dbinom(0, 1, param[1])
      p2 <- dbinom(1, 1, param[1]) * dpois(0, lambda[i])
      ll <- ll + log(p1 + p2)
    }
  }
  return(ll)
}

setup <- createBayesianSetup(loglik, prior = NULL,
                             lower = c(0, -100, -100),
                             upper = c(1, 100, 100))
out <- runMCMC(setup, sampler = "DEzs",
               settings = list(iterations = 61500,
                               burnin = 1500,
                               thin = 5))

log(p1 + p2)のところは、logSumExp()をつかうとなぜかうまくいきませんでした。

結果です。

> gelmanDiagnostics(out)
Potential scale reduction factors:

      Point est. Upper C.I.
par 1          1       1.00
par 2          1       1.01
par 3          1       1.00

Multivariate psrf

1
> summary(out)
# # # # # # # # # # # # # # # # # # # # # # # # # 
## MCMC chain summary ## 
# # # # # # # # # # # # # # # # # # # # # # # # # 
 
# MCMC sampler:  DEzs 
# Nr. Chains:  3 
# Iterations per chain:  4001 
# Rejection rate:  0.204 
# Effective sample size:  6734 
# Runtime:  15.16  sec. 
 
# Parameter       MAP      2.5%    median   97.5% 
#  par 1 :        0.402    0.247    0.424    0.646 
#  par 2 :       -0.016   -1.276   -0.071    1.020 
#  par 3 :        0.505   -0.333    0.498    1.290 

## DIC:  4.429086e+40 
## Convergence 
 Gelman Rubin multivariate psrf:   
 
## Correlations 
       par 1  par 2  par 3
par 1  1.000 -0.092  0.014
par 2 -0.092  1.000 -0.309
par 3  0.014 -0.309  1.000

タグ：R MCMC