Kosugitti 先生の世界一簡単なrstanコードに触発されて、
二番目に簡単な rstan コードを書いてみました。
library(rstan)
N <- 1000
x <- rnorm(N, mean = 50, sd = 10)
y <- 10 + 0.8 * x + rnorm(N, mean =0, sd = 7)
stancode <- '
data{
int<lower=0> N;
real x[N];
real y[N];
}
parameters {
real alpha;
real beta;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i] ~ normal(alpha + beta * x[i], s);
alpha ~ normal(0, 100);
beta ~ normal(0, 100);
s ~ inv_gamma(0.001, 0.001);
}
'
datastan <- list(N=N, x=x, y=y)
fit <- stan(model_code=stancode, data=datastan, iter=1000, chain=4)
traceplot(fit, ask=T)
print(fit, digit=2)
結果
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 の推定がうまくいってないような気がするんですが。。。
ご査収ください。