読者です 読者をやめる 読者になる 読者になる

廿TT

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

久保(2012)ベイズモデル入門のための例題

R Stan

概要

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

『データ解析のための統計モデリング入門』第 8 章の例題を RStan でやる.

R: RStanをためす (2):Taglibro de H:So-netブログ でも同じことがやられている.

説明は省略してコードと実行結果だけ羅列します.

glm(一般化線形モデル)

ポアソン回帰を行う. データは 生態学データ解析 - 本/データ解析のための統計モデリング入門 よりダウンロードできる.

load("~/Downloads/d.RData")
fit1 <- glm(y~x,family=poisson,data=d) #(link="identity")
plot(d)
curve(exp(coef(fit1)[1]+coef(fit1)[2]*x),add=TRUE)

f:id:abrahamcow:20150729004158p:plain
黒い実線が推定されたポアソン分布の平均.

> summary(fit1)

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

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.52168  -0.53195   0.06417   0.40797   1.57939  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.56606    0.35995   4.351 1.36e-05 ***
x            0.08334    0.06838   1.219    0.223    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 15.660  on 19  degrees of freedom
Residual deviance: 14.171  on 18  degrees of freedom
AIC: 94.036

Number of Fisher Scoring iterations: 4

RStan(ベイズモデル)

Stan のコードはこう.

data {
  int<lower=0> N;
  real<lower=0> x[N];
  int<lower=0> y[N];
}
parameters {
  real b0;
  real b1;
}
transformed parameters {
  real lambda[N];
  for (i in 1:N) {
      lambda[i] <- exp(b0 + b1*x[i]); #リンク関数
      }
}
model {
  for (i in 1:N) {
      y[i] ~ poisson( lambda[i] );
      }
      b0 ~ normal(0, 100); #無情報事前分布
      b1 ~ normal(0, 100);
}

これを R から動かす.

library(rstan)
N <- dim(d)[1]
dat <-list(N=N,x=d$x,y=d$y)
fit <- stan(file="poissonreg.stan",model_name="pois",par=c("b0","b1"), data=dat)
traceplot(fit)

f:id:abrahamcow:20150729005109p:plain

> ex <-extract(fit)
> summary(ex$b0)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.351   1.337   1.571   1.566   1.821   2.667 
> summary(ex$b1)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.12100  0.03569  0.08158  0.08208  0.12730  0.30800 

おおよそ glm とおなじような値が求まった.

ベイズモデルでは母数を確率変数だと考える. 係数 b0, b1 の分布は下図の通り.

f:id:abrahamcow:20150729005337p:plain

f:id:abrahamcow:20150729005344p:plain

MCMCpack(おまけ)

func <- function(beta) {
  link  <- exp(beta[1] + beta[2] * d$x)
  sum(dpois(d$y,link,log=TRUE))
}
optim(par=c(1,1),fn=func,control=list(fnscale=-1))
library(MCMCpack)
d.res <- MCMCmetrop1R(func, theta.init = c(0, 0))
plot(d.res)

f:id:abrahamcow:20150729005936p:plain

> summary(d.res)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

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

        Mean      SD  Naive SE Time-series SE
[1,] 1.56271 0.35362 0.0025004       0.007641
[2,] 0.08313 0.06696 0.0004735       0.001445

2. Quantiles for each variable:

         2.5%     25%     50%    75% 97.5%
var1  0.83922 1.32887 1.56691 1.8066 2.233
var2 -0.04749 0.03717 0.08335 0.1277 0.217

関連エントリ

abrahamcow.hatenablog.com