1. hoxo_m

    Posted

    hoxo_m
Changes in title
+二番目に簡単な rstan コード
Changes in tags
+R
3.0.3
+RStan
2.3.0
Changes in body
Source | HTML | Preview

世界一簡単なrstanコード
の次に簡単な rstan コードを書いてみました。

library(rstan) 

N <- 1000
y1 <- rnorm(N, mean = 50, sd = 10)
y2 <- 10 + 0.8 * y1 + rnorm(N, mean =0, sd = 7)

stancode <- '
  data{
    int<lower=0> N;
    real y1[N];
    real y2[N];
  }
  parameters {
    real alpha;
    real<lower=0> s;
    real beta;
  }
  model{
    for(i in 1:N)
      y2[i] ~ normal(alpha + beta * y1[i], s);
    s ~ inv_gamma(0.001, 0.001);
  }
'
datastan <- list(N=N, y1=y1, y2=y2)
fit <- stan(model_code=stancode, data=datastan, iter=1000, chain=4)
traceplot(fit, ask=T)
print(fit, digit=3)

Rplot.png

結果
           mean se_mean    sd      2.5%       25%       50%       75%     97.5% n_eff  Rhat
alpha     9.161   0.047 1.111     7.105     8.393     9.140     9.831    11.413   550 1.000
beta      0.819   0.001 0.022     0.774     0.806     0.820     0.834     0.861   555 1.000
s         7.024   0.006 0.152     6.735     6.917     7.023     7.125     7.327   715 1.000
lp__  -2448.119   0.064 1.224 -2451.173 -2448.605 -2447.812 -2447.263 -2446.777   360 1.009

alpha の推定がうまくいってないような気がするんですが。。。

ご査収ください。