Stan で新しい値に対する予測区間を出す #rstatsj

######################################## パッケージの読み込み
library(rstan)
library(ggplot2)
library(dplyr)

######################################## モデル作成用のデータを生成
set.seed(314)

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

ggplot(data.frame(x=x, y=y), aes(x=x, y=y)) + geom_point()

######################################## 新しい値
x_new <- seq(10, 90, by=10)

######################################## パラメータと予測値を推定する Stan コード
data{
  int<lower=0> N;
  real x[N];
  real y[N];
  int<lower=0> N_new;
  real x_new[N_new];
}
parameters {
  real alpha;
  real beta;
  real<lower=0> s;
  real y_pred[N_new]; # 予測値
}
model{
  for(i in 1:N)
    y[i] ~ normal(alpha + beta * x[i], s);
  for(i in 1:N_new)
    y_pred[i] ~ normal(alpha + beta * x_new[i], s);
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  s ~ inv_gamma(0.001, 0.001);
}

######################################## Stan により、パラメータと予測値を推定
datastan <- list(N=N, x=x, y=y, N_new=length(x_new), x_new=x_new)
fit <- stan(file="predict_for_new_value.stan", data=datastan, iter=1000, chain=4)
# traceplot(fit, ask=F)
print(fit, digit=2)

              mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
alpha         8.81    0.03 1.18     6.50     7.99     8.82     9.60    11.07  1269    1
beta          0.83    0.00 0.02     0.78     0.81     0.83     0.84     0.87  1273    1
s             7.03    0.00 0.15     6.75     6.92     7.02     7.13     7.34  2000    1
y_pred[1]    16.88    0.16 7.17     2.79    12.30    16.84    21.58    31.28  2000    1
y_pred[2]    25.43    0.16 7.12    11.40    20.48    25.41    30.54    39.10  2000    1
y_pred[3]    33.32    0.15 6.89    20.08    28.59    33.53    37.95    46.45  2000    1
y_pred[4]    41.76    0.15 6.86    28.19    37.17    41.75    46.51    54.95  2000    1
y_pred[5]    50.03    0.16 7.25    36.09    45.27    50.01    54.69    64.29  2000    1
y_pred[6]    58.50    0.16 7.26    44.25    53.70    58.45    63.19    72.69  2000    1
y_pred[7]    66.58    0.15 6.82    53.36    61.92    66.58    71.16    79.91  2000    1
y_pred[8]    74.79    0.15 6.89    61.09    70.11    74.98    79.50    87.86  2000    1
y_pred[9]    83.26    0.17 7.39    68.90    78.35    83.30    88.35    98.26  2000    1
lp__      -2471.61    0.09 2.51 -2477.60 -2472.99 -2471.31 -2469.80 -2467.78   733    1

######################################## 予測値の抽出
y_pred <- data.frame(rstan::extract(fit)$y_pred)
colnames(y_pred) <- paste0("x_", 1:9, rep(0,9))

######################################## 予測範囲を描画
data <- tidyr::gather(y_pred, x, y, seq_len(ncol(y_pred)))
ggplot(data, aes(x=x, y=y)) + geom_violin()

######################################## 予測範囲の計算
data %>% group_by(x) %>% 
  summarise(mean=mean(y), lower=quantile(y, probs=0.025), upper=quantile(y, probs=0.975))

     x     mean     lower    upper
1 x_10 16.88353  2.786534 31.28435
2 x_20 25.43116 11.404681 39.09586
3 x_30 33.31677 20.082840 46.45167
4 x_40 41.76491 28.190227 54.95213
5 x_50 50.03347 36.090850 64.29103
6 x_60 58.50089 44.254372 72.69129
7 x_70 66.58487 53.355930 79.91333
8 x_80 74.79041 61.085447 87.86033
9 x_90 83.25508 68.898657 98.25819

######################################## (参考) lm() による予測範囲の計算
predict(lm(y ~ x), newdata=data.frame(x=x_new), interval = "prediction")

       fit       lwr      upr
1 17.03153  3.131567 30.93149
2 25.29911 11.449794 39.14843
3 33.56670 19.753380 47.38002
4 41.83428 28.042209 55.62636
5 50.10187 36.316215 63.88753
6 58.36946 44.575376 72.16354
7 66.63704 52.819719 80.45437
8 74.90463 61.049319 88.75994
9 83.17222 69.264297 97.08013

# x <- rnorm(N, mean = 50, sd = 10)
x <- rnorm(N, mean = 50, sd = 1)