ボートレースの勝率を推定する
UPDATE: 2024-01-20 16:22:39.827141
はじめに
このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回は下記のブログに記載されているテニス選手の強さを推定する内容を参考にさせていただき、写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
第8章の「動的一般化線形モデル：二項分布を仮定した例」を写経する。
8.2 GLMの復習
一般化線形モデルは、確率分布、線形予測子、リンク関数を部品とするモデルでこれらの3つを考えてモデリングを行う。ここで扱うデータはボートレースの勝ち負け(0/1)が記録されているデータ。
- 確率分布: 二項分布
- 線形予測子: 真の状態\(\mu[t]\)
- リンク関数: ロジット関数(ロジスティック関数の逆関数)を利用。二項分布のパラメタは0-1の範囲しか取れないため。
動的一般化線形モデルは、線形予測子が動的なものに対応している一般化線形モデルの拡張とも言える。ここでは下記の二項分布を仮定した動的一般化線形モデルを利用する。
\[ \begin{eqnarray} \mu[t] &=& \mu[t-1] + \epsilon_{\mu} &\quad \epsilon_{\mu} \sim Normal(0, \sigma_{\mu}) \\ y[t] &\sim& Bernoulli(logistic(\mu[t]) \\ \text{Simple Version} \\ \mu[t] &\sim& Normal(\mu[t-1], \sigma_{\mu}) \\ y[t] &\sim& Bernoulli(logistic(\mu[t]) \\ \end{eqnarray} \]
説明変数などがあれば下記のように時変係数モデルを参考に拡張する。
\[ \begin{eqnarray} \mu[t] &=& \mu[t-1] + \epsilon_{\mu} &\quad \epsilon_{\mu} \sim Normal(0, \sigma_{\mu}) \\ \beta[t] &=& \beta[t-1] + \epsilon_{\beta} &\quad \epsilon_{\beta} \sim Normal(0, \sigma_{\beta}) \\ \theta[t] &=& logistic(\mu[t] + \beta[t] \cdot x[t]) \\ y[t] &\sim& Bernoulli(\theta[t]) \end{eqnarray} \]
ここではKFASパッケージに含まれるboatデータを利用する。boatデータはオックスフォード大学(が勝利した場合は0)とケンブリッジ大学(が勝利した場合は1)の間で毎年行われるボートレースの結果を記録しており、183個のうち28個が欠損値という時系列データ。
library(tidyverse)
library(rstan)
#library(brms)
#library(bayesplot)
#library(patchwork)
library(KFAS)
options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
data('boat')
boat## Time Series:
## Start = 1829
## End = 2011
## Frequency = 1
## [1] 0 NA NA NA NA NA NA 1 NA NA 1 1 1 0 NA NA 1 1 NA NA 1 NA NA 0 NA
## [26] 0 NA 1 0 1 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 NA 0
## [51] 1 0 0 0 0 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 0 1 1
## [76] 1 0 1 1 1 0 0 0 0 0 1 NA NA NA NA NA 1 1 1 0 1 1 1 1 1
## [101] 1 1 1 1 1 1 1 1 0 0 1 NA NA NA NA NA NA 0 1 1 1 1 1 0 1
## [126] 0 1 1 1 1 0 0 1 1 0 1 0 0 0 1 1 1 1 1 1 0 1 0 0 0
## [151] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 0 1 0 0
## [176] 1 0 0 1 0 0 1 0データを見るとわかるが、2個目から7個目までは欠損値で、それ以降もところどころ欠損値がある。
which(!is.na(boat)) ## [1] 1 8 11 12 13 14 17 18 21 24 26 28 29 30 31 32 33 34
## [19] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 50 51 52 53
## [37] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [55] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 92 93 94
## [73] 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 118
## [91] 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
## [109] 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
## [127] 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
## [145] 173 174 175 176 177 178 179 180 181 182 183データを準備する。欠損値があるデータでも扱えるのが状態空間モデルの強み。データを渡す場合は、欠損値を除いたデータy、長さlen_obs、インデックスobs_noが必要になる。
boat_omit_NA <- na.omit(as.numeric(boat))
# データの準備
data_list <- list(
T = length(boat),
len_obs = length(boat_omit_NA), # 観測値があるデータのみ
y = boat_omit_NA, # 欠損値を除外したデータのみ
obs_no = which(!is.na(boat)) # 観測値があるデータのインデックス
)
data_list## $T
## [1] 183
##
## $len_obs
## [1] 155
##
## $y
## [1] 0 1 1 1 1 0 1 1 1 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 1 0 0 0
## [38] 0 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 1 1 1 0 0 0 0 0 1 1 1 1 0 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 1 1 1 1 1 0 1 0 1 1 1 1 0 0 1 1 0 1 0 0 0
## [112] 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 0 1 0 0 1
## [149] 0 0 1 0 0 1 0
## attr(,"na.action")
## [1] 2 3 4 5 6 7 9 10 15 16 19 20 22 23 25 27 49 87 88
## [20] 89 90 91 112 113 114 115 116 117
## attr(,"class")
## [1] "omit"
##
## $obs_no
## [1] 1 8 11 12 13 14 17 18 21 24 26 28 29 30 31 32 33 34
## [19] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 50 51 52 53
## [37] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [55] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 92 93 94
## [73] 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 118
## [91] 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
## [109] 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
## [127] 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
## [145] 173 174 175 176 177 178 179 180 181 182 183Stanモデルはこちら。
data {
int T; // データ取得期間の長さ
int len_obs; // 観測値が得られた個数
int y[len_obs]; // 観測値
int obs_no[len_obs]; // 観測値が得られた時点
}
parameters {
real mu[T]; // 状態の推定値
real<lower=0> s_mu; // 過程誤差の標準偏差
}
model {
s_mu ~ student_t(3, 0, 10);
// 状態方程式に従い、状態が遷移する
for(t in 2:T) {
mu[t] ~ normal(mu[t-1], s_mu);
}
// 観測方程式に従い、観測値が得られるが「観測値が得られた時点」でのみ実行
for(t in 1:len_obs) {
y[t] ~ bernoulli_logit(mu[obs_no[t]]);
}
}
generated quantities{
real probs[T]; // 推定された勝率
probs = inv_logit(mu);
}
モデルの部分が気になるので、深掘りしておく。観測方程式に少し違和感があるので、boatデータと照らし合わせて深掘りしておく。
model {
// 状態方程式に従い、状態が遷移する
for(t in 2:T) {
mu[t] ~ normal(mu[t-1], s_mu);
}
// 観測方程式に従い、観測値が得られるが「観測値が得られた時点」でのみ実行
for(t in 1:len_obs) {
y[t] ~ bernoulli_logit(mu[obs_no[t]]);
}
}
// 状態方程式
// t=2: mu[2] ~ normal(mu[1], s_mu);
// t=3: mu[3] ~ normal(mu[2], s_mu);
// t=4: mu[4] ~ normal(mu[3], s_mu);
// t=5: mu[5] ~ normal(mu[4], s_mu);
// 状態方程式は欠損値があった時点であっても、前の時点の真の状態から生成される
// 観測方程式
// t=2: y[2] ~ bernoulli_logit(mu[obs_no[2]]) -> y[2] ~ bernoulli_logit(mu[8])
// t=3: y[3] ~ bernoulli_logit(mu[obs_no[3]]) -> y[3] ~ bernoulli_logit(mu[11])
// t=4: y[4] ~ bernoulli_logit(mu[obs_no[4]]) -> y[4] ~ bernoulli_logit(mu[12])
// t=5: y[5] ~ bernoulli_logit(mu[obs_no[5]]) -> y[5] ~ bernoulli_logit(mu[13])
// y = boat_omit_NA より「欠損値を除外したデータ」なので、
// boatのデータ
// boat[1] - y[1]: 0 <- boat[1]を意味する
// boat[2] ------: NA
// boat[3] ------: NA
// boat[4] ------: NA
// boat[5] ------: NA
// boat[6] ------: NA
// boat[7] ------: NA
// boat[8] - y[2]: 1 <- boat[8]を意味し、y[2]はmu[8]から生成されている関係性になる。ここでは、stan_model()関数で最初にコンパイルしておいてから、
model <- stan_model('model-ts.stan')sampling()関数でサンプリングする。
fit <- sampling(object = model, data = data_list, seed = 1989)推定結果を確認する。
print(fit, prob = c(0.025, 0.5, 0.975), digits_summary = 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
## mu[1] -0.5 0.2 2.0 -5.3 -0.2 2.7 171 1.0
## mu[2] -0.2 0.1 2.0 -4.6 0.0 3.0 263 1.0
## mu[3] 0.0 0.1 1.9 -4.2 0.2 3.3 709 1.0
## mu[4] 0.2 0.0 1.9 -3.7 0.3 3.6 1736 1.0
## mu[5] 0.5 0.0 1.8 -3.0 0.5 3.9 1978 1.0
## mu[6] 0.7 0.0 1.7 -2.4 0.7 4.2 1948 1.0
## mu[7] 1.0 0.0 1.6 -1.9 0.9 4.3 1898 1.0
## mu[8] 1.2 0.1 1.5 -1.2 1.0 4.5 272 1.0
## mu[9] 1.3 0.1 1.5 -1.0 1.1 4.7 225 1.0
## mu[10] 1.4 0.1 1.4 -1.0 1.2 4.7 154 1.0
## mu[11] 1.5 0.1 1.4 -0.7 1.3 4.6 97 1.0
## mu[12] 1.4 0.1 1.3 -0.6 1.3 4.4 108 1.0
## mu[13] 1.2 0.1 1.2 -0.8 1.1 3.9 208 1.0
## mu[14] 0.9 0.0 1.1 -1.1 0.9 3.2 1287 1.0
## mu[15] 1.0 0.0 1.2 -1.2 0.9 3.7 604 1.0
## mu[16] 1.2 0.1 1.3 -1.1 1.0 4.0 222 1.0
## mu[17] 1.3 0.1 1.3 -0.8 1.2 4.2 95 1.0
## mu[18] 1.3 0.1 1.3 -0.8 1.1 4.4 90 1.0
## mu[19] 1.1 0.1 1.4 -1.0 0.9 4.3 139 1.0
## mu[20] 1.0 0.1 1.3 -1.2 0.8 4.0 187 1.0
## mu[21] 0.8 0.1 1.2 -1.3 0.7 3.7 279 1.0
## mu[22] 0.4 0.0 1.2 -1.9 0.4 3.1 2035 1.0
## mu[23] 0.1 0.0 1.2 -2.4 0.1 2.5 2397 1.0
## mu[24] -0.3 0.1 1.2 -2.9 -0.2 1.7 375 1.0
## mu[25] -0.4 0.1 1.2 -2.9 -0.3 1.7 371 1.0
## mu[26] -0.5 0.1 1.1 -3.0 -0.4 1.4 320 1.0
## mu[27] -0.3 0.0 1.1 -2.6 -0.3 1.8 2756 1.0
## mu[28] -0.2 0.0 1.0 -2.1 -0.2 1.9 1766 1.0
## mu[29] -0.4 0.0 0.9 -2.4 -0.4 1.4 3114 1.0
## mu[30] -0.3 0.0 1.0 -2.2 -0.3 1.6 2272 1.0
## mu[31] -0.6 0.0 0.9 -2.6 -0.6 1.2 3107 1.0
## mu[32] -0.7 0.0 0.9 -2.7 -0.6 1.2 3305 1.0
## mu[33] -1.1 0.1 1.1 -3.5 -1.0 0.6 94 1.0
## mu[34] -1.5 0.2 1.2 -4.4 -1.3 0.3 63 1.1
## mu[35] -1.6 0.2 1.3 -4.7 -1.4 0.2 49 1.1
## mu[36] -1.7 0.2 1.4 -5.2 -1.5 0.2 45 1.1
## mu[37] -1.7 0.2 1.4 -5.0 -1.4 0.2 42 1.1
## mu[38] -1.6 0.2 1.3 -4.8 -1.3 0.3 46 1.1
## mu[39] -1.4 0.2 1.3 -4.4 -1.2 0.4 64 1.1
## mu[40] -1.1 0.1 1.1 -3.7 -0.9 0.7 111 1.0
## mu[41] -0.6 0.0 1.0 -2.7 -0.5 1.1 834 1.0
## mu[42] 0.1 0.1 1.0 -1.6 0.0 2.4 314 1.0
## mu[43] 0.5 0.1 1.1 -1.2 0.4 3.1 89 1.0
## mu[44] 0.7 0.2 1.2 -1.0 0.6 3.6 60 1.1
## mu[45] 0.8 0.2 1.2 -1.0 0.6 3.5 55 1.1
## mu[46] 0.6 0.1 1.1 -1.1 0.4 3.1 126 1.0
## mu[47] 0.2 0.0 0.9 -1.6 0.2 2.2 1936 1.0
## mu[48] 0.2 0.0 1.0 -1.5 0.1 2.4 594 1.0
## mu[49] 0.0 0.0 1.0 -2.0 -0.1 2.0 4659 1.0
## mu[50] -0.3 0.0 0.9 -2.4 -0.3 1.5 3555 1.0
## mu[51] -0.3 0.0 0.9 -2.2 -0.3 1.5 4034 1.0
## mu[52] -0.7 0.1 1.0 -3.0 -0.6 1.0 230 1.0
## mu[53] -0.9 0.1 1.1 -3.4 -0.7 0.8 103 1.0
## mu[54] -0.8 0.1 1.0 -3.3 -0.7 0.8 149 1.0
## mu[55] -0.6 0.1 0.9 -2.7 -0.5 1.0 359 1.0
## mu[56] -0.2 0.0 0.9 -1.9 -0.2 1.6 3731 1.0
## mu[57] -0.2 0.0 0.9 -1.9 -0.1 1.6 4252 1.0
## mu[58] 0.3 0.1 1.0 -1.4 0.2 2.4 191 1.0
## mu[59] 0.4 0.1 1.1 -1.3 0.3 3.0 93 1.0
## mu[60] 0.3 0.1 1.1 -1.4 0.2 2.7 123 1.0
## mu[61] -0.1 0.0 1.0 -1.8 -0.1 2.1 410 1.0
## mu[62] -0.7 0.0 1.0 -2.9 -0.6 1.1 514 1.0
## mu[63] -1.1 0.1 1.1 -3.8 -0.9 0.7 127 1.0
## mu[64] -1.4 0.1 1.3 -4.4 -1.2 0.4 76 1.1
## mu[65] -1.6 0.2 1.4 -4.9 -1.3 0.3 42 1.1
## mu[66] -1.6 0.2 1.4 -5.2 -1.3 0.3 44 1.1
## mu[67] -1.6 0.2 1.4 -4.9 -1.3 0.3 42 1.1
## mu[68] -1.4 0.2 1.3 -4.6 -1.2 0.4 53 1.1
## mu[69] -1.2 0.1 1.2 -4.1 -1.0 0.6 77 1.0
## mu[70] -0.7 0.1 1.0 -3.0 -0.6 0.9 289 1.0
## mu[71] -0.1 0.0 0.9 -1.9 -0.1 1.8 1778 1.0
## mu[72] 0.2 0.0 0.9 -1.6 0.1 2.2 745 1.0
## mu[73] 0.2 0.0 0.9 -1.6 0.1 2.0 3210 1.0
## mu[74] 0.5 0.1 1.0 -1.1 0.4 2.9 207 1.0
## mu[75] 0.7 0.1 1.1 -1.0 0.6 3.2 124 1.0
## mu[76] 0.7 0.1 1.0 -1.0 0.6 3.0 194 1.0
## mu[77] 0.4 0.0 0.9 -1.3 0.4 2.4 3725 1.0
## mu[78] 0.6 0.1 1.0 -1.1 0.5 2.8 318 1.0
## mu[79] 0.5 0.1 1.0 -1.2 0.4 2.8 270 1.0
## mu[80] 0.3 0.0 0.9 -1.5 0.2 2.2 3752 1.0
## mu[81] -0.3 0.1 1.0 -2.6 -0.2 1.4 189 1.0
## mu[82] -0.6 0.1 1.1 -3.2 -0.4 1.1 75 1.1
## mu[83] -0.7 0.2 1.2 -3.4 -0.5 1.1 55 1.1
## mu[84] -0.6 0.1 1.2 -3.3 -0.5 1.2 70 1.1
## mu[85] -0.3 0.1 1.1 -2.7 -0.2 1.4 106 1.0
## mu[86] 0.2 0.0 1.1 -1.9 0.3 2.4 4069 1.0
## mu[87] 0.5 0.0 1.3 -2.1 0.5 3.1 3361 1.0
## mu[88] 0.7 0.0 1.4 -1.8 0.6 3.7 2186 1.0
## mu[89] 0.9 0.0 1.5 -1.7 0.8 4.1 1655 1.0
## mu[90] 1.2 0.1 1.4 -1.4 1.0 4.4 660 1.0
## mu[91] 1.4 0.1 1.4 -0.9 1.2 4.6 400 1.0
## mu[92] 1.6 0.1 1.3 -0.4 1.4 4.8 212 1.0
## mu[93] 1.7 0.1 1.2 -0.3 1.5 4.7 145 1.0
## mu[94] 1.7 0.1 1.2 -0.2 1.5 4.4 337 1.0
## mu[95] 1.5 0.0 1.1 -0.5 1.4 3.8 2221 1.0
## mu[96] 1.9 0.1 1.2 0.0 1.7 4.5 268 1.0
## mu[97] 2.2 0.1 1.3 0.2 2.0 5.3 139 1.0
## mu[98] 2.4 0.2 1.4 0.4 2.2 5.7 70 1.0
## mu[99] 2.5 0.2 1.5 0.4 2.3 6.2 60 1.1
## mu[100] 2.7 0.2 1.6 0.5 2.4 6.4 65 1.1
## mu[101] 2.7 0.2 1.6 0.5 2.4 6.6 59 1.1
## mu[102] 2.7 0.2 1.6 0.5 2.4 6.7 62 1.1
## mu[103] 2.6 0.2 1.6 0.5 2.3 6.6 66 1.1
## mu[104] 2.5 0.2 1.5 0.4 2.3 6.3 67 1.1
## mu[105] 2.4 0.1 1.4 0.3 2.1 5.9 99 1.0
## mu[106] 2.1 0.1 1.3 0.2 1.9 5.2 120 1.0
## mu[107] 1.8 0.1 1.2 0.0 1.7 4.7 192 1.0
## mu[108] 1.4 0.0 1.0 -0.4 1.4 3.6 2676 1.0
## mu[109] 0.8 0.1 1.0 -1.4 0.9 2.7 283 1.0
## mu[110] 0.7 0.1 1.1 -1.8 0.8 2.6 207 1.0
## mu[111] 0.9 0.0 1.1 -1.3 0.9 3.2 2536 1.0
## mu[112] 0.9 0.0 1.2 -1.7 0.9 3.4 2661 1.0
## mu[113] 0.9 0.0 1.3 -1.9 0.9 3.5 2399 1.0
## mu[114] 0.9 0.0 1.4 -2.1 0.9 3.6 2384 1.0
## mu[115] 0.8 0.0 1.4 -2.1 0.9 3.5 2196 1.0
## mu[116] 0.8 0.0 1.4 -2.1 0.9 3.4 2590 1.0
## mu[117] 0.8 0.0 1.3 -1.8 0.9 3.3 2603 1.0
## mu[118] 0.8 0.0 1.1 -1.5 0.8 2.9 1337 1.0
## mu[119] 1.2 0.0 1.0 -0.7 1.1 3.5 1081 1.0
## mu[120] 1.4 0.1 1.1 -0.4 1.3 3.9 262 1.0
## mu[121] 1.5 0.1 1.1 -0.3 1.3 4.0 155 1.0
## mu[122] 1.4 0.1 1.1 -0.4 1.3 3.8 161 1.0
## mu[123] 1.2 0.0 1.0 -0.6 1.1 3.4 430 1.0
## mu[124] 0.8 0.0 0.9 -1.0 0.8 2.7 4325 1.0
## mu[125] 0.9 0.0 0.9 -0.9 0.9 2.9 4668 1.0
## mu[126] 0.7 0.0 0.9 -1.1 0.8 2.6 3131 1.0
## mu[127] 1.0 0.0 0.9 -0.6 1.0 3.1 735 1.0
## mu[128] 1.1 0.1 1.0 -0.5 1.0 3.5 249 1.0
## mu[129] 1.1 0.1 1.0 -0.6 0.9 3.5 213 1.0
## mu[130] 0.8 0.0 0.9 -0.9 0.7 2.9 1162 1.0
## mu[131] 0.3 0.0 0.9 -1.7 0.3 2.0 888 1.0
## mu[132] 0.2 0.0 1.0 -1.8 0.2 2.0 789 1.0
## mu[133] 0.4 0.0 0.9 -1.3 0.4 2.3 4510 1.0
## mu[134] 0.4 0.0 0.9 -1.4 0.4 2.4 3126 1.0
## mu[135] 0.1 0.0 0.9 -1.9 0.1 1.8 2514 1.0
## mu[136] 0.1 0.0 0.9 -1.8 0.1 1.9 4692 1.0
## mu[137] -0.2 0.1 1.0 -2.5 -0.1 1.5 296 1.0
## mu[138] -0.2 0.1 1.0 -2.5 -0.1 1.5 188 1.0
## mu[139] 0.0 0.0 0.9 -2.0 0.1 1.7 696 1.0
## mu[140] 0.5 0.0 1.0 -1.2 0.5 2.7 664 1.0
## mu[141] 0.9 0.1 1.1 -0.8 0.7 3.4 114 1.0
## mu[142] 1.0 0.1 1.2 -0.8 0.8 3.8 71 1.1
## mu[143] 1.0 0.2 1.2 -0.8 0.8 3.9 55 1.1
## mu[144] 0.8 0.2 1.2 -0.9 0.6 3.6 60 1.1
## mu[145] 0.4 0.1 1.0 -1.2 0.3 2.8 173 1.0
## mu[146] -0.1 0.0 0.9 -1.9 -0.1 1.8 4158 1.0
## mu[147] -0.4 0.0 0.9 -2.2 -0.4 1.5 4139 1.0
## mu[148] -1.0 0.1 1.0 -3.4 -0.9 0.7 278 1.0
## mu[149] -1.4 0.1 1.2 -4.4 -1.2 0.5 87 1.0
## mu[150] -1.7 0.2 1.4 -4.9 -1.5 0.1 67 1.1
## mu[151] -2.0 0.2 1.5 -5.5 -1.7 0.0 53 1.1
## mu[152] -2.1 0.2 1.5 -5.9 -1.8 0.0 50 1.1
## mu[153] -2.2 0.2 1.5 -5.9 -1.9 -0.1 49 1.1
## mu[154] -2.1 0.2 1.5 -5.7 -1.9 -0.1 45 1.1
## mu[155] -2.1 0.2 1.4 -5.7 -1.8 0.0 47 1.1
## mu[156] -1.9 0.2 1.3 -5.1 -1.7 0.0 70 1.1
## mu[157] -1.7 0.1 1.2 -4.4 -1.6 0.2 108 1.0
## mu[158] -1.4 0.0 1.1 -3.8 -1.3 0.4 457 1.0
## mu[159] -1.5 0.1 1.1 -4.1 -1.4 0.2 101 1.0
## mu[160] -1.6 0.1 1.2 -4.5 -1.4 0.2 67 1.1
## mu[161] -1.5 0.1 1.2 -4.6 -1.3 0.3 75 1.1
## mu[162] -1.3 0.1 1.2 -4.2 -1.1 0.4 95 1.0
## mu[163] -1.0 0.1 1.1 -3.7 -0.8 0.8 153 1.0
## mu[164] -0.5 0.0 1.0 -2.7 -0.4 1.2 703 1.0
## mu[165] 0.2 0.1 1.0 -1.5 0.1 2.6 242 1.0
## mu[166] 0.7 0.1 1.2 -1.0 0.5 3.4 80 1.1
## mu[167] 1.0 0.2 1.3 -0.9 0.8 4.0 49 1.1
## mu[168] 1.1 0.2 1.4 -0.8 0.9 4.4 40 1.1
## mu[169] 1.1 0.2 1.3 -0.8 0.9 4.2 35 1.1
## mu[170] 1.0 0.2 1.2 -0.9 0.8 3.9 44 1.1
## mu[171] 0.7 0.1 1.1 -1.0 0.6 3.2 90 1.0
## mu[172] 0.2 0.0 0.9 -1.6 0.2 2.2 3739 1.0
## mu[173] 0.1 0.0 0.9 -1.7 0.1 2.1 3296 1.0
## mu[174] -0.3 0.0 0.9 -2.3 -0.3 1.4 822 1.0
## mu[175] -0.5 0.0 1.0 -2.6 -0.4 1.2 517 1.0
## mu[176] -0.3 0.0 0.9 -2.3 -0.3 1.4 4119 1.0
## mu[177] -0.6 0.0 1.0 -2.7 -0.5 1.2 523 1.0
## mu[178] -0.7 0.0 1.0 -2.9 -0.6 1.1 441 1.0
## mu[179] -0.5 0.0 1.0 -2.6 -0.5 1.3 4119 1.0
## mu[180] -0.8 0.1 1.0 -3.0 -0.7 1.0 371 1.0
## mu[181] -0.8 0.1 1.1 -3.1 -0.7 1.0 449 1.0
## mu[182] -0.6 0.0 1.1 -2.9 -0.6 1.6 3966 1.0
## mu[183] -0.8 0.0 1.2 -3.5 -0.7 1.4 696 1.0
## s_mu 0.7 0.1 0.4 0.2 0.6 1.8 19 1.2
## probs[1] 0.4 0.0 0.3 0.0 0.4 0.9 289 1.0
## probs[2] 0.5 0.0 0.3 0.0 0.5 1.0 503 1.0
## probs[3] 0.5 0.0 0.3 0.0 0.5 1.0 1121 1.0
## probs[4] 0.5 0.0 0.3 0.0 0.6 1.0 1633 1.0
## probs[5] 0.6 0.0 0.3 0.0 0.6 1.0 2273 1.0
## probs[6] 0.6 0.0 0.3 0.1 0.7 1.0 2224 1.0
## probs[7] 0.7 0.0 0.2 0.1 0.7 1.0 2251 1.0
## probs[8] 0.7 0.0 0.2 0.2 0.7 1.0 563 1.0
## probs[9] 0.7 0.0 0.2 0.3 0.8 1.0 431 1.0
## probs[10] 0.7 0.0 0.2 0.3 0.8 1.0 266 1.0
## probs[11] 0.7 0.0 0.2 0.3 0.8 1.0 145 1.0
## probs[12] 0.8 0.0 0.2 0.3 0.8 1.0 146 1.0
## probs[13] 0.7 0.0 0.2 0.3 0.8 1.0 302 1.0
## probs[14] 0.7 0.0 0.2 0.2 0.7 1.0 1116 1.0
## probs[15] 0.7 0.0 0.2 0.2 0.7 1.0 810 1.0
## probs[16] 0.7 0.0 0.2 0.3 0.7 1.0 436 1.0
## probs[17] 0.7 0.0 0.2 0.3 0.8 1.0 161 1.0
## probs[18] 0.7 0.0 0.2 0.3 0.8 1.0 151 1.0
## probs[19] 0.7 0.0 0.2 0.3 0.7 1.0 292 1.0
## probs[20] 0.7 0.0 0.2 0.2 0.7 1.0 379 1.0
## probs[21] 0.6 0.0 0.2 0.2 0.7 1.0 446 1.0
## probs[22] 0.6 0.0 0.2 0.1 0.6 1.0 2477 1.0
## probs[23] 0.5 0.0 0.2 0.1 0.5 0.9 2285 1.0
## probs[24] 0.5 0.0 0.2 0.1 0.5 0.8 702 1.0
## probs[25] 0.4 0.0 0.2 0.1 0.4 0.8 765 1.0
## probs[26] 0.4 0.0 0.2 0.0 0.4 0.8 550 1.0
## probs[27] 0.4 0.0 0.2 0.1 0.4 0.9 3080 1.0
## probs[28] 0.5 0.0 0.2 0.1 0.4 0.9 2202 1.0
## probs[29] 0.4 0.0 0.2 0.1 0.4 0.8 3117 1.0
## probs[30] 0.4 0.0 0.2 0.1 0.4 0.8 2008 1.0
## probs[31] 0.4 0.0 0.2 0.1 0.4 0.8 3586 1.0
## probs[32] 0.4 0.0 0.2 0.1 0.3 0.8 3395 1.0
## probs[33] 0.3 0.0 0.2 0.0 0.3 0.7 209 1.0
## probs[34] 0.2 0.0 0.2 0.0 0.2 0.6 107 1.0
## probs[35] 0.2 0.0 0.2 0.0 0.2 0.6 73 1.1
## probs[36] 0.2 0.0 0.2 0.0 0.2 0.6 53 1.1
## probs[37] 0.2 0.0 0.2 0.0 0.2 0.5 52 1.1
## probs[38] 0.2 0.0 0.2 0.0 0.2 0.6 54 1.1
## probs[39] 0.3 0.0 0.2 0.0 0.2 0.6 97 1.0
## probs[40] 0.3 0.0 0.2 0.0 0.3 0.7 177 1.0
## probs[41] 0.4 0.0 0.2 0.1 0.4 0.7 1985 1.0
## probs[42] 0.5 0.0 0.2 0.2 0.5 0.9 366 1.0
## probs[43] 0.6 0.0 0.2 0.2 0.6 1.0 95 1.0
## probs[44] 0.6 0.0 0.2 0.3 0.6 1.0 59 1.1
## probs[45] 0.6 0.0 0.2 0.3 0.6 1.0 56 1.1
## probs[46] 0.6 0.0 0.2 0.3 0.6 1.0 112 1.0
## probs[47] 0.5 0.0 0.2 0.2 0.5 0.9 1849 1.0
## probs[48] 0.5 0.0 0.2 0.2 0.5 0.9 854 1.0
## probs[49] 0.5 0.0 0.2 0.1 0.5 0.9 4498 1.0
## probs[50] 0.4 0.0 0.2 0.1 0.4 0.8 4135 1.0
## probs[51] 0.4 0.0 0.2 0.1 0.4 0.8 4063 1.0
## probs[52] 0.4 0.0 0.2 0.0 0.4 0.7 373 1.0
## probs[53] 0.3 0.0 0.2 0.0 0.3 0.7 144 1.0
## probs[54] 0.3 0.0 0.2 0.0 0.3 0.7 218 1.0
## probs[55] 0.4 0.0 0.2 0.1 0.4 0.7 601 1.0
## probs[56] 0.5 0.0 0.2 0.1 0.4 0.8 3550 1.0
## probs[57] 0.5 0.0 0.2 0.1 0.5 0.8 4240 1.0
## probs[58] 0.5 0.0 0.2 0.2 0.5 0.9 221 1.0
## probs[59] 0.6 0.0 0.2 0.2 0.6 1.0 113 1.0
## probs[60] 0.6 0.0 0.2 0.2 0.5 0.9 140 1.0
## probs[61] 0.5 0.0 0.2 0.1 0.5 0.9 470 1.0
## probs[62] 0.4 0.0 0.2 0.1 0.3 0.7 1235 1.0
## probs[63] 0.3 0.0 0.2 0.0 0.3 0.7 225 1.0
## probs[64] 0.3 0.0 0.2 0.0 0.2 0.6 117 1.0
## probs[65] 0.2 0.0 0.2 0.0 0.2 0.6 42 1.1
## probs[66] 0.2 0.0 0.2 0.0 0.2 0.6 43 1.1
## probs[67] 0.2 0.0 0.2 0.0 0.2 0.6 44 1.1
## probs[68] 0.3 0.0 0.2 0.0 0.2 0.6 59 1.1
## probs[69] 0.3 0.0 0.2 0.0 0.3 0.6 112 1.0
## probs[70] 0.4 0.0 0.2 0.0 0.4 0.7 443 1.0
## probs[71] 0.5 0.0 0.2 0.1 0.5 0.9 2347 1.0
## probs[72] 0.5 0.0 0.2 0.2 0.5 0.9 945 1.0
## probs[73] 0.5 0.0 0.2 0.2 0.5 0.9 3633 1.0
## probs[74] 0.6 0.0 0.2 0.2 0.6 0.9 273 1.0
## probs[75] 0.6 0.0 0.2 0.3 0.6 1.0 167 1.0
## probs[76] 0.6 0.0 0.2 0.3 0.6 1.0 264 1.0
## probs[77] 0.6 0.0 0.2 0.2 0.6 0.9 3761 1.0
## probs[78] 0.6 0.0 0.2 0.2 0.6 0.9 483 1.0
## probs[79] 0.6 0.0 0.2 0.2 0.6 0.9 422 1.0
## probs[80] 0.6 0.0 0.2 0.2 0.6 0.9 4184 1.0
## probs[81] 0.4 0.0 0.2 0.1 0.4 0.8 245 1.0
## probs[82] 0.4 0.0 0.2 0.0 0.4 0.8 99 1.0
## probs[83] 0.4 0.0 0.2 0.0 0.4 0.7 69 1.1
## probs[84] 0.4 0.0 0.2 0.0 0.4 0.8 86 1.0
## probs[85] 0.4 0.0 0.2 0.1 0.5 0.8 148 1.0
## probs[86] 0.5 0.0 0.2 0.1 0.6 0.9 4038 1.0
## probs[87] 0.6 0.0 0.2 0.1 0.6 1.0 3517 1.0
## probs[88] 0.6 0.0 0.2 0.1 0.7 1.0 3204 1.0
## probs[89] 0.7 0.0 0.2 0.1 0.7 1.0 3456 1.0
## probs[90] 0.7 0.0 0.2 0.2 0.7 1.0 3249 1.0
## probs[91] 0.7 0.0 0.2 0.3 0.8 1.0 2428 1.0
## probs[92] 0.8 0.0 0.2 0.4 0.8 1.0 645 1.0
## probs[93] 0.8 0.0 0.1 0.4 0.8 1.0 455 1.0
## probs[94] 0.8 0.0 0.1 0.5 0.8 1.0 1157 1.0
## probs[95] 0.8 0.0 0.2 0.4 0.8 1.0 1674 1.0
## probs[96] 0.8 0.0 0.1 0.5 0.9 1.0 638 1.0
## probs[97] 0.9 0.0 0.1 0.6 0.9 1.0 282 1.0
## probs[98] 0.9 0.0 0.1 0.6 0.9 1.0 178 1.0
## probs[99] 0.9 0.0 0.1 0.6 0.9 1.0 130 1.0
## probs[100] 0.9 0.0 0.1 0.6 0.9 1.0 98 1.0
## probs[101] 0.9 0.0 0.1 0.6 0.9 1.0 83 1.0
## probs[102] 0.9 0.0 0.1 0.6 0.9 1.0 96 1.0
## probs[103] 0.9 0.0 0.1 0.6 0.9 1.0 108 1.0
## probs[104] 0.9 0.0 0.1 0.6 0.9 1.0 81 1.0
## probs[105] 0.9 0.0 0.1 0.6 0.9 1.0 237 1.0
## probs[106] 0.8 0.0 0.1 0.6 0.9 1.0 303 1.0
## probs[107] 0.8 0.0 0.1 0.5 0.8 1.0 894 1.0
## probs[108] 0.8 0.0 0.2 0.4 0.8 1.0 4210 1.0
## probs[109] 0.7 0.0 0.2 0.2 0.7 0.9 209 1.0
## probs[110] 0.6 0.0 0.2 0.1 0.7 0.9 177 1.0
## probs[111] 0.7 0.0 0.2 0.2 0.7 1.0 1150 1.0
## probs[112] 0.7 0.0 0.2 0.2 0.7 1.0 1051 1.0
## probs[113] 0.7 0.0 0.2 0.1 0.7 1.0 1425 1.0
## probs[114] 0.7 0.0 0.2 0.1 0.7 1.0 1359 1.0
## probs[115] 0.7 0.0 0.2 0.1 0.7 1.0 1254 1.0
## probs[116] 0.7 0.0 0.2 0.1 0.7 1.0 1234 1.0
## probs[117] 0.7 0.0 0.2 0.1 0.7 1.0 989 1.0
## probs[118] 0.7 0.0 0.2 0.2 0.7 0.9 744 1.0
## probs[119] 0.7 0.0 0.2 0.3 0.8 1.0 3775 1.0
## probs[120] 0.8 0.0 0.2 0.4 0.8 1.0 884 1.0
## probs[121] 0.8 0.0 0.2 0.4 0.8 1.0 398 1.0
## probs[122] 0.8 0.0 0.2 0.4 0.8 1.0 413 1.0
## probs[123] 0.7 0.0 0.2 0.4 0.8 1.0 1149 1.0
## probs[124] 0.7 0.0 0.2 0.3 0.7 0.9 2718 1.0
## probs[125] 0.7 0.0 0.2 0.3 0.7 0.9 4585 1.0
## probs[126] 0.7 0.0 0.2 0.2 0.7 0.9 2597 1.0
## probs[127] 0.7 0.0 0.2 0.3 0.7 1.0 2598 1.0
## probs[128] 0.7 0.0 0.2 0.4 0.7 1.0 475 1.0
## probs[129] 0.7 0.0 0.2 0.4 0.7 1.0 382 1.0
## probs[130] 0.7 0.0 0.2 0.3 0.7 0.9 3316 1.0
## probs[131] 0.6 0.0 0.2 0.2 0.6 0.9 896 1.0
## probs[132] 0.5 0.0 0.2 0.1 0.6 0.9 820 1.0
## probs[133] 0.6 0.0 0.2 0.2 0.6 0.9 4783 1.0
## probs[134] 0.6 0.0 0.2 0.2 0.6 0.9 4443 1.0
## probs[135] 0.5 0.0 0.2 0.1 0.5 0.9 2117 1.0
## probs[136] 0.5 0.0 0.2 0.1 0.5 0.9 4450 1.0
## probs[137] 0.5 0.0 0.2 0.1 0.5 0.8 362 1.0
## probs[138] 0.5 0.0 0.2 0.1 0.5 0.8 256 1.0
## probs[139] 0.5 0.0 0.2 0.1 0.5 0.8 947 1.0
## probs[140] 0.6 0.0 0.2 0.2 0.6 0.9 1150 1.0
## probs[141] 0.7 0.0 0.2 0.3 0.7 1.0 148 1.0
## probs[142] 0.7 0.0 0.2 0.3 0.7 1.0 103 1.0
## probs[143] 0.7 0.0 0.2 0.3 0.7 1.0 67 1.1
## probs[144] 0.6 0.0 0.2 0.3 0.7 1.0 81 1.0
## probs[145] 0.6 0.0 0.2 0.2 0.6 0.9 232 1.0
## probs[146] 0.5 0.0 0.2 0.1 0.5 0.9 3993 1.0
## probs[147] 0.4 0.0 0.2 0.1 0.4 0.8 4049 1.0
## probs[148] 0.3 0.0 0.2 0.0 0.3 0.7 477 1.0
## probs[149] 0.2 0.0 0.2 0.0 0.2 0.6 128 1.0
## probs[150] 0.2 0.0 0.1 0.0 0.2 0.5 82 1.1
## probs[151] 0.2 0.0 0.1 0.0 0.2 0.5 59 1.1
## probs[152] 0.2 0.0 0.1 0.0 0.1 0.5 57 1.1
## probs[153] 0.2 0.0 0.1 0.0 0.1 0.5 51 1.1
## probs[154] 0.2 0.0 0.1 0.0 0.1 0.5 49 1.1
## probs[155] 0.2 0.0 0.1 0.0 0.1 0.5 59 1.1
## probs[156] 0.2 0.0 0.1 0.0 0.2 0.5 86 1.0
## probs[157] 0.2 0.0 0.1 0.0 0.2 0.5 156 1.0
## probs[158] 0.2 0.0 0.2 0.0 0.2 0.6 636 1.0
## probs[159] 0.2 0.0 0.2 0.0 0.2 0.6 142 1.0
## probs[160] 0.2 0.0 0.1 0.0 0.2 0.6 95 1.0
## probs[161] 0.2 0.0 0.2 0.0 0.2 0.6 102 1.0
## probs[162] 0.3 0.0 0.2 0.0 0.2 0.6 157 1.0
## probs[163] 0.3 0.0 0.2 0.0 0.3 0.7 291 1.0
## probs[164] 0.4 0.0 0.2 0.1 0.4 0.8 2290 1.0
## probs[165] 0.5 0.0 0.2 0.2 0.5 0.9 299 1.0
## probs[166] 0.6 0.0 0.2 0.3 0.6 1.0 92 1.0
## probs[167] 0.7 0.0 0.2 0.3 0.7 1.0 53 1.1
## probs[168] 0.7 0.0 0.2 0.3 0.7 1.0 42 1.1
## probs[169] 0.7 0.0 0.2 0.3 0.7 1.0 36 1.1
## probs[170] 0.7 0.0 0.2 0.3 0.7 1.0 42 1.1
## probs[171] 0.6 0.0 0.2 0.3 0.6 1.0 101 1.0
## probs[172] 0.5 0.0 0.2 0.2 0.5 0.9 3301 1.0
## probs[173] 0.5 0.0 0.2 0.2 0.5 0.9 3358 1.0
## probs[174] 0.4 0.0 0.2 0.1 0.4 0.8 1362 1.0
## probs[175] 0.4 0.0 0.2 0.1 0.4 0.8 1215 1.0
## probs[176] 0.4 0.0 0.2 0.1 0.4 0.8 4390 1.0
## probs[177] 0.4 0.0 0.2 0.1 0.4 0.8 913 1.0
## probs[178] 0.4 0.0 0.2 0.1 0.4 0.7 806 1.0
## probs[179] 0.4 0.0 0.2 0.1 0.4 0.8 4000 1.0
## probs[180] 0.3 0.0 0.2 0.0 0.3 0.7 682 1.0
## probs[181] 0.3 0.0 0.2 0.0 0.3 0.7 833 1.0
## probs[182] 0.4 0.0 0.2 0.1 0.4 0.8 4203 1.0
## probs[183] 0.4 0.0 0.2 0.0 0.3 0.8 2139 1.0
## lp__ -90.0 24.3 92.4 -259.8 -94.7 109.1 14 1.2
##
## Samples were drawn using NUTS(diag_e) at Sat Jan 20 16:22:47 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).単純に勝率x%とせずに状態空間モデルを利用すれば、勝率の推移を計算できる、時点によっての各チームの勝率を推定できる。
- 1829 - 1849: ケンブリッジが優勢
- 1850 - 1914: オックスフォードが優勢もケンブリッジが勝つときもある
- 1915 - 1970: 再びケンブリッジが優勢
- 1971 - 1995: オックスフォードが優勢
years <- seq(from = as.POSIXct("1829-01-01"), by = "1 year", len = length(boat))
ms <- rstan::extract(fit)
d_est <- data.frame(t(apply(ms$probs, 2, quantile, probs = c(0.025, 0.5, 0.975))))
colnames(d_est) <- c("lwr", "fit", "upr")
d_est <- cbind(data.frame(years, game = as.numeric(boat)), d_est)
ggplot(data = d_est, aes(x = years)) +
theme_bw(base_size = 15) +
geom_point(aes(x = years, y = game), alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit)) +
geom_ribbon(aes(x = years, y = fit, ymin = lwr, ymax = upr), alpha = 0.3) +
labs(x = 'years', y = 'win rate', title = 'Cambridge University\'s Winning Percentage') +
scale_x_datetime(breaks = scales::date_breaks("10 year"), date_labels = "%Y") +
theme_bw()