1. hoxo_m

    No comment

    hoxo_m
Changes in body
Source | HTML | Preview

Kosugitti 先生の世界一簡単な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 beta;
    real<lower=0> s;
  }
  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=2)

Rplot01.png

結果
          mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
alpha    11.60    0.06 1.11     9.29    10.90    11.60    12.33    13.73   396 1.01
beta      0.77    0.00 0.02     0.73     0.76     0.77     0.79     0.82   404 1.01
s         7.01    0.01 0.15     6.71     6.90     7.01     7.10     7.32   438 1.00
lp__  -2447.58    0.06 1.26 -2450.88 -2448.09 -2447.25 -2446.70 -2446.22   451 1.00

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

ご査収ください。