ベイズ統計04-普及モデル-
UPDATE: 2024-01-08 20:14:49.408378
はじめに
このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回は「社会科学のためのベイズ統計モデリング」の第8章を写経していく。基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
8.4.1 携帯電話の普及メカニズム
携帯電話の普及メカニズムに関するモデルを構築するのが今回の目的。携帯電話は、1985年に登場し、90年代に急速に普及した。その普及率の生成過程をモデル化する。
library(dplyr)
library(rstan)
library(ggplot2)
options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
d <- read.csv('https://raw.githubusercontent.com/HiroshiHamada/BMS/master/ch08/mobile-rate.csv')
d$time <- seq(1993, 2016, 1)
ggplot(d, aes(time, Y)) +
theme_bw(base_size = 15) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = seq(1993, 2016, 2)) +
labs(title = 'Mechanism of cell phone diffusion')
「曲線のあてはめ」と「データ生成過程からの関数系の導出」は異なる作業。普及の推移は関数\(f\)に似ているので、\(y=f(t)\)型で表現できそうだと考えるのが「あてはめ」という作業。あてはめはただ形が似ている以外の根拠はない。仮定を置くことで導出する。
- 契約数\(y\)は時間\(t\)の経過とともに継続的に増加する
- 契約数には上限\(m\)がある。\(k\)を増加の早さを決める。
- 未契約者はランダムに契約者と接触し、未契約者の一部が新たな契約者となる
3つ目は普及メカニズムを示している。最初の契約者は少ないが、契約者が未契約者に接触することで一定の割合で契約者に転換し、その契約者がさらに未契約者に転換することで一定の割合で契約者になっていくというメカニズムを想像している。ある時点での契約者の増え方は、その瞬間の契約者と契約の充足率に比例すると考える。契約者\(y\)の瞬間的な増分は、下記の微分方程式で表すことができる。
\[ \begin{eqnarray} \frac{dy}{dt} = ky \left( 1 - \frac{y}{m}\right) \end{eqnarray} \]
イメージしやすいように可視化しておく。\(k=1, m=100\)として、\(y\)と括弧の関係を見る。契約者が少ない間は増加しやすいが、充足上限に近づくと契約者の増加しにくくなる。
- \(y\)が0だとカッコ内は0.0となって0が返る
- \(y\)が10だとカッコ内は0.9となって9が返る
- \(y\)が20だとカッコ内は0.8となって16が返る
- \(y\)が50だとカッコ内は0.5となって25が返る
- \(y\)が90だとカッコ内は0.1となって9が返る
- \(y\)が100だとカッコ内は0.0となって0が返る
f <- function(y, k, m){
res <- k * y * (1 - (y/m))
return(res)
}
plot(
1:100,
f(y = 1:100, k = 1, m = 100),
type = 'l'
)
細かい導出は参考書を見るとして、この微分方程式から最終的に下記の関数が導出できる。
\[ \begin{eqnarray} Y &=& \frac{my_{0}}{(m-y_{0})e^{-kt}+y_{0}} + \epsilon \\ \end{eqnarray} \]
定数部分を\(\mu\)とすると、下記の通りモデル化できる。
\[ \begin{eqnarray} Y &=& \mu + \epsilon , \ \ \epsilon \sim Normal(0, \sigma)\\ \end{eqnarray} \]
下記のように書き直すとモデルがわかりやすくなる。つまり、平均\(\frac{my_{0}}{(m-y_{0})e^{-kt}+y_{0}}\)、標準偏差\(\sigma\)から\(Y\)が生成されると仮定している。
\[ \begin{eqnarray} Y &\sim& Normal \left( \frac{my_{0}}{(m-y_{0})e^{-kt}+y_{0}}, \sigma \right)\\ \end{eqnarray} \]
このようなデータ生成過程を考えず、単純に線形モデルを当てはめることもできる。
\[ \begin{eqnarray} Y &\sim& Normal \left( \sum_{i=0}^{k} \beta_{i}x_{i}, \sigma \right)\\ \end{eqnarray} \]
8.4.2 データとの対応
Stanでサンプリングするために必要なデータを作成する。
#予測分布用
n_pred <- 60
#予測分布用に説明変数を作成
t_pred <- seq(from = 1, to = 24, length = n_pred)
#契約者数の初期値
y0 <- 213.1
data <-
list(
n = nrow(d),
t = d$X,
Y = d$Y,
t_pred = t_pred,
n_pred = n_pred,
y0 = y0
)
data## $n
## [1] 24
##
## $t
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
##
## $Y
## [1] 213.1 433.1 1020.4 2087.7 3152.7 4153.0 5113.9 6094.2 6912.1
## [10] 7565.7 8152.0 8699.8 9179.2 9671.8 10272.5 10748.7 11218.3 11953.5
## [19] 12820.5 13604.0 14188.0 14879.0 15654.0 16344.0
##
## $t_pred
## [1] 1.000000 1.389831 1.779661 2.169492 2.559322 2.949153 3.338983
## [8] 3.728814 4.118644 4.508475 4.898305 5.288136 5.677966 6.067797
## [15] 6.457627 6.847458 7.237288 7.627119 8.016949 8.406780 8.796610
## [22] 9.186441 9.576271 9.966102 10.355932 10.745763 11.135593 11.525424
## [29] 11.915254 12.305085 12.694915 13.084746 13.474576 13.864407 14.254237
## [36] 14.644068 15.033898 15.423729 15.813559 16.203390 16.593220 16.983051
## [43] 17.372881 17.762712 18.152542 18.542373 18.932203 19.322034 19.711864
## [50] 20.101695 20.491525 20.881356 21.271186 21.661017 22.050847 22.440678
## [57] 22.830508 23.220339 23.610169 24.000000
##
## $n_pred
## [1] 60
##
## $y0
## [1] 213.1Stanのモデルは下記の通り。
data {
int n;
int t[n];
int n_pred;
real Y[n]; real y0;
real t_pred[n_pred];
}
parameters {
real <lower=16344,upper=20000> m;
real <lower=0,upper=5> k;
real<lower=0> sigma;
}
transformed parameters{
real mu[n];
for (i in 1:n) {
mu[i] = (m*y0)/((m-y0)*exp(-k*t[i]) + y0);
}
}
model{
for (i in 1:n) {
Y[i] ~ normal(mu[i], sigma);
}
}
generated quantities{
real y_pred[n_pred];
real log_lik[n];
for (i in 1:n_pred){
y_pred[i] = normal_rng((m*y0)/((m-y0)*exp(-k*t_pred[i])+y0), sigma);
}
for(i in 1:n){
log_lik[i] = normal_lpdf(Y[i]|mu[i], sigma);
}
}ここでは、stan_model()関数で最初にコンパイルしておいてから、
model <- stan_model('note_bayes04−001.stan')sampling()関数でサンプリングする。
fit <- sampling(object = model, data = data, seed = 1989)推定結果を確認する。
print(fit)## 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% 25% 50% 75%
## m 16757.27 9.16 437.52 16353.91 16461.46 16623.97 16898.79
## k 0.34 0.00 0.02 0.31 0.33 0.34 0.35
## sigma 2065.81 7.05 349.89 1511.72 1822.92 2019.12 2259.55
## mu[1] 298.18 0.12 5.33 288.12 294.63 297.93 301.56
## mu[2] 416.50 0.32 14.81 388.95 406.58 415.67 425.78
## mu[3] 580.27 0.66 30.70 523.97 559.60 578.33 599.28
## mu[4] 805.45 1.21 56.12 703.90 767.57 801.44 839.78
## mu[5] 1112.18 2.06 95.11 941.85 1047.80 1104.71 1169.86
## mu[6] 1524.62 3.29 152.35 1253.85 1421.14 1511.80 1616.45
## mu[7] 2069.61 5.02 232.33 1659.01 1911.72 2049.48 2209.52
## mu[8] 2773.16 7.29 337.56 2177.89 2543.38 2743.98 2975.40
## mu[9] 3654.63 10.05 466.05 2828.78 3336.45 3615.95 3936.39
## mu[10] 4718.60 13.10 608.77 3626.02 4303.86 4671.31 5093.27
## mu[11] 5947.04 16.06 748.98 4582.28 5438.24 5894.96 6415.63
## mu[12] 7295.60 18.48 865.09 5685.64 6707.14 7248.49 7843.78
## mu[13] 8697.66 19.92 937.07 6913.53 8066.34 8667.31 9306.85
## mu[14] 10076.59 20.14 953.38 8220.40 9441.45 10061.88 10704.79
## mu[15] 11361.64 19.15 914.74 9535.49 10757.02 11370.53 11969.68
## mu[16] 12501.09 17.20 833.02 10801.35 11963.60 12526.93 13062.34
## mu[17] 13467.99 14.69 726.49 11970.89 13012.05 13502.22 13958.10
## mu[18] 14258.41 12.02 614.42 12977.08 13879.14 14293.96 14673.94
## mu[19] 14885.03 9.58 513.38 13833.36 14566.24 14908.78 15232.93
## mu[20] 15369.73 7.66 435.09 14500.95 15095.92 15376.57 15653.06
## mu[21] 15737.55 6.48 385.04 15009.32 15493.57 15723.85 15958.55
## mu[22] 16012.60 6.06 361.51 15393.52 15784.36 15973.14 16188.85
## mu[23] 16216.02 6.20 357.26 15679.86 15988.03 16150.19 16368.92
## mu[24] 16365.21 6.59 363.81 15891.59 16128.72 16280.17 16509.56
## y_pred[1] 326.62 33.68 2120.19 -3789.28 -1056.71 336.03 1690.68
## y_pred[2] 390.33 33.88 2092.68 -3745.93 -936.17 373.55 1738.18
## y_pred[3] 369.47 33.67 2132.76 -4068.78 -995.63 334.98 1772.73
## y_pred[4] 474.24 33.07 2131.29 -3749.00 -895.00 403.29 1891.27
## y_pred[5] 486.93 33.81 2112.17 -3509.33 -937.71 484.73 1824.24
## y_pred[6] 574.58 33.02 2051.44 -3469.01 -755.06 585.26 1864.87
## y_pred[7] 653.67 34.42 2119.75 -3494.27 -700.28 646.07 1997.55
## y_pred[8] 772.83 33.82 2111.72 -3273.81 -605.07 771.42 2177.07
## y_pred[9] 821.78 33.06 2089.91 -3339.35 -518.68 810.29 2163.43
## y_pred[10] 935.94 33.10 2055.76 -3038.08 -397.80 899.99 2245.03
## y_pred[11] 1087.77 33.75 2097.65 -3054.45 -254.71 1155.53 2465.54
## y_pred[12] 1202.50 34.20 2084.03 -3013.79 -140.70 1247.36 2511.75
## y_pred[13] 1430.84 32.92 2087.29 -2701.44 121.43 1396.25 2808.97
## y_pred[14] 1585.13 35.31 2096.12 -2631.68 231.83 1600.05 2960.84
## y_pred[15] 1824.47 32.87 2059.14 -2130.68 472.80 1779.07 3128.32
## y_pred[16] 1956.97 33.59 2143.41 -2309.70 562.74 1990.54 3366.26
## y_pred[17] 2261.29 36.87 2126.29 -1891.25 868.17 2240.32 3677.41
## y_pred[18] 2522.90 32.86 2103.49 -1608.57 1144.83 2492.93 3901.93
## y_pred[19] 2773.84 35.12 2081.51 -1403.26 1416.35 2778.05 4153.21
## y_pred[20] 3122.50 34.13 2157.65 -1116.74 1736.08 3079.02 4462.58
## y_pred[21] 3450.06 35.82 2154.93 -827.57 2023.35 3486.56 4859.89
## y_pred[22] 3819.32 33.45 2102.70 -345.45 2405.57 3824.75 5224.02
## y_pred[23] 4253.67 36.03 2171.00 -49.81 2880.25 4183.26 5671.94
## y_pred[24] 4655.42 34.51 2219.50 309.80 3190.17 4659.69 6083.47
## y_pred[25] 5095.92 36.49 2248.17 697.85 3671.50 5109.23 6516.23
## y_pred[26] 5644.46 35.55 2183.77 1318.23 4284.00 5646.01 7060.20
## y_pred[27] 6036.20 36.95 2259.94 1598.03 4519.77 6047.83 7531.80
## y_pred[28] 6629.94 37.49 2261.45 2177.16 5192.35 6609.92 8032.08
## y_pred[29] 7238.74 35.97 2211.28 2954.31 5772.49 7261.58 8679.28
## y_pred[30] 7666.97 37.72 2287.51 3283.72 6168.27 7660.50 9136.40
## y_pred[31] 8289.82 38.75 2322.90 3752.63 6814.10 8220.31 9746.69
## y_pred[32] 8785.58 37.76 2288.26 4397.08 7293.01 8788.20 10286.14
## y_pred[33] 9333.14 39.91 2276.18 4732.35 7871.74 9375.66 10795.17
## y_pred[34] 9894.73 37.41 2286.62 5298.71 8450.53 9882.94 11404.00
## y_pred[35] 10365.09 40.07 2320.35 5711.37 8848.87 10360.62 11891.57
## y_pred[36] 10930.64 38.64 2301.50 6410.51 9405.04 10891.69 12445.74
## y_pred[37] 11372.25 37.97 2285.71 6869.53 9915.09 11363.32 12845.75
## y_pred[38] 11779.22 36.19 2259.20 7334.65 10300.27 11746.90 13312.06
## y_pred[39] 12308.21 37.98 2235.71 7853.73 10846.33 12309.37 13767.95
## y_pred[40] 12716.05 36.61 2215.10 8300.52 11240.12 12749.70 14167.80
## y_pred[41] 13089.40 36.05 2176.52 8750.92 11698.01 13124.00 14492.36
## y_pred[42] 13431.45 36.08 2232.02 8942.88 11994.64 13413.04 14878.80
## y_pred[43] 13777.07 37.73 2188.84 9425.28 12281.13 13815.77 15233.59
## y_pred[44] 14129.17 38.35 2192.13 9747.33 12724.68 14164.13 15577.76
## y_pred[45] 14381.19 36.54 2223.65 10089.85 12928.34 14426.23 15822.85
## y_pred[46] 14651.23 35.99 2169.71 10290.53 13277.40 14649.26 16051.70
## y_pred[47] 14853.13 33.90 2148.23 10637.45 13392.41 14831.01 16299.56
## y_pred[48] 15089.06 34.63 2132.31 10909.93 13714.12 15067.52 16449.98
## y_pred[49] 15285.67 33.16 2120.49 11047.93 13912.49 15311.83 16648.32
## y_pred[50] 15418.26 34.60 2163.83 11044.78 14025.93 15420.89 16841.99
## y_pred[51] 15552.41 33.53 2143.30 11389.48 14142.66 15504.62 16973.85
## y_pred[52] 15745.10 33.79 2151.80 11541.25 14327.78 15743.76 17135.29
## y_pred[53] 15836.72 34.26 2148.62 11583.79 14424.45 15826.05 17244.54
## y_pred[54] 15917.29 35.33 2153.49 11568.74 14540.45 15898.08 17324.07
## y_pred[55] 16002.38 33.08 2090.06 11834.91 14612.09 16035.05 17373.65
## y_pred[56] 16116.25 33.41 2102.84 12064.80 14742.94 16074.34 17443.78
## y_pred[57] 16130.80 33.93 2142.86 11886.53 14766.77 16143.93 17532.61
## y_pred[58] 16253.12 34.19 2068.01 12211.95 14873.57 16217.74 17606.31
## y_pred[59] 16253.31 32.14 2097.54 12199.65 14855.05 16252.88 17623.69
## y_pred[60] 16428.45 33.25 2108.52 12370.82 15059.71 16412.29 17739.50
## log_lik[1] -8.54 0.00 0.16 -8.89 -8.64 -8.53 -8.43
## log_lik[2] -8.54 0.00 0.16 -8.89 -8.64 -8.53 -8.43
## log_lik[3] -8.56 0.00 0.16 -8.90 -8.66 -8.55 -8.46
## log_lik[4] -8.75 0.00 0.10 -8.99 -8.81 -8.73 -8.67
## log_lik[5] -9.07 0.00 0.06 -9.20 -9.10 -9.06 -9.03
## log_lik[6] -9.42 0.00 0.15 -9.79 -9.49 -9.39 -9.31
## log_lik[7] -9.72 0.01 0.27 -10.36 -9.87 -9.68 -9.53
## log_lik[8] -9.95 0.01 0.38 -10.83 -10.16 -9.90 -9.68
## log_lik[9] -9.90 0.01 0.43 -10.89 -10.15 -9.86 -9.60
## log_lik[10] -9.60 0.01 0.43 -10.55 -9.86 -9.58 -9.30
## log_lik[11] -9.21 0.01 0.38 -10.06 -9.45 -9.17 -8.93
## log_lik[12] -8.87 0.01 0.30 -9.57 -9.04 -8.82 -8.64
## log_lik[13] -8.67 0.01 0.23 -9.21 -8.79 -8.64 -8.51
## log_lik[14] -8.66 0.01 0.24 -9.20 -8.79 -8.63 -8.51
## log_lik[15] -8.79 0.01 0.30 -9.50 -8.93 -8.73 -8.59
## log_lik[16] -9.01 0.01 0.37 -9.88 -9.21 -8.95 -8.74
## log_lik[17] -9.24 0.01 0.40 -10.14 -9.48 -9.20 -8.95
## log_lik[18] -9.25 0.01 0.34 -9.99 -9.47 -9.22 -9.01
## log_lik[19] -9.10 0.00 0.25 -9.63 -9.27 -9.09 -8.93
## log_lik[20] -8.95 0.00 0.19 -9.34 -9.07 -8.94 -8.82
## log_lik[21] -8.85 0.00 0.17 -9.22 -8.95 -8.84 -8.74
## log_lik[22] -8.71 0.00 0.17 -9.09 -8.81 -8.69 -8.59
## log_lik[23] -8.59 0.00 0.18 -8.96 -8.70 -8.57 -8.47
## log_lik[24] -8.55 0.00 0.17 -8.91 -8.66 -8.54 -8.44
## lp__ -182.68 0.03 1.34 -186.18 -183.29 -182.34 -181.70
## 97.5% n_eff Rhat
## m 17958.13 2281 1
## k 0.38 2140 1
## sigma 2889.87 2460 1
## mu[1] 309.15 2143 1
## mu[2] 447.29 2142 1
## mu[3] 644.67 2141 1
## mu[4] 924.16 2140 1
## mu[5] 1315.08 2140 1
## mu[6] 1851.86 2140 1
## mu[7] 2569.45 2142 1
## mu[8] 3499.28 2145 1
## mu[9] 4656.71 2152 1
## mu[10] 6015.47 2161 1
## mu[11] 7527.48 2174 1
## mu[12] 9099.38 2191 1
## mu[13] 10614.43 2212 1
## mu[14] 11969.12 2241 1
## mu[15] 13130.98 2282 1
## mu[16] 14055.72 2344 1
## mu[17] 14799.67 2446 1
## mu[18] 15386.05 2613 1
## mu[19] 15851.78 2875 1
## mu[20] 16228.86 3227 1
## mu[21] 16564.89 3528 1
## mu[22] 16860.76 3553 1
## mu[23] 17134.93 3324 1
## mu[24] 17314.95 3045 1
## y_pred[1] 4596.17 3964 1
## y_pred[2] 4717.67 3814 1
## y_pred[3] 4578.32 4012 1
## y_pred[4] 4677.66 4153 1
## y_pred[5] 4748.60 3904 1
## y_pred[6] 4753.53 3859 1
## y_pred[7] 4855.92 3793 1
## y_pred[8] 4900.89 3898 1
## y_pred[9] 5052.86 3996 1
## y_pred[10] 5057.30 3857 1
## y_pred[11] 5229.09 3862 1
## y_pred[12] 5291.78 3713 1
## y_pred[13] 5548.31 4021 1
## y_pred[14] 5593.28 3523 1
## y_pred[15] 6055.39 3924 1
## y_pred[16] 6264.55 4071 1
## y_pred[17] 6535.19 3325 1
## y_pred[18] 6776.35 4097 1
## y_pred[19] 6882.23 3512 1
## y_pred[20] 7476.98 3998 1
## y_pred[21] 7748.83 3620 1
## y_pred[22] 7983.00 3951 1
## y_pred[23] 8540.33 3631 1
## y_pred[24] 9053.16 4137 1
## y_pred[25] 9651.50 3796 1
## y_pred[26] 9962.07 3774 1
## y_pred[27] 10456.32 3740 1
## y_pred[28] 11146.84 3638 1
## y_pred[29] 11591.57 3780 1
## y_pred[30] 12252.14 3678 1
## y_pred[31] 13030.13 3593 1
## y_pred[32] 13349.12 3672 1
## y_pred[33] 13771.66 3253 1
## y_pred[34] 14310.76 3735 1
## y_pred[35] 14892.07 3354 1
## y_pred[36] 15408.86 3547 1
## y_pred[37] 15910.36 3623 1
## y_pred[38] 16091.45 3897 1
## y_pred[39] 16649.21 3466 1
## y_pred[40] 17086.24 3660 1
## y_pred[41] 17279.70 3646 1
## y_pred[42] 17946.65 3827 1
## y_pred[43] 18147.49 3365 1
## y_pred[44] 18407.97 3267 1
## y_pred[45] 18737.81 3704 1
## y_pred[46] 18849.21 3634 1
## y_pred[47] 19014.61 4016 1
## y_pred[48] 19363.46 3791 1
## y_pred[49] 19461.73 4089 1
## y_pred[50] 19699.29 3911 1
## y_pred[51] 19807.45 4086 1
## y_pred[52] 20069.35 4056 1
## y_pred[53] 20047.09 3933 1
## y_pred[54] 20176.69 3715 1
## y_pred[55] 20166.45 3993 1
## y_pred[56] 20328.22 3960 1
## y_pred[57] 20353.68 3987 1
## y_pred[58] 20581.61 3659 1
## y_pred[59] 20451.22 4259 1
## y_pred[60] 20729.24 4021 1
## log_lik[1] -8.24 2465 1
## log_lik[2] -8.24 2464 1
## log_lik[3] -8.28 2459 1
## log_lik[4] -8.59 2403 1
## log_lik[5] -8.96 2184 1
## log_lik[6] -9.18 2443 1
## log_lik[7] -9.30 2491 1
## log_lik[8] -9.33 2519 1
## log_lik[9] -9.16 2550 1
## log_lik[10] -8.88 2570 1
## log_lik[11] -8.58 2511 1
## log_lik[12] -8.40 2261 1
## log_lik[13] -8.30 1832 1
## log_lik[14] -8.31 1878 1
## log_lik[15] -8.36 2239 1
## log_lik[16] -8.47 2533 1
## log_lik[17] -8.62 2706 1
## log_lik[18] -8.68 2870 1
## log_lik[19] -8.65 3111 1
## log_lik[20] -8.60 3219 1
## log_lik[21] -8.55 3043 1
## log_lik[22] -8.42 2584 1
## log_lik[23] -8.27 2334 1
## log_lik[24] -8.25 2328 1
## lp__ -181.15 1484 1
##
## Samples were drawn using NUTS(diag_e) at Mon Jan 8 20:15:15 2024.
## 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).事後分布を可視化しておく。
stan_plot(
fit,
pars = c('m', 'k', 'sigma'),
point_est = 'mean',
ci_level = 0.95,
outer_level = 1.00,
show_density = TRUE,
fill_color = 'grey') +
theme_bw()
平均と標準偏差の事後分布を使って、予測分布を可視化する。
ms <- rstan::extract(fit)
qua <- apply(ms$y_pred, 2, quantile, prob = c(0.025, 0.25, 0.5, 0.75, 0.975))
d_est <- data.frame(X = t_pred, t(qua), check.names = FALSE)
ggplot() +
theme_bw(base_size = 15) +
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 = X, y = Y), size = 3) +
labs(x = 'Time (hour)', y = 'Y') +
scale_x_continuous(labels = seq(1993, 2016, 2), breaks = seq(1, 24, 2)) +
labs(title = 'Mechanism of cell phone diffusion')
参考書に書かれており通り、普及メカニズムを考慮したモデルと線形モデルでは、線形モデルのほうが予測精度が良い。ただ、データ生成過程を考慮してモデルを作ったのであれば、モデルを改良することができるため、\(m,k\)などを時点で変化するモデルなどに改良すれば、線形モデルでも優れたモデルに改良できる。
線形モデルの推定を可視化を行っておく。ここでは、stan_model()関数で最初にコンパイルしておいてから、
model2 <- stan_model('note_bayes04−002.stan')sampling()関数でサンプリングして、線形モデルでも予測分布を可視化しておく。
fit2 <- sampling(object = model2, data = data, seed = 1989)
ms2 <- rstan::extract(fit2)
qua2 <- apply(ms2$y_pred, 2, quantile, prob = c(0.025, 0.25, 0.5, 0.75, 0.975))
d_est2 <- data.frame(X = t_pred, t(qua2), check.names = FALSE)
ggplot() +
theme_bw(base_size = 15) +
geom_ribbon(data = d_est2, aes(x = X, ymin = `2.5%`, ymax = `97.5%`), fill = 'black', alpha = 1/6) +
geom_ribbon(data = d_est2, aes(x = X, ymin = `25%`, ymax = `75%`), fill = 'black', alpha = 2/6) +
geom_line(data = d_est2, aes(x = X, y = `50%`), size = 0.5) +
geom_point(data = d, aes(x = X, y = Y), size = 3) +
labs(x = 'Time (hour)', y = 'Y') +
scale_x_continuous(labels = seq(1993, 2016, 2), breaks = seq(1, 24, 2)) +
labs(title = 'Mechanism of cell phone diffusion')