rastudylife

機械学習や統計についての勉強ログ

Stanを使ってベイズモデル化の練習

最近は,『みどりぼん』も割と読み進めてきたので『あひるぼん』を読みながら,実際にStanでMCMC実践してました.

 
↓『みどりぼん 』

 ↓『あひるぼん』

『あひるぼん』を参考に『みどりぼん』9章に出てくる種子数のベイズモデル化をStanを使って書いてみました.

  • 『みどりぼん』サポートページの第9章の場所においてある "d.RData"を使います.

生態学データ解析 - 本/データ解析のための統計モデリング入門

> green <- (load("~/d.RData"))
> MyData <- get(green)
> MyData
          x  y
1  3.000000  5
2  3.210526  3
3  3.421053  6
4  3.631579  7
5  3.842105  7
#以下略
#個体サイズxに対しての種子数y
> #glmで
> glm(y ~ x, family=poisson, data=MyData)

Call:  glm(formula = y ~ x, family = poisson, data = MyData)

Coefficients:
(Intercept)            x  
    1.56606      0.08334  

Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
Null Deviance:	    15.66 
Residual Deviance: 14.17 	AIC: 94.04

上記の結果より,平均種子数は \lambda=\exp(1.566 + 0.08x)と予測される.

data{
  int N;
  real X[N];
  int Y[N];
}

parameters{
  real b[2];
}

transformed parameters{
  real lambda[N];
  for (n in 1:N)
    lambda[n] = exp(b[1] + b[2]*X[n]);
}

model{
  for (n in 1:N)
    Y[n] ~ poisson(lambda[n]);
}

generated quantities {
  int m_pred[N];
  for (n in 1:N)
    m_pred[n] = poisson_rng(lambda[n]);
}
  • RでStanをfitすると
> library(rstan)
> data <- list(N=nrow(MyData), X=MyData$x, Y=MyData$y)
> fit <- stan(file='green_model.stan', data=data, seed=1234)
> fit
Inference for Stan model: green_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

             mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
b[1]         1.56    0.01 0.36   0.84   1.31   1.56   1.80   2.27   795    1
b[2]         0.08    0.00 0.07  -0.05   0.04   0.08   0.13   0.22   786    1
lambda[1]    6.19    0.03 1.04   4.38   5.45   6.12   6.84   8.46   923    1
# 以下略

b_1=1.56 b_2=0.08ということがわかり,glm関数で求めた数値とほぼ同じ結果が得られたことがわかった.