MCMCmetrop1Rに事前分布を入れてみる

もう4年ほど前にMCMCmetrop1Rを試し始めたのだが、そのころは事前分布はあまり考えていなかった。というか、MCMCmetrop1Rのhelpのサンプルコードにもないのだが、これは要するに一様分布を仮定しているということになるのだろうか。

というようなことをふと思いついたので、ではどうやって事前分布を入れるのか、と試してみた。

試行錯誤して、MCMCmetrop1Rに渡す関数の返り値をこういうふうにした。beta[1] = mean, beta[2] = sdである。

    l <- sum(log(dnorm(x, mean = beta[1], sd = beta[2]))) +
         log(dnorm(beta[1],
                   mean = prior.m.mean, sd = prior.m.sd)) +
         log(dinvgamma(beta[2]^2,
                       shape = prior.s.shape, scale = prior.s.scale))

実際に走らせてみる。

library(MCMCpack)

set.seed(1)
m <- 10
s <- 3
x <- rnorm(10, mean = m, sd = s)

# priors
prior.m.mean <- 0
prior.m.sd <- 10
prior.s.shape <- 0.01
prior.s.scale <- 0.01

## with prior
llnormfun <- function(beta, x) {
#  beta[1] -> mean
#  beta[2] -> s.d.
  if (beta[2] < 0.0) {
    l <- -Inf
  } else {
    l <- sum(log(dnorm(x, mean = beta[1], sd = beta[2]))) +
         log(dnorm(beta[1],
                   mean = prior.m.mean, sd = prior.m.sd)) +
         log(dinvgamma(beta[2]^2,
                       shape = prior.s.shape, scale = prior.s.scale))
  }
  return(l)
}

post <- MCMCmetrop1R(llnormfun, theta.init = c(5, 5), 
                     x = x, 
                     thin = 10, mcmc = 11000, burnin = 1000, 
                     tune = 0.8, 
                     verbose = 1000,
                     logfun = TRUE)

結果

> summary(post)

Iterations = 1001:11991
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1100 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean     SD Naive SE Time-series SE
[1,] 10.320 0.8113  0.02446        0.03694
[2,]  2.448 0.5995  0.01808        0.03343

2. Quantiles for each variable:

      2.5%   25%    50%    75%  97.5%
var1 8.548 9.852 10.347 10.803 11.867
var2 1.574 2.017  2.361  2.753  3.916

JAGSでも同様にやってみて、結果を比較する。

library(R2WinBUGS)
library(rjags)

model <- function() {
  for (i in 1:N) {
    x[i] ~ dnorm(m, tau);
  }
  s <- 1/sqrt(tau);
  m ~ dnorm(0, 1/(10^2));
  tau ~ dgamma(0.01, 0.01);
    ## dgamma(shape, rate)
    ## rate in dgamma is equivalant to scale in dinvgamma
}
write.model(model, "model.txt")

jags.mod <- jags.model("model.txt",
                       data = list(x = x, N = length(x)),
                       n.chains = 1,
                       n.adapt = 1000)
jags.post <- coda.samples(jags.mod,
                          variable.names = c("m", "s"),
                          n.iter = 10000,
                          thin = 10)
summary(jags.post)

結果

Iterations = 10:10000
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

    Mean     SD Naive SE Time-series SE
m 10.331 0.8291  0.02622        0.02794
s  2.516 0.7028  0.02222        0.02159

2. Quantiles for each variable:

   2.5%   25%    50%    75%  97.5%
m 8.583 9.796 10.357 10.846 12.028
s 1.601 2.040  2.382  2.801  4.337

割といい感じである。

とはいえ、これだけではまだなんとも、なのでもうちょっと試してみる。

prior.m.mean <- -20
prior.m.sd <- 1
prior.s.shape <- 10
prior.s.scale <- 1

とすると、MCMCmetrop1Rでは

Iterations = 1001:11991
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1100 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean    SD Naive SE Time-series SE
[1,] -18.96 1.006  0.03035        0.03641
[2,]  17.26 2.367  0.07135        0.08975

2. Quantiles for each variable:

       2.5%    25%    50%    75%  97.5%
var1 -20.87 -19.68 -18.95 -18.24 -16.99
var2  13.59  15.59  17.04  18.55  22.76

JAGSでは

Iterations = 10:10000
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

    Mean    SD Naive SE Time-series SE
m -19.04 1.008  0.03186        0.03349
s  17.49 2.296  0.07261        0.07239

2. Quantiles for each variable:

    2.5%    25%    50%    75%  97.5%
m -21.10 -19.71 -19.06 -18.35 -17.12
s  13.66  15.84  17.27  18.80  22.43

つづいて、

prior.m.mean <- 20
prior.m.sd <- 3
prior.s.shape <- 5
prior.s.scale <- 5

MCMCmetrop1Rでは、

Iterations = 1001:11991
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1100 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean     SD Naive SE Time-series SE
[1,] 10.739 0.5808 0.017512        0.02410
[2,]  1.819 0.3129 0.009434        0.01413

2. Quantiles for each variable:

      2.5%    25%   50%    75%  97.5%
var1 9.661 10.344 10.72 11.129 11.882
var2 1.331  1.598  1.78  1.982  2.525

JAGSでは

Iterations = 10:10000
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

    Mean    SD Naive SE Time-series SE
m 10.729 0.601 0.019006        0.01942
s  1.850 0.312 0.009866        0.01050

2. Quantiles for each variable:

   2.5%    25%    50%    75%  97.5%
m 9.601 10.302 10.698 11.142 11.967
s 1.372  1.628  1.813  2.041  2.572

まあ なんだかこれでよさそうな気がする。


と、年休消化しつつ、Rをいじる日。

タグ：MCMC R