偉大なるごみ箱師に Stan で新しい値に対する予測区間を出す方法を教わったので、やってみる(ずいぶん放置していた)。
やり方は一言で言えば欠測値推定と同じやり方と言えばわかる人にはわかるだろう。
R
######################################## パッケージの読み込み
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()
R
######################################## 新しい値
x_new <- seq(10, 90, by=10)
predict_for_new_value.stan
######################################## パラメータと予測値を推定する 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);
}
R
######################################## 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
R
######################################## 予測値の抽出
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()
R
######################################## 予測範囲の計算
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
R
######################################## (参考) 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
ちなみに、外挿に対して予測区間が広がることは、データ生成を次のように一行変えてみて確認した。
R
# x <- rnorm(N, mean = 50, sd = 10)
x <- rnorm(N, mean = 50, sd = 1)
まとめ
この方法だと、パラメータ推定と予測値推定を同時にやってるので、同じモデルで予測を繰り返したい場合になんか無駄が多い気がするのでいい方法知っている方はご教授ください。
私がごみ箱師の説明を勘違いしている可能性もあるので、その場合もご指摘いただけるとありがたいです。