StanとRでベイズ統計モデリング04
UPDATE: 2023-12-16 11:39:26.455421
はじめに
このノートは「StanとRでベイズ統計モデリング」の内容を写経することで、ベイズ統計への理解を深めていくために作成している。
基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、「StanとRでベイズ統計モデリング」を読み進めるための自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
今回は第7章「回帰分析の悩みどころ」のチャプターを写経していく。
7.2 対数を取るか否か
対数を取らない通常の回帰の例は書籍を見るとして、ここでは対数を取る場合のメモを残しておく。対数変換を取る場合のモデル式は下記の通り。
\(log_{10} (Area[n])\)から\(\mu\)が決まり、\(\mu\)にノイズが乗って\(log_{10} (Y[n])\)が生成される。
\[ \begin{eqnarray} \mu[n] &=& b_{1} + b_{2} log_{10} (Area[n]) \\ log_{10} (Y[n]) &\sim& Normal(\mu[n], \sigma) \end{eqnarray} \]
対数を取るので、Stanに渡すデータを対数変換(Area = log10(d$X), Y = log10(d$Y), Area_new = log10(X_new))してから渡している。こちらで対数変換せずに渡すのであれば、Stanのモデル式で変換を施せば良いと思われる。
options(max.print = 999999)
library(dplyr)
library(ggplot2)
library(rstan)
d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap07/input/data-rental.txt')
colnames(d) <- c('Y', 'X')
X_new <- seq(from = 10, to = 120, length = 50)
data <- list(N = nrow(d), Area = log10(d$X), Y = log10(d$Y), N_new = 50, Area_new = log10(X_new))
data## $N
## [1] 100
##
## $Area
## [1] 1.609167 1.791971 1.741388 1.754501 1.754501 1.190332 1.505150 1.521269
## [9] 1.805025 1.301030 1.301030 1.918083 1.423082 1.373464 1.351796 1.297323
## [17] 1.364363 1.302980 1.370513 1.618884 1.643453 1.290035 1.632356 1.725422
## [25] 1.488974 1.808481 1.368659 1.600319 1.608526 1.483587 1.491362 1.733839
## [33] 2.052386 1.322219 1.192289 1.322219 1.322219 1.569725 1.363236 1.756636
## [41] 1.524136 1.396548 1.065953 1.499137 1.550717 1.462398 1.860697 1.653309
## [49] 1.385606 1.322219 1.563006 1.737829 1.542701 1.462398 1.505150 1.269513
## [57] 1.669317 1.687975 1.415974 1.464788 1.255273 1.230449 1.703291 1.498035
## [65] 1.260071 1.431685 1.276462 1.710117 1.484585 1.662380 1.227630 1.361917
## [73] 1.480582 1.329398 1.298198 1.653213 1.910944 1.582631 1.646208 1.220370
## [81] 1.834039 1.691435 1.577607 1.356408 1.421933 2.041630 1.545307 1.715586
## [89] 1.371437 1.582291 1.506505 1.518514 1.704837 1.496376 1.397940 1.287802
## [97] 1.448242 1.776047 1.516932 1.614897
##
## $Y
## [1] 2.420058 2.715167 2.826981 2.674889 2.658584 2.186108 2.369216 2.422393
## [9] 2.812913 2.348694 2.206286 2.878522 2.344196 2.304491 2.390935 2.207365
## [17] 2.301464 2.292699 2.260787 2.695919 2.464454 2.198657 2.675063 2.802363
## [25] 2.492201 2.938520 2.486615 2.637690 2.434569 2.414137 2.192010 2.545307
## [33] 3.271958 2.265290 2.253822 2.265290 2.317227 2.561698 2.371068 2.635484
## [41] 2.263873 2.473633 2.095518 2.450588 2.493597 2.392697 3.091597 2.622421
## [49] 2.268363 2.265290 2.541330 2.946845 2.656098 2.265996 2.369216 2.281601
## [57] 2.494850 2.777064 2.502973 2.295585 2.162564 2.169674 2.773494 2.616370
## [65] 2.246006 2.398981 2.287802 2.603253 2.358316 2.678354 2.165541 2.215648
## [73] 2.322374 2.315970 2.252853 2.547036 3.022428 2.414973 2.591065 2.189939
## [81] 2.803662 2.662002 2.892762 2.254790 2.381296 3.244153 2.328991 2.702172
## [89] 2.289589 2.337060 2.421604 2.344392 2.448397 2.590730 2.293031 2.321201
## [97] 2.344392 2.781396 2.511883 2.626340
##
## $N_new
## [1] 50
##
## $Area_new
## [1] 1.000000 1.087955 1.161062 1.223618 1.278287 1.326837 1.370502 1.410174
## [9] 1.446524 1.480066 1.511201 1.540253 1.567482 1.593105 1.617300 1.640218
## [17] 1.661986 1.682716 1.702501 1.721424 1.739556 1.756962 1.773697 1.789811
## [25] 1.805348 1.820349 1.834849 1.848880 1.862472 1.875652 1.888443 1.900869
## [33] 1.912948 1.924701 1.936144 1.947294 1.958164 1.968769 1.979121 1.989232
## [41] 1.999113 2.008774 2.018225 2.027474 2.036531 2.045403 2.054097 2.062620
## [49] 2.070980 2.079181モデルは下記の通り。
data {
int N;
real Area[N];
real Y[N];
int N_new;
real Area_new[N_new];
}
parameters {
real b1;
real b2;
real<lower=0> sigma;
}
transformed parameters {
real mu[N];
for (n in 1:N){
mu[n] = b1 + b2*Area[n];
}
}
model {
for (n in 1:N)
Y[n] ~ normal(mu[n], sigma);
}
generated quantities {
real y_pred[N];
real y_new[N_new];
for (n in 1:N){
y_pred[n] = normal_rng(mu[n], sigma);
}
for (n in 1:N_new){
y_new[n] = normal_rng(b1 + b2*Area_new[n], sigma);
}
}ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果はこちら。y_pred[n]が1つの物件ごとの物件価格の予測分布で、y_new[n]が各エリアサイズごとの予測区間を計算したもの。
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## b1 0.8 0 0.1 0.6 0.8 1.0 1168 1
## b2 1.1 0 0.1 1.0 1.1 1.2 1168 1
## sigma 0.1 0 0.0 0.1 0.1 0.1 1676 1
## mu[1] 2.6 0 0.0 2.6 2.6 2.6 3166 1
## mu[2] 2.8 0 0.0 2.7 2.8 2.8 1642 1
## mu[3] 2.7 0 0.0 2.7 2.7 2.8 1830 1
## mu[4] 2.7 0 0.0 2.7 2.7 2.8 1772 1
## mu[5] 2.7 0 0.0 2.7 2.7 2.8 1772 1
## mu[6] 2.1 0 0.0 2.1 2.1 2.2 1333 1
## mu[7] 2.5 0 0.0 2.4 2.5 2.5 4586 1
## mu[8] 2.5 0 0.0 2.5 2.5 2.5 4657 1
## mu[9] 2.8 0 0.0 2.8 2.8 2.8 1607 1
## mu[10] 2.2 0 0.0 2.2 2.2 2.3 1510 1
## mu[11] 2.2 0 0.0 2.2 2.2 2.3 1510 1
## mu[12] 2.9 0 0.0 2.9 2.9 3.0 1418 1
## mu[13] 2.4 0 0.0 2.4 2.4 2.4 2534 1
## mu[14] 2.3 0 0.0 2.3 2.3 2.3 1914 1
## mu[15] 2.3 0 0.0 2.3 2.3 2.3 1753 1
## mu[16] 2.2 0 0.0 2.2 2.2 2.3 1497 1
## mu[17] 2.3 0 0.0 2.3 2.3 2.3 1841 1
## mu[18] 2.2 0 0.0 2.2 2.2 2.3 1516 1
## mu[19] 2.3 0 0.0 2.3 2.3 2.3 1889 1
## mu[20] 2.6 0 0.0 2.6 2.6 2.6 2990 1
## mu[21] 2.6 0 0.0 2.6 2.6 2.6 2610 1
## mu[22] 2.2 0 0.0 2.2 2.2 2.3 1473 1
## mu[23] 2.6 0 0.0 2.6 2.6 2.6 2770 1
## mu[24] 2.7 0 0.0 2.7 2.7 2.7 1913 1
## mu[25] 2.4 0 0.0 2.4 2.4 2.5 4428 1
## mu[26] 2.8 0 0.0 2.8 2.8 2.8 1598 1
## mu[27] 2.3 0 0.0 2.3 2.3 2.3 1874 1
## mu[28] 2.6 0 0.0 2.5 2.6 2.6 3355 1
## mu[29] 2.6 0 0.0 2.6 2.6 2.6 3178 1
## mu[30] 2.4 0 0.0 2.4 2.4 2.5 4360 1
## mu[31] 2.5 0 0.0 2.4 2.5 2.5 4456 1
## mu[32] 2.7 0 0.0 2.7 2.7 2.8 1868 1
## mu[33] 3.1 0 0.0 3.0 3.1 3.1 1321 1
## mu[34] 2.3 0 0.0 2.2 2.3 2.3 1594 1
## mu[35] 2.1 0 0.0 2.1 2.1 2.2 1335 1
## mu[36] 2.3 0 0.0 2.2 2.3 2.3 1594 1
## mu[37] 2.3 0 0.0 2.2 2.3 2.3 1594 1
## mu[38] 2.5 0 0.0 2.5 2.5 2.6 4324 1
## mu[39] 2.3 0 0.0 2.3 2.3 2.3 1833 1
## mu[40] 2.7 0 0.0 2.7 2.7 2.8 1763 1
## mu[41] 2.5 0 0.0 2.5 2.5 2.5 4660 1
## mu[42] 2.3 0 0.0 2.3 2.3 2.4 2144 1
## mu[43] 2.0 0 0.0 1.9 2.0 2.0 1268 1
## mu[44] 2.5 0 0.0 2.4 2.5 2.5 4536 1
## mu[45] 2.5 0 0.0 2.5 2.5 2.5 4537 1
## mu[46] 2.4 0 0.0 2.4 2.4 2.4 3535 1
## mu[47] 2.9 0 0.0 2.8 2.9 2.9 1493 1
## mu[48] 2.6 0 0.0 2.6 2.6 2.7 2472 1
## mu[49] 2.3 0 0.0 2.3 2.3 2.4 2026 1
## mu[50] 2.3 0 0.0 2.2 2.3 2.3 1594 1
## mu[51] 2.5 0 0.0 2.5 2.5 2.6 4408 1
## mu[52] 2.7 0 0.0 2.7 2.7 2.8 1848 1
## mu[53] 2.5 0 0.0 2.5 2.5 2.5 4599 1
## mu[54] 2.4 0 0.0 2.4 2.4 2.4 3535 1
## mu[55] 2.5 0 0.0 2.4 2.5 2.5 4586 1
## mu[56] 2.2 0 0.0 2.2 2.2 2.2 1431 1
## mu[57] 2.6 0 0.0 2.6 2.6 2.7 2293 1
## mu[58] 2.7 0 0.0 2.6 2.7 2.7 2163 1
## mu[59] 2.4 0 0.0 2.3 2.4 2.4 2408 1
## mu[60] 2.4 0 0.0 2.4 2.4 2.4 3648 1
## mu[61] 2.2 0 0.0 2.2 2.2 2.2 1407 1
## mu[62] 2.2 0 0.0 2.1 2.2 2.2 1373 1
## mu[63] 2.7 0 0.0 2.7 2.7 2.7 2055 1
## mu[64] 2.5 0 0.0 2.4 2.5 2.5 4526 1
## mu[65] 2.2 0 0.0 2.2 2.2 2.2 1415 1
## mu[66] 2.4 0 0.0 2.4 2.4 2.4 2704 1
## mu[67] 2.2 0 0.0 2.2 2.2 2.2 1444 1
## mu[68] 2.7 0 0.0 2.7 2.7 2.7 2008 1
## mu[69] 2.4 0 0.0 2.4 2.4 2.5 4373 1
## mu[70] 2.6 0 0.0 2.6 2.6 2.7 2365 1
## mu[71] 2.2 0 0.0 2.1 2.2 2.2 1370 1
## mu[72] 2.3 0 0.0 2.3 2.3 2.3 1823 1
## mu[73] 2.4 0 0.0 2.4 2.4 2.5 4319 1
## mu[74] 2.3 0 0.0 2.2 2.3 2.3 1627 1
## mu[75] 2.2 0 0.0 2.2 2.2 2.3 1500 1
## mu[76] 2.6 0 0.0 2.6 2.6 2.7 2473 1
## mu[77] 2.9 0 0.0 2.9 2.9 3.0 1426 1
## mu[78] 2.6 0 0.0 2.5 2.6 2.6 4077 1
## mu[79] 2.6 0 0.0 2.6 2.6 2.6 2571 1
## mu[80] 2.2 0 0.0 2.1 2.2 2.2 1361 1
## mu[81] 2.8 0 0.0 2.8 2.8 2.9 1541 1
## mu[82] 2.7 0 0.0 2.6 2.7 2.7 2141 1
## mu[83] 2.5 0 0.0 2.5 2.5 2.6 4218 1
## mu[84] 2.3 0 0.0 2.3 2.3 2.3 1784 1
## mu[85] 2.4 0 0.0 2.4 2.4 2.4 2513 1
## mu[86] 3.1 0 0.0 3.0 3.1 3.1 1326 1
## mu[87] 2.5 0 0.0 2.5 2.5 2.5 4581 1
## mu[88] 2.7 0 0.0 2.7 2.7 2.7 1972 1
## mu[89] 2.3 0 0.0 2.3 2.3 2.3 1897 1
## mu[90] 2.6 0 0.0 2.5 2.6 2.6 4094 1
## mu[91] 2.5 0 0.0 2.4 2.5 2.5 4596 1
## mu[92] 2.5 0 0.0 2.5 2.5 2.5 4652 1
## mu[93] 2.7 0 0.0 2.7 2.7 2.7 2044 1
## mu[94] 2.5 0 0.0 2.4 2.5 2.5 4510 1
## mu[95] 2.3 0 0.0 2.3 2.3 2.4 2160 1
## mu[96] 2.2 0 0.0 2.2 2.2 2.3 1467 1
## mu[97] 2.4 0 0.0 2.4 2.4 2.4 3088 1
## mu[98] 2.8 0 0.0 2.7 2.8 2.8 1692 1
## mu[99] 2.5 0 0.0 2.5 2.5 2.5 4648 1
## mu[100] 2.6 0 0.0 2.6 2.6 2.6 3060 1
## y_pred[1] 2.6 0 0.1 2.4 2.6 2.8 3897 1
## y_pred[2] 2.8 0 0.1 2.6 2.8 3.0 3548 1
## y_pred[3] 2.7 0 0.1 2.5 2.7 3.0 3668 1
## y_pred[4] 2.7 0 0.1 2.5 2.7 3.0 3726 1
## y_pred[5] 2.7 0 0.1 2.5 2.7 3.0 3734 1
## y_pred[6] 2.1 0 0.1 1.9 2.1 2.3 3428 1
## y_pred[7] 2.5 0 0.1 2.2 2.5 2.7 4038 1
## y_pred[8] 2.5 0 0.1 2.3 2.5 2.7 4236 1
## y_pred[9] 2.8 0 0.1 2.6 2.8 3.0 3627 1
## y_pred[10] 2.2 0 0.1 2.0 2.2 2.5 4119 1
## y_pred[11] 2.2 0 0.1 2.0 2.2 2.5 4140 1
## y_pred[12] 2.9 0 0.1 2.7 2.9 3.2 3868 1
## y_pred[13] 2.4 0 0.1 2.2 2.4 2.6 3939 1
## y_pred[14] 2.3 0 0.1 2.1 2.3 2.5 3654 1
## y_pred[15] 2.3 0 0.1 2.1 2.3 2.5 3841 1
## y_pred[16] 2.2 0 0.1 2.0 2.2 2.5 3634 1
## y_pred[17] 2.3 0 0.1 2.1 2.3 2.5 4073 1
## y_pred[18] 2.2 0 0.1 2.0 2.2 2.5 4215 1
## y_pred[19] 2.3 0 0.1 2.1 2.3 2.5 4050 1
## y_pred[20] 2.6 0 0.1 2.4 2.6 2.8 4218 1
## y_pred[21] 2.6 0 0.1 2.4 2.6 2.8 4117 1
## y_pred[22] 2.2 0 0.1 2.0 2.2 2.5 3965 1
## y_pred[23] 2.6 0 0.1 2.4 2.6 2.8 3879 1
## y_pred[24] 2.7 0 0.1 2.5 2.7 2.9 3849 1
## y_pred[25] 2.5 0 0.1 2.2 2.5 2.7 4116 1
## y_pred[26] 2.8 0 0.1 2.6 2.8 3.0 3956 1
## y_pred[27] 2.3 0 0.1 2.1 2.3 2.5 3823 1
## y_pred[28] 2.6 0 0.1 2.3 2.6 2.8 4061 1
## y_pred[29] 2.6 0 0.1 2.4 2.6 2.8 4002 1
## y_pred[30] 2.4 0 0.1 2.2 2.4 2.7 3835 1
## y_pred[31] 2.4 0 0.1 2.2 2.4 2.7 3772 1
## y_pred[32] 2.7 0 0.1 2.5 2.7 2.9 3828 1
## y_pred[33] 3.1 0 0.1 2.8 3.1 3.3 3449 1
## y_pred[34] 2.3 0 0.1 2.0 2.3 2.5 3926 1
## y_pred[35] 2.1 0 0.1 1.9 2.1 2.3 3534 1
## y_pred[36] 2.3 0 0.1 2.0 2.3 2.5 3996 1
## y_pred[37] 2.3 0 0.1 2.1 2.3 2.5 3853 1
## y_pred[38] 2.5 0 0.1 2.3 2.5 2.8 4163 1
## y_pred[39] 2.3 0 0.1 2.1 2.3 2.5 3925 1
## y_pred[40] 2.7 0 0.1 2.5 2.7 3.0 3882 1
## y_pred[41] 2.5 0 0.1 2.3 2.5 2.7 3689 1
## y_pred[42] 2.3 0 0.1 2.1 2.3 2.6 3812 1
## y_pred[43] 2.0 0 0.1 1.8 2.0 2.2 3242 1
## y_pred[44] 2.5 0 0.1 2.2 2.5 2.7 3562 1
## y_pred[45] 2.5 0 0.1 2.3 2.5 2.7 3925 1
## y_pred[46] 2.4 0 0.1 2.2 2.4 2.6 3922 1
## y_pred[47] 2.9 0 0.1 2.6 2.9 3.1 3892 1
## y_pred[48] 2.6 0 0.1 2.4 2.6 2.9 3977 1
## y_pred[49] 2.3 0 0.1 2.1 2.3 2.6 3974 1
## y_pred[50] 2.3 0 0.1 2.0 2.3 2.5 3952 1
## y_pred[51] 2.5 0 0.1 2.3 2.5 2.7 4151 1
## y_pred[52] 2.7 0 0.1 2.5 2.7 2.9 3612 1
## y_pred[53] 2.5 0 0.1 2.3 2.5 2.7 3942 1
## y_pred[54] 2.4 0 0.1 2.2 2.4 2.6 3942 1
## y_pred[55] 2.5 0 0.1 2.2 2.5 2.7 4305 1
## y_pred[56] 2.2 0 0.1 2.0 2.2 2.4 3880 1
## y_pred[57] 2.6 0 0.1 2.4 2.6 2.9 4039 1
## y_pred[58] 2.7 0 0.1 2.5 2.7 2.9 3598 1
## y_pred[59] 2.4 0 0.1 2.1 2.4 2.6 3699 1
## y_pred[60] 2.4 0 0.1 2.2 2.4 2.6 3754 1
## y_pred[61] 2.2 0 0.1 2.0 2.2 2.4 3561 1
## y_pred[62] 2.2 0 0.1 1.9 2.2 2.4 3697 1
## y_pred[63] 2.7 0 0.1 2.5 2.7 2.9 3859 1
## y_pred[64] 2.5 0 0.1 2.2 2.5 2.7 4072 1
## y_pred[65] 2.2 0 0.1 2.0 2.2 2.4 3373 1
## y_pred[66] 2.4 0 0.1 2.2 2.4 2.6 3919 1
## y_pred[67] 2.2 0 0.1 2.0 2.2 2.4 3384 1
## y_pred[68] 2.7 0 0.1 2.5 2.7 2.9 3838 1
## y_pred[69] 2.4 0 0.1 2.2 2.4 2.7 3933 1
## y_pred[70] 2.6 0 0.1 2.4 2.6 2.9 3940 1
## y_pred[71] 2.2 0 0.1 1.9 2.2 2.4 3720 1
## y_pred[72] 2.3 0 0.1 2.1 2.3 2.5 3969 1
## y_pred[73] 2.4 0 0.1 2.2 2.4 2.7 3769 1
## y_pred[74] 2.3 0 0.1 2.1 2.3 2.5 3825 1
## y_pred[75] 2.2 0 0.1 2.0 2.2 2.5 4010 1
## y_pred[76] 2.6 0 0.1 2.4 2.6 2.9 3772 1
## y_pred[77] 2.9 0 0.1 2.7 2.9 3.1 3805 1
## y_pred[78] 2.6 0 0.1 2.3 2.6 2.8 4058 1
## y_pred[79] 2.6 0 0.1 2.4 2.6 2.8 3904 1
## y_pred[80] 2.2 0 0.1 1.9 2.1 2.4 3585 1
## y_pred[81] 2.8 0 0.1 2.6 2.8 3.1 3645 1
## y_pred[82] 2.7 0 0.1 2.4 2.7 2.9 3953 1
## y_pred[83] 2.5 0 0.1 2.3 2.5 2.8 4124 1
## y_pred[84] 2.3 0 0.1 2.1 2.3 2.5 3678 1
## y_pred[85] 2.4 0 0.1 2.2 2.4 2.6 4120 1
## y_pred[86] 3.1 0 0.1 2.8 3.1 3.3 3126 1
## y_pred[87] 2.5 0 0.1 2.3 2.5 2.7 3972 1
## y_pred[88] 2.7 0 0.1 2.5 2.7 2.9 3761 1
## y_pred[89] 2.3 0 0.1 2.1 2.3 2.5 3943 1
## y_pred[90] 2.6 0 0.1 2.3 2.6 2.8 3228 1
## y_pred[91] 2.5 0 0.1 2.2 2.5 2.7 4091 1
## y_pred[92] 2.5 0 0.1 2.3 2.5 2.7 3987 1
## y_pred[93] 2.7 0 0.1 2.5 2.7 2.9 3788 1
## y_pred[94] 2.5 0 0.1 2.2 2.5 2.7 4060 1
## y_pred[95] 2.3 0 0.1 2.1 2.3 2.6 4191 1
## y_pred[96] 2.2 0 0.1 2.0 2.2 2.5 3831 1
## y_pred[97] 2.4 0 0.1 2.2 2.4 2.6 3707 1
## y_pred[98] 2.8 0 0.1 2.5 2.8 3.0 4090 1
## y_pred[99] 2.5 0 0.1 2.3 2.5 2.7 3825 1
## y_pred[100] 2.6 0 0.1 2.4 2.6 2.8 3976 1
## y_new[1] 1.9 0 0.1 1.7 1.9 2.1 3495 1
## y_new[2] 2.0 0 0.1 1.8 2.0 2.2 3606 1
## y_new[3] 2.1 0 0.1 1.9 2.1 2.3 3797 1
## y_new[4] 2.2 0 0.1 1.9 2.2 2.4 4069 1
## y_new[5] 2.2 0 0.1 2.0 2.2 2.4 3948 1
## y_new[6] 2.3 0 0.1 2.0 2.3 2.5 3856 1
## y_new[7] 2.3 0 0.1 2.1 2.3 2.5 3855 1
## y_new[8] 2.4 0 0.1 2.1 2.4 2.6 3857 1
## y_new[9] 2.4 0 0.1 2.2 2.4 2.6 4067 1
## y_new[10] 2.4 0 0.1 2.2 2.4 2.7 3930 1
## y_new[11] 2.5 0 0.1 2.2 2.5 2.7 4112 1
## y_new[12] 2.5 0 0.1 2.3 2.5 2.7 3812 1
## y_new[13] 2.5 0 0.1 2.3 2.5 2.8 4171 1
## y_new[14] 2.6 0 0.1 2.3 2.6 2.8 4085 1
## y_new[15] 2.6 0 0.1 2.4 2.6 2.8 3831 1
## y_new[16] 2.6 0 0.1 2.4 2.6 2.8 3552 1
## y_new[17] 2.6 0 0.1 2.4 2.6 2.9 4031 1
## y_new[18] 2.7 0 0.1 2.4 2.7 2.9 3897 1
## y_new[19] 2.7 0 0.1 2.5 2.7 2.9 3682 1
## y_new[20] 2.7 0 0.1 2.5 2.7 2.9 3742 1
## y_new[21] 2.7 0 0.1 2.5 2.7 2.9 3862 1
## y_new[22] 2.7 0 0.1 2.5 2.7 3.0 3686 1
## y_new[23] 2.8 0 0.1 2.5 2.8 3.0 3910 1
## y_new[24] 2.8 0 0.1 2.6 2.8 3.0 3538 1
## y_new[25] 2.8 0 0.1 2.6 2.8 3.0 4123 1
## y_new[26] 2.8 0 0.1 2.6 2.8 3.0 3774 1
## y_new[27] 2.8 0 0.1 2.6 2.8 3.1 3977 1
## y_new[28] 2.8 0 0.1 2.6 2.8 3.1 3302 1
## y_new[29] 2.9 0 0.1 2.6 2.9 3.1 3871 1
## y_new[30] 2.9 0 0.1 2.7 2.9 3.1 3468 1
## y_new[31] 2.9 0 0.1 2.7 2.9 3.1 3671 1
## y_new[32] 2.9 0 0.1 2.7 2.9 3.1 3808 1
## y_new[33] 2.9 0 0.1 2.7 2.9 3.1 3726 1
## y_new[34] 2.9 0 0.1 2.7 2.9 3.2 3586 1
## y_new[35] 2.9 0 0.1 2.7 2.9 3.2 3668 1
## y_new[36] 3.0 0 0.1 2.7 3.0 3.2 3980 1
## y_new[37] 3.0 0 0.1 2.7 3.0 3.2 3222 1
## y_new[38] 3.0 0 0.1 2.8 3.0 3.2 3704 1
## y_new[39] 3.0 0 0.1 2.8 3.0 3.2 3658 1
## y_new[40] 3.0 0 0.1 2.8 3.0 3.2 3256 1
## y_new[41] 3.0 0 0.1 2.8 3.0 3.2 3758 1
## y_new[42] 3.0 0 0.1 2.8 3.0 3.3 3068 1
## y_new[43] 3.0 0 0.1 2.8 3.0 3.3 3491 1
## y_new[44] 3.0 0 0.1 2.8 3.0 3.3 3476 1
## y_new[45] 3.1 0 0.1 2.8 3.1 3.3 3135 1
## y_new[46] 3.1 0 0.1 2.8 3.1 3.3 3692 1
## y_new[47] 3.1 0 0.1 2.8 3.1 3.3 3435 1
## y_new[48] 3.1 0 0.1 2.9 3.1 3.3 3362 1
## y_new[49] 3.1 0 0.1 2.9 3.1 3.3 3779 1
## y_new[50] 3.1 0 0.1 2.9 3.1 3.3 2822 1
## lp__ 167.6 0 1.3 164.1 168.0 169.1 1140 1
##
## Samples were drawn using NUTS(diag_e) at Sat Dec 16 11:39:50 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).まずは各エリアサイズごとの予測区間を可視化する。ms$y_newはMCMCサンプル4000行、エリアサイズの範囲(10-120)を50分割した50列が保存されている。ms$y_newの値は対数のままなので、元に戻すために10^ms$y_newとしている。
ms <- rstan::extract(fit)
qua <- apply(10^ms$y_new, 2, quantile, probs = c(0.1, 0.25, 0.50, 0.75, 0.9))
d_est <- data.frame(X = X_new, t(qua), check.names = FALSE)
ggplot() +
theme_bw(base_size = 18) +
geom_ribbon(data = d_est, aes(x = X, ymin = `10%`, ymax = `90%`), fill = 'black', alpha = 1/6) +
geom_ribbon(data = d_est, aes(x = X, ymin = `25%`, ymax = `75%`), fill = 'black', alpha = 2/6) +
geom_line(data = d_est, aes(x = X, y = `50%`), size = 1) +
geom_point(data = d, aes(x = X, y = Y), shape = 1, size = 2) +
coord_cartesian(xlim = c(11, 118), ylim = c(-50, 1900)) +
labs(x = 'Area', y = 'Y')
次は観測値と予測分布の関係に関する散布図。ms$y_predの値は対数のままなので、元に戻すために10^ms$y_predとしている。
qua <- apply(10^ms$y_pred, 2, quantile, probs = c(0.1, 0.25, 0.50, 0.75, 0.9))
d_est <- data.frame(X = d$Y, t(qua), check.names = FALSE)
ggplot(data = d_est, aes(x = X, y = `50%`)) +
theme_bw(base_size = 18) +
coord_fixed(ratio = 1, xlim = c(-50, 1900), ylim = c(-50, 1900)) +
geom_pointrange(aes(ymin = `10%`, ymax = `90%`), col = 'grey5', fill = 'grey95', shape = 21) +
geom_abline(aes(slope = 1, intercept = 0), col = 'black', alpha = 3/5, linetype = 'dashed') +
labs(x = 'Observed', y = 'Predicted')
これら2つの可視化から、物件価格がマイナスを取らないようになったことやモデルの当てはまりが改善されていることがわかる。ただ、書籍に書かれている通り、逆にエリアサイズが大きい部分での予測は悪くなっている。
7.3 非線形な関係
Stanは自由にモデリングできるのも特徴の一つ。下記の通り、非線形なモデリングも計算可能である。このモデルでは、1時点ごとの\(Y[t]\)は平均\(a \{1 - exp(-b \ Time[t]) \}\)、標準偏差\(\sigma\)から生成されると仮定している(はず)。(はず)と書いたのは、書籍では微分方程式の解として出てくるモデルであり、経過時間と血中濃度の関係を例として書かれている。経過時点ごとに血中濃度は独立として良いのだろうかという疑問があったため・・・このモデルは独立に生成されていると解釈するほうが誤りなのだろうか。そのあたりのよくわかってない。
\[ \begin{eqnarray} Y[t] &\sim& Normal(a \{1 - exp(-b \ Time[t]) \}, \sigma) \end{eqnarray} \]
ここで\(a\)は\(Y\)の上限を決めるパラメタであり、\(b\)は\(Y\)の曲線の立ち上がりの強さを決めるパラメタ。増減表を書けばそうなっている(矢印は上に凸な矢印)。
\[ \begin{array}{c|ccccc} x & 0 & \cdots & + \infty \\ \hline f’(x) & ab & \cdots & 0 \\ \hline f’’(x) & -ab^{2} & \cdots & 0 \\ \hline f(x) & 0 & \curvearrowright & a \end{array} \]
時点ごとの予測区間を算出するため、Time_newを定義している。
d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap07/input/data-conc.txt')
T_new <- 60
Time_new <- seq(from = 0, to = 24, length = T_new)
data <- list(T = nrow(d), Time = d$Time, Y = d$Y, T_new = T_new, Time_new = Time_new)
data## $T
## [1] 6
##
## $Time
## [1] 1 2 4 8 12 24
##
## $Y
## [1] 2.4 5.0 7.5 11.9 12.5 12.7
##
## $T_new
## [1] 60
##
## $Time_new
## [1] 0.0000000 0.4067797 0.8135593 1.2203390 1.6271186 2.0338983
## [7] 2.4406780 2.8474576 3.2542373 3.6610169 4.0677966 4.4745763
## [13] 4.8813559 5.2881356 5.6949153 6.1016949 6.5084746 6.9152542
## [19] 7.3220339 7.7288136 8.1355932 8.5423729 8.9491525 9.3559322
## [25] 9.7627119 10.1694915 10.5762712 10.9830508 11.3898305 11.7966102
## [31] 12.2033898 12.6101695 13.0169492 13.4237288 13.8305085 14.2372881
## [37] 14.6440678 15.0508475 15.4576271 15.8644068 16.2711864 16.6779661
## [43] 17.0847458 17.4915254 17.8983051 18.3050847 18.7118644 19.1186441
## [49] 19.5254237 19.9322034 20.3389831 20.7457627 21.1525424 21.5593220
## [55] 21.9661017 22.3728814 22.7796610 23.1864407 23.5932203 24.0000000モデル式は下記の通り。
data {
int T;
real Time[T];
real Y[T];
int T_new;
real Time_new[T_new];
}
parameters {
real<lower=0, upper=100> a; // 制限がないと収束しない
real<lower=0, upper=5> b; // 制限がないと収束しない
real<lower=0> s_Y;
}
model {
for (t in 1:T){
Y[t] ~ normal(a*(1 - exp(-b*Time[t])), s_Y);
}
}
generated quantities {
real y_new[T_new];
for (t in 1:T_new){
y_new[t] = normal_rng(a*(1 - exp(-b*Time_new[t])), s_Y);
}
}ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果はこちら。y_new[n]が1時点ごとの血中濃度の予測分布。
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## a 13.1 0.0 1.0 11.3 13.1 14.7 806 1
## b 0.3 0.0 0.1 0.2 0.2 0.4 423 1
## s_Y 0.9 0.0 0.7 0.3 0.7 2.7 472 1
## y_new[1] 0.0 0.0 1.1 -2.2 0.0 2.0 3943 1
## y_new[2] 1.2 0.0 1.2 -0.8 1.2 3.3 2173 1
## y_new[3] 2.3 0.0 1.3 0.1 2.3 4.6 2577 1
## y_new[4] 3.4 0.0 1.2 1.4 3.3 5.9 2199 1
## y_new[5] 4.2 0.0 1.3 1.9 4.2 6.7 2884 1
## y_new[6] 5.1 0.0 1.3 2.9 5.0 7.6 2622 1
## y_new[7] 5.8 0.0 1.3 3.5 5.8 8.4 1771 1
## y_new[8] 6.5 0.0 1.3 4.0 6.4 9.1 2093 1
## y_new[9] 7.1 0.0 1.3 4.8 7.0 9.7 2727 1
## y_new[10] 7.6 0.0 1.3 5.1 7.6 10.1 2868 1
## y_new[11] 8.1 0.0 1.2 5.8 8.1 10.6 2420 1
## y_new[12] 8.6 0.0 1.2 6.2 8.6 11.1 3809 1
## y_new[13] 9.0 0.0 1.3 6.6 9.0 11.6 3334 1
## y_new[14] 9.4 0.0 1.3 6.9 9.3 11.7 3753 1
## y_new[15] 9.7 0.0 1.3 7.3 9.7 12.1 3684 1
## y_new[16] 10.0 0.0 1.3 7.5 10.0 12.5 4026 1
## y_new[17] 10.3 0.0 1.3 7.7 10.3 12.6 3726 1
## y_new[18] 10.5 0.0 1.3 8.0 10.5 12.9 3719 1
## y_new[19] 10.7 0.0 1.2 8.3 10.8 13.0 3256 1
## y_new[20] 11.0 0.0 1.3 8.6 11.0 13.4 3913 1
## y_new[21] 11.1 0.0 1.2 8.7 11.1 13.4 2913 1
## y_new[22] 11.3 0.0 1.3 9.0 11.3 13.6 3418 1
## y_new[23] 11.4 0.0 1.2 9.1 11.5 13.6 2596 1
## y_new[24] 11.6 0.0 1.4 9.0 11.6 13.8 3406 1
## y_new[25] 11.7 0.0 1.2 9.2 11.8 14.0 2421 1
## y_new[26] 11.9 0.0 1.3 9.4 11.9 14.0 3654 1
## y_new[27] 12.0 0.0 1.3 9.4 12.0 14.2 2736 1
## y_new[28] 12.0 0.0 1.3 9.5 12.1 14.2 1976 1
## y_new[29] 12.2 0.0 1.3 9.7 12.2 14.4 1863 1
## y_new[30] 12.2 0.0 1.2 9.7 12.3 14.6 2049 1
## y_new[31] 12.3 0.0 1.3 9.7 12.4 14.7 1990 1
## y_new[32] 12.3 0.0 1.3 9.7 12.4 14.5 1679 1
## y_new[33] 12.4 0.0 1.3 9.8 12.5 14.6 1363 1
## y_new[34] 12.5 0.0 1.3 9.9 12.5 15.3 1614 1
## y_new[35] 12.5 0.0 1.3 10.0 12.6 15.0 1311 1
## y_new[36] 12.6 0.0 1.3 9.9 12.6 14.8 2385 1
## y_new[37] 12.6 0.0 1.3 10.0 12.7 15.0 2154 1
## y_new[38] 12.7 0.0 1.2 10.1 12.7 15.0 1766 1
## y_new[39] 12.7 0.0 1.3 9.9 12.8 15.0 1914 1
## y_new[40] 12.7 0.0 1.4 9.9 12.8 15.1 1755 1
## y_new[41] 12.8 0.0 1.4 10.2 12.8 14.9 1608 1
## y_new[42] 12.8 0.0 1.3 10.1 12.8 15.1 2016 1
## y_new[43] 12.8 0.0 1.3 10.2 12.9 15.1 1990 1
## y_new[44] 12.8 0.0 1.4 9.9 12.9 15.1 2076 1
## y_new[45] 12.9 0.0 1.3 10.2 12.9 15.2 1886 1
## y_new[46] 12.9 0.0 1.3 10.1 12.9 15.2 1656 1
## y_new[47] 12.9 0.0 1.3 10.5 13.0 15.3 2342 1
## y_new[48] 12.9 0.0 1.3 10.3 13.0 15.5 1553 1
## y_new[49] 12.9 0.0 1.4 10.2 13.0 15.2 1380 1
## y_new[50] 12.9 0.0 1.4 10.1 13.0 15.3 1331 1
## y_new[51] 12.9 0.0 1.3 10.2 13.0 15.3 1550 1
## y_new[52] 13.0 0.0 1.3 10.1 13.0 15.4 1321 1
## y_new[53] 13.0 0.0 1.4 10.2 13.0 15.5 1368 1
## y_new[54] 13.0 0.0 1.4 10.2 13.0 15.5 1912 1
## y_new[55] 13.0 0.0 1.3 10.3 13.0 15.5 2054 1
## y_new[56] 13.0 0.0 1.3 10.3 13.0 15.4 1469 1
## y_new[57] 13.0 0.0 1.4 10.2 13.0 15.4 1817 1
## y_new[58] 13.0 0.0 1.4 10.2 13.1 15.4 1772 1
## y_new[59] 13.0 0.0 1.4 10.3 13.0 15.6 1347 1
## y_new[60] 13.0 0.0 1.4 10.1 13.0 15.6 1735 1
## lp__ -0.2 0.1 2.0 -5.5 0.4 2.0 405 1
##
## Samples were drawn using NUTS(diag_e) at Sat Dec 16 11:40:13 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).予測分布の可視化を行いたいので、これまで通り、MCMCサンプルを取り出して、パーセンタイルを計算する。そのあとは可視化するだけ。
ms <- rstan::extract(fit)
qua <- apply(ms$y_new, 2, quantile, prob = c(0.025, 0.25, 0.5, 0.75, 0.975))
d_est <- data.frame(X = Time_new, t(qua), check.names = FALSE)
ggplot() +
theme_bw(base_size = 18) +
geom_ribbon(data = d_est, aes(x = X, ymin = `2.5%`, ymax = `97.5%`), fill = 'black', alpha = 1/6) +
geom_ribbon(data = d_est, aes(x = X, ymin = `25%`, ymax = `75%`), fill = 'black', alpha = 2/6) +
geom_line(data = d_est, aes(x = X, y = `50%`), size = 0.5) +
geom_point(data = d, aes(x = Time, y = Y), size = 3) +
labs(x = 'Time (hour)', y = 'Y') +
scale_x_continuous(breaks = d$Time, limit = c(0, 24)) +
ylim(-2.5, 16)
7.8 打ち切り
打ち切りデータと聞くと「生存時間分析」が思い浮かぶが、センサーや機器、時間の上限によって打ち切られているデータを分析することもベイズでは融通が効く。
書籍のデータにあるように検出限界値が25で、それより小さい場合は「不等号つき<25」で表されるようなデータがあったとする。このようなデータから平均と標準偏差を計算したい。不等号つきデータを削除して分析すれば、平均は過大評価され、不等号を外して25として扱えば、平均は過小評価される可能性がある。
d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap07/input/data-protein.txt')
d## Y
## 1 <25
## 2 32.3
## 3 <25
## 4 28.3
## 5 30.8
## 6 35.2モデルのメカニズムとして、平均\(\mu\)、標準偏差\(\sigma\)から潜在的な値\(y\)が生成されたと考える。閾値\(L\)とすると、\(y>L\)であれば、そのまま測定値\(Y\)となるが、\(y<L\)であれば測定値は\(<L\)となる。
\[ \begin{eqnarray} Not \ Censor \ \ \ \ Y[n] &\sim& Normal(\mu, \sigma_Y) \\ Censor \ \ \ \ y[n] &\sim& Normal(\mu, \sigma_Y) \ \ \ \ \ ※y[n] < L \end{eqnarray} \]
サンプリングするためには尤度が必要なので、尤度を考える。打ち切りがない場合は\(Normal(Y|\mu,\sigma_Y)\)でよく、打ち切りがある場合の尤度は\(y < L\)である確率\(Prob[y < L]\)である。尤度\(Prob[y < L]\)はよくよく見ると、これは分布関数なので、分布関数を尤度に利用すればよい。
\[ \begin{eqnarray} Prob[y < L] &=& \int_{-\infty}^{L} Normal(Y|\mu,\sigma_Y) \\ &=& \int_{-\infty}^{L} \frac{1}{\sqrt{2 \pi} \sigma} exp \left[ -\frac{1}{2}\left(\frac{y - \mu}{\sigma} \right)^2\right] dy \\ &=& \int_{-\infty}^{\frac{L-\mu}{\sigma}} \frac{1}{\sqrt{2 \pi}}exp \left[-\frac{z^{2}}{2} \right] dz \\ &=& \Phi \left( \frac{L-\mu}{\sigma} \right) \end{eqnarray} \]
2行目から3行目にかけて、積分変数を\(y\)から\(z\)に変更するために、下記の変数変換を行っている。積分範囲は\(y\)が\(-\infty \rightarrow L\)のとき、\(z\)は\(-\infty \rightarrow \frac{L-\mu}{\sigma}\)となる。
\[ \begin{eqnarray} z &=& \frac{y - \mu}{\sigma} \\ dy &=& \sigma dz \\ \end{eqnarray} \]
モデル式は下記の通り。打切りがある場合は、対数を取ってtargetに足し込めば良い。こうすることで、__lpに打ち切り分の対数尤度を加えることが可能になる。normal_lcdfは\(log \Phi \left( \frac{L-\mu}{\sigma}
\right)\)を意味する。
data {
int N_obs;
int N_cens;
real Y_obs[N_obs];
real L;
}
parameters {
real mu;
real<lower=0> s_Y;
}
model {
for (n in 1:N_obs){
Y_obs[n] ~ normal(mu, s_Y);
}
for (n in 1:N_cens){
target += normal_lcdf(L | mu, s_Y);
}
// 上3行のかわりに下記だけでも良い。N_censの分だけnormal_lcdfを加える
// target += N_cens * normal_lcdf(L | mu, s_Y);
}
データは打ち切りがある場合とない場合に分けて渡せばよい。
idx <- grep('<', d$Y)
Y_obs <- as.numeric(d[-idx, ])
L <- as.numeric(sub('<', '', d[idx,]))[1]
data <- list(N_obs = length(Y_obs), N_cens = length(idx), Y_obs = Y_obs, L = L)
data## $N_obs
## [1] 4
##
## $N_cens
## [1] 2
##
## $Y_obs
## [1] 32.3 28.3 30.8 35.2
##
## $L
## [1] 25ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果はこちら。サンプルが少ないので、平均\(\mu\)の信用区間がマイナスとなっているが、今回は問題ない。
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## mu 24.6 0.6 12.0 -6.1 26.9 38.7 383 1
## s_Y 16.2 1.0 21.1 4.2 10.5 59.3 454 1
## lp__ -10.4 0.1 1.5 -14.5 -9.9 -8.9 328 1
##
## Samples were drawn using NUTS(diag_e) at Sat Dec 16 11:40:33 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).