UPDATE: 2024-01-23 16:05:44.967038
このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回はRasmus Bååth先生のブログに記載されているサッカーチームの強さを推定する内容を参考にさせていただき、写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
こちらのブログではJAGSを利用しており、再現はStanで行うため、結果は完全に一致しない。また、JAGSとStanでは正規分布のパラメタの与え方が異なるので、そのあたりは修正を行っているが、ブログの数式はそのまま表示しているので、数式とStanモデルでは差異がある点は注意が必要である。
この分析は、
これが本分析の目的とのこと。まずは必要なデータやパッケージを準備しておく。ブログで使用されている通り、使用するデータはラ・リーガの5つのシーズン(2008/09, 2009/10, 2010/11, 2011/12, 2012/13
)の試合結果が記録されているデータを利用する。このデータはseR
2013データ分析コンテストの一環として提供されたリーガ・エスパニョーラのデータセット。
library(tidyverse)
library(rstan)
library(patchwork)
library(xtable)
source("plotPost.R")
options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
set.seed(12345)
# Convenience function to generate the type of column names Jags outputs.
col_name <- function(name, ...) {
paste0(name, "[", paste(..., sep=",") , "]")
}
load("laliga.RData")
head(laliga)
## Season Week HomeTeam AwayTeam HomeGoals
## 1 2008/09 1 Athletic Club Bilbao Union Deportiva Almeria 1
## 2 2008/09 1 Atlético Madrid Málaga CF 4
## 3 2008/09 1 Betis Sevilla Real Club Recreativo Huelva 0
## 4 2008/09 1 CA Osasuna Villarreal CF 1
## 5 2008/09 1 CD Numancia FC Barcelona 1
## 6 2008/09 1 Deportivo de La Coruña Real Madrid CF 2
## AwayGoals
## 1 3
## 2 0
## 3 1
## 4 1
## 5 0
## 6 1
モデリングや可視化で必要になるカラムを作成しておく。2012/13
シーズンのいくつかのレコードには欠損値が含まれているので、そのレコードは除外しておく。
# -1 Away win, 0 Draw, 1 Home win
laliga$MatchResult <- sign(laliga$HomeGoals - laliga$AwayGoals)
# Creating a data frame d with only the complete match results
# laliga %>% filter(is.na(HomeGoals))
d <- na.omit(laliga)
ここではいくつかの特徴(昇格、降格)があるチームをサンプリングしてモデリングを行う。すべてのシーズン、すべてのチームのデータを利用して推定したほうが良いと思うが、あくまでも学習用なので、ここでは時間を考慮して必要なチームのレコードを利用する。
# Sampling
t <- c(
'Real Valladolid',
'Atlético Madrid',
'FC Barcelona',
'Villarreal CF',
'FC Sevilla',
'Real Madrid CF',
'FC Valencia',
'Real Zaragoza',
'CD Tenerife'
)
d <- d %>%
filter(HomeTeam %in% t & AwayTeam %in% t)
dim(d)
## [1] 235 7
特徴に関して、下記の通り説明を記載しておく。
Villarreal CF
:
2012/13
シーズンは1部リーグではないReal Zaragoza
:
2008/09
シーズンは1部リーグではないReal Valladolid
:
2010/11,2011/12
シーズンは1部リーグではないCD Tenerife
:
2009/10
シーズンのみ1部リーグplt_particicipation <- d %>%
distinct(Season, HomeTeam) %>%
ggplot(., aes(Season, HomeTeam)) +
theme_bw(base_size = 15) +
geom_point(size = 5) +
labs(xlab = 'Season', y = 'Team', title = 'Particicipation by Season')
plt_particicipation
モデリングに関しては3段階に分けて行い、徐々に複雑にしていく。最終的なモデルでもそのままデータが利用できるように、ここでStanに渡すデータを作成する。
teams <- unique(c(d$HomeTeam, d$AwayTeam))
seasons <- unique(d$Season)
data_list <- list(
HomeGoals = d$HomeGoals,
AwayGoals = d$AwayGoals,
HomeTeam = as.numeric(factor(d$HomeTeam, levels = teams)),
AwayTeam = as.numeric(factor(d$AwayTeam, levels = teams)),
Season = as.numeric(factor(d$Season, levels=seasons)),
n_teams = length(teams),
n_games = nrow(d),
n_seasons = length(seasons)
)
data_list
## $HomeGoals
## [1] 2 0 6 0 1 4 3 6 1 0 0 0 4 3 2 1 1 1 3 1 3 1 1 4 1 1 3 4 0 2 3 0 4 2 2 2 0
## [38] 3 3 1 3 3 2 1 4 2 5 1 0 2 3 2 2 0 4 6 5 2 3 3 1 3 1 3 1 1 2 1 6 2 1 0 0 4
## [75] 4 0 0 4 1 2 2 2 2 1 6 4 1 3 1 3 1 2 2 3 3 3 3 0 0 0 2 2 2 3 1 3 1 3 4 2 2
## [112] 4 1 3 1 1 1 3 2 0 2 5 1 2 2 3 1 1 5 0 2 1 1 1 4 3 1 0 2 0 3 1 1 4 1 0 3 0
## [149] 5 1 1 3 3 2 2 1 5 0 1 2 2 5 1 0 2 0 3 0 3 4 2 4 0 1 2 2 3 3 0 0 1 5 1 1 1
## [186] 2 0 1 1 1 0 1 1 3 1 3 1 0 1 0 1 1 2 2 2 2 2 2 4 1 3 4 1 2 2 4 1 2 2 0 2 1
## [223] 4 0 2 2 2 4 1 1 2 0 0 1 1
##
## $AwayGoals
## [1] 1 1 1 1 2 4 2 0 0 0 3 3 0 4 0 0 0 2 1 0 3 2 0 3 2 1 2 1 1 0 1 0 0 4 2 6 2
## [38] 3 0 0 2 1 0 0 1 4 2 2 2 2 0 1 1 0 2 1 0 3 2 1 2 1 0 3 1 2 3 1 0 1 1 5 4 0
## [75] 1 0 3 2 3 1 0 1 1 1 2 1 5 2 1 0 4 4 2 2 0 0 0 2 0 2 1 0 1 1 2 1 4 1 1 3 0
## [112] 0 0 3 2 1 0 1 1 2 0 0 1 0 0 1 2 1 0 3 0 0 3 0 2 0 2 1 2 1 1 0 1 0 2 1 1 1
## [149] 0 1 0 2 6 3 6 3 0 6 0 2 2 0 0 0 2 0 0 1 1 0 3 1 1 3 6 2 0 1 0 0 2 1 2 2 1
## [186] 1 2 2 1 0 0 4 4 0 2 0 0 1 1 1 0 0 1 3 0 2 1 0 0 1 1 0 2 0 3 1 3 0 0 5 0 1
## [223] 1 3 1 2 1 0 1 1 1 3 1 2 1
##
## $HomeTeam
## [1] 1 2 3 1 2 4 1 3 1 5 4 5 3 6 3 5 6 4 7 6 7 2 5 2 7 6 2 5 1 6 7 1 3 5 7 6 4
## [38] 3 7 2 4 4 7 8 5 1 3 8 4 7 6 2 5 7 6 3 4 2 5 7 9 4 3 1 5 2 7 9 6 2 3 9 1 3
## [75] 7 9 1 4 9 5 7 8 2 1 6 2 9 6 8 3 1 8 9 6 8 4 5 1 9 6 1 6 4 2 8 5 4 2 3 5 4
## [112] 3 7 8 2 7 2 5 3 8 4 3 7 5 6 3 8 4 3 8 6 4 8 6 6 3 2 8 2 7 2 3 5 8 2 7 5 4
## [149] 7 6 4 5 7 6 5 4 3 8 7 4 7 3 5 2 4 3 6 8 2 3 7 6 8 6 5 4 2 6 4 2 5 3 2 7 5
## [186] 8 5 7 4 8 6 8 2 5 3 6 7 4 6 8 3 5 2 5 7 3 8 7 6 1 3 2 5 6 1 3 1 2 7 7 1 7
## [223] 6 1 3 8 6 5 2 8 7 8 5 2 1
##
## $AwayTeam
## [1] 2 5 2 7 6 2 5 1 6 7 1 3 7 5 6 4 7 3 2 4 4 1 2 3 1 2 4 1 3 1 5 4 5 6 3 3 5
## [38] 4 6 7 6 7 5 9 8 7 2 1 6 2 9 8 6 3 1 8 9 6 4 8 5 1 6 9 1 4 6 2 8 5 4 3 2 5
## [75] 4 7 3 8 8 7 1 5 3 8 4 7 6 5 2 7 6 3 4 2 7 5 9 4 1 3 5 7 2 9 6 2 3 1 9 3 7
## [112] 1 9 4 3 2 8 2 7 3 2 5 8 7 2 4 5 7 6 4 7 5 6 5 4 2 7 2 5 3 4 8 3 7 6 5 8 3
## [149] 4 3 8 4 6 8 6 6 4 6 2 5 3 2 7 5 8 5 4 7 8 8 6 2 5 3 6 7 4 8 3 7 4 7 3 5 2
## [186] 4 3 8 6 2 7 3 6 8 6 5 4 2 7 1 7 6 1 3 8 6 5 2 8 7 8 5 1 2 6 2 3 8 5 6 8 3
## [223] 5 2 5 7 3 8 7 6 1 3 2 6 5
##
## $Season
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [149] 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [186] 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [223] 5 5 5 5 5 5 5 5 5 5 5 5 5
##
## $n_teams
## [1] 9
##
## $n_games
## [1] 235
##
## $n_seasons
## [1] 5
数字とチーム名の対応は下記の通りである。
data.frame(
no = as.numeric(factor(unique(d$HomeTeam), levels=teams)),
teams
)
## no teams
## 1 1 Real Valladolid
## 2 2 Atlético Madrid
## 3 3 FC Barcelona
## 4 4 Villarreal CF
## 5 5 FC Sevilla
## 6 6 Real Madrid CF
## 7 7 FC Valencia
## 8 8 Real Zaragoza
## 9 9 CD Tenerife
サッカーの試合はポアソン分布に従うと一般的に知られている。サッカーの試合はアディショナルタイムがあれど90分という単位時間があり、そこで得点という事象が発生する回数を表す離散確率分布がポアソン分布。今回のデータも例に漏れず、ポアソン分布に従っていると仮定しても問題はなさそうである。青棒が観測値、赤線がホームゴールの平均値を\(\lambda\)とするポアソン分布。可視化はしていないが、ホームゴールだけでなく、アウェイゴールもポアソン分布に従っている。
df_homegoals <- d %>%
group_by(HomeGoals) %>%
count() %>%
ungroup() %>%
mutate(
freq = n/sum(n),
prob = dpois(HomeGoals, mean(d$HomeGoals))
)
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = df_homegoals, aes(HomeGoals, freq), stat = 'identity', fill = 'lightsteelblue', alpha = 1/2) +
geom_line(data = df_homegoals, aes(HomeGoals, prob), stat = 'identity', col = 'coral') +
scale_x_continuous(breaks = 0:10) +
labs(title = 'HomeGoals')
ブログで紹介されている1つ目のモデルを利用する。全てのチームが同じように優れているわけではなく、全てのチームが潜在的なスキル(=強さ)を持っていると仮定する。ここで、あるチームのスキルから相手チームのスキルを引いたものが試合の予測となるとする。ゴール数はポアソン分布と仮定できる。チーム\(i\)がチーム\(j\)と対戦したときののゴール数の分布は次のようになる。ここでベースラインは、両チームの実力が等しい場合の平均ゴール数である。
\[ \begin{eqnarray} Goals &\sim& \text{Poisson}(\lambda) \\ \log(\lambda) &=& \text{baseline} + \text{skill}_i - \text{skill}_j \\ \end{eqnarray} \]
ホームチーム\(i\)とアウェイチーム\(j\)の試合のゴール結果は次のようにモデル化される。
\[ \begin{eqnarray} HomeGoals_{i,j} &\sim& \text{Poison}(\lambda_{\text{home},i,j}) \\ AwayGoals_{i,j} &\sim& \text{Poison}(\lambda_{\text{away},i,j}) \\ \log(\lambda_{\text{home},i,j}) &=& \text{baseline} + \text{skill}_i - \text{skill}_j \\ \log(\lambda_{\text{away},i,j}) &=& \text{baseline} + \text{skill}_j - \text{skill}_i \end{eqnarray} \] これにいくつかの事前分布を加える。ベースラインと全\(n\)チームのスキルに関する事前分布は下記の通り。
\[ \begin{eqnarray} \text{baseline} &\sim& \text{Normal}(0, 4^2) \\ \text{skill}_{1 \ldots n} &\sim& \text{Normal}(\mu_\text{teams}, \sigma_{\text{teams}}^2) \\ \mu_\text{teams} &\sim& \text{Normal}(0, 4^2) \\ \sigma_\text{teams} &\sim& \text{Uniform}(0, 3) \\ \end{eqnarray} \] 今回ここで推定する際は、\(\mu_{teams}\)は利用しておらず、また、JAGSは分散の代わりに精度(分散の逆数)で正規分布をパラメタ化するが、ここでは、Stanを利用していることもあり、平均0、標準偏差\(\sigma_{teams}\)の正規分布を代わりに利用している。さらに、ブログでは1つのチームのスキルを定数0に固定しているが、ここでは階層モデルのようにチームの強さは特定の正規分布から生成されるゆるい制約を仮定し、変更している。
parameters {
real baseline; // ベースラインのパラメータ
real<lower=0> group_sigma; // グループの標準偏差
// real group_skill; // グループのスキル
real skill[n_teams]; // 各チームのスキル
}
transformed parameters {
matrix[n_teams, n_teams] lambda_home; // ホームチームの得点率行列
matrix[n_teams, n_teams] lambda_away; // アウェイチームの得点率行列
for (home_i in 1:n_teams) {
for (away_i in 1:n_teams) {
lambda_home[home_i, away_i] = exp(baseline + skill[home_i] - skill[away_i]);
lambda_away[home_i, away_i] = exp(baseline + skill[away_i] - skill[home_i]);
}
}
}
model {
// group_skill ~ normal(0, 4);
group_sigma ~ uniform(0, 3);
baseline ~ normal(0, 4);
for (n in 1:n_teams) {
skill[n] ~ normal(0, group_sigma);
}
for (i in 1:n_games) {
HomeGoals[i] ~ poisson(lambda_home[HomeTeam[i], AwayTeam[i]]);
AwayGoals[i] ~ poisson(lambda_away[HomeTeam[i], AwayTeam[i]]);
}
}
model1 <- stan_model('model1.stan')
sampling()
関数でサンプリングする。
fit1 <- sampling(object = model1, data = data_list, seed = 1989, chains = 4, iter = 10000, thin = 2)
推定結果を確認する。
print(fit1, prob = c(0.025, 0.5, 0.975), digits_summary = 2)
## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=5000; thin=2;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## baseline 0.41 0.00 0.04 0.33 0.41 0.48 8643 1
## group_sigma 0.39 0.00 0.13 0.22 0.36 0.72 6887 1
## skill[1] -0.19 0.00 0.16 -0.51 -0.19 0.12 4172 1
## skill[2] -0.01 0.00 0.15 -0.31 -0.01 0.29 4005 1
## skill[3] 0.59 0.00 0.15 0.29 0.58 0.90 4023 1
## skill[4] -0.05 0.00 0.15 -0.35 -0.05 0.25 4051 1
## skill[5] -0.02 0.00 0.15 -0.32 -0.01 0.28 3995 1
## skill[6] 0.38 0.00 0.15 0.09 0.38 0.68 4034 1
## skill[7] 0.00 0.00 0.15 -0.30 0.00 0.30 3996 1
## skill[8] -0.26 0.00 0.15 -0.57 -0.26 0.04 4136 1
## skill[9] -0.45 0.00 0.18 -0.81 -0.45 -0.11 4809 1
## lambda_home[1,1] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[1,2] 1.26 0.00 0.15 0.99 1.25 1.57 7729 1
## lambda_home[1,3] 0.69 0.00 0.09 0.54 0.69 0.88 8333 1
## lambda_home[1,4] 1.31 0.00 0.16 1.03 1.30 1.64 8295 1
## lambda_home[1,5] 1.27 0.00 0.15 1.00 1.26 1.57 8277 1
## lambda_home[1,6] 0.85 0.00 0.10 0.66 0.85 1.08 8165 1
## lambda_home[1,7] 1.25 0.00 0.14 0.98 1.24 1.55 7796 1
## lambda_home[1,8] 1.62 0.00 0.19 1.28 1.61 2.02 8039 1
## lambda_home[1,9] 1.96 0.00 0.30 1.44 1.94 2.60 8307 1
## lambda_home[2,1] 1.82 0.00 0.21 1.44 1.81 2.25 8212 1
## lambda_home[2,2] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[2,3] 0.83 0.00 0.09 0.67 0.83 1.02 8549 1
## lambda_home[2,4] 1.57 0.00 0.16 1.27 1.56 1.92 8387 1
## lambda_home[2,5] 1.52 0.00 0.15 1.24 1.51 1.84 8523 1
## lambda_home[2,6] 1.02 0.00 0.11 0.83 1.02 1.25 8647 1
## lambda_home[2,7] 1.49 0.00 0.15 1.21 1.49 1.81 8457 1
## lambda_home[2,8] 1.94 0.00 0.20 1.58 1.93 2.36 8532 1
## lambda_home[2,9] 2.36 0.00 0.34 1.76 2.33 3.09 8113 1
## lambda_home[3,1] 3.29 0.00 0.35 2.65 3.28 4.03 8234 1
## lambda_home[3,2] 2.74 0.00 0.25 2.27 2.73 3.25 8629 1
## lambda_home[3,3] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[3,4] 2.84 0.00 0.27 2.35 2.83 3.40 8182 1
## lambda_home[3,5] 2.76 0.00 0.25 2.30 2.74 3.26 8575 1
## lambda_home[3,6] 1.86 0.00 0.18 1.53 1.85 2.22 8763 1
## lambda_home[3,7] 2.71 0.00 0.24 2.25 2.70 3.22 8565 1
## lambda_home[3,8] 3.52 0.00 0.33 2.92 3.51 4.19 8486 1
## lambda_home[3,9] 4.27 0.01 0.60 3.21 4.23 5.56 8294 1
## lambda_home[4,1] 1.75 0.00 0.21 1.38 1.74 2.17 8129 1
## lambda_home[4,2] 1.45 0.00 0.15 1.17 1.45 1.78 7816 1
## lambda_home[4,3] 0.80 0.00 0.09 0.64 0.80 0.99 7939 1
## lambda_home[4,4] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[4,5] 1.46 0.00 0.15 1.19 1.46 1.79 8321 1
## lambda_home[4,6] 0.99 0.00 0.11 0.78 0.98 1.22 8175 1
## lambda_home[4,7] 1.44 0.00 0.15 1.16 1.43 1.75 7888 1
## lambda_home[4,8] 1.87 0.00 0.20 1.51 1.86 2.29 8061 1
## lambda_home[4,9] 2.27 0.00 0.33 1.68 2.25 2.99 7982 1
## lambda_home[5,1] 1.80 0.00 0.20 1.44 1.79 2.22 8037 1
## lambda_home[5,2] 1.50 0.00 0.15 1.22 1.49 1.81 8064 1
## lambda_home[5,3] 0.83 0.00 0.09 0.67 0.82 1.01 8229 1
## lambda_home[5,4] 1.56 0.00 0.16 1.26 1.55 1.89 8353 1
## lambda_home[5,5] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[5,6] 1.02 0.00 0.11 0.82 1.01 1.24 8234 1
## lambda_home[5,7] 1.48 0.00 0.15 1.21 1.48 1.79 8169 1
## lambda_home[5,8] 1.93 0.00 0.19 1.57 1.92 2.33 8109 1
## lambda_home[5,9] 2.34 0.00 0.34 1.74 2.32 3.07 8035 1
## lambda_home[6,1] 2.68 0.00 0.29 2.15 2.66 3.28 7998 1
## lambda_home[6,2] 2.22 0.00 0.21 1.84 2.22 2.66 8540 1
## lambda_home[6,3] 1.23 0.00 0.12 1.00 1.22 1.48 8234 1
## lambda_home[6,4] 2.31 0.00 0.23 1.89 2.30 2.79 8186 1
## lambda_home[6,5] 2.24 0.00 0.21 1.84 2.23 2.68 8406 1
## lambda_home[6,6] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[6,7] 2.20 0.00 0.21 1.82 2.19 2.62 8292 1
## lambda_home[6,8] 2.86 0.00 0.28 2.35 2.86 3.43 8548 1
## lambda_home[6,9] 3.47 0.01 0.49 2.59 3.44 4.52 8030 1
## lambda_home[7,1] 1.83 0.00 0.20 1.46 1.82 2.27 7837 1
## lambda_home[7,2] 1.52 0.00 0.15 1.24 1.52 1.84 8283 1
## lambda_home[7,3] 0.84 0.00 0.09 0.68 0.84 1.03 8577 1
## lambda_home[7,4] 1.58 0.00 0.16 1.28 1.58 1.93 8119 1
## lambda_home[7,5] 1.54 0.00 0.15 1.25 1.53 1.86 8376 1
## lambda_home[7,6] 1.03 0.00 0.11 0.83 1.03 1.26 8579 1
## lambda_home[7,7] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[7,8] 1.96 0.00 0.20 1.60 1.96 2.38 8101 1
## lambda_home[7,9] 2.38 0.00 0.34 1.77 2.35 3.11 7902 1
## lambda_home[8,1] 1.41 0.00 0.17 1.11 1.40 1.76 8007 1
## lambda_home[8,2] 1.17 0.00 0.13 0.94 1.17 1.43 8315 1
## lambda_home[8,3] 0.65 0.00 0.08 0.51 0.64 0.80 8257 1
## lambda_home[8,4] 1.22 0.00 0.14 0.97 1.21 1.50 8233 1
## lambda_home[8,5] 1.18 0.00 0.13 0.95 1.18 1.45 8243 1
## lambda_home[8,6] 0.80 0.00 0.09 0.63 0.79 0.99 8469 1
## lambda_home[8,7] 1.16 0.00 0.13 0.93 1.16 1.42 7986 1
## lambda_home[8,8] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_home[8,9] 1.83 0.00 0.27 1.36 1.81 2.41 8459 1
## lambda_home[9,1] 1.18 0.00 0.19 0.86 1.16 1.57 8935 1
## lambda_home[9,2] 0.98 0.00 0.15 0.71 0.97 1.30 8561 1
## lambda_home[9,3] 0.54 0.00 0.09 0.39 0.53 0.73 8718 1
## lambda_home[9,4] 1.02 0.00 0.16 0.74 1.00 1.36 8989 1
## lambda_home[9,5] 0.99 0.00 0.15 0.73 0.97 1.31 8821 1
## lambda_home[9,6] 0.66 0.00 0.11 0.48 0.66 0.90 8730 1
## lambda_home[9,7] 0.97 0.00 0.15 0.71 0.96 1.28 8477 1
## lambda_home[9,8] 1.26 0.00 0.19 0.93 1.25 1.67 8943 1
## lambda_home[9,9] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[1,1] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[1,2] 1.82 0.00 0.21 1.44 1.81 2.25 8212 1
## lambda_away[1,3] 3.29 0.00 0.35 2.65 3.28 4.03 8234 1
## lambda_away[1,4] 1.75 0.00 0.21 1.38 1.74 2.17 8129 1
## lambda_away[1,5] 1.80 0.00 0.20 1.44 1.79 2.22 8037 1
## lambda_away[1,6] 2.68 0.00 0.29 2.15 2.66 3.28 7998 1
## lambda_away[1,7] 1.83 0.00 0.20 1.46 1.82 2.27 7837 1
## lambda_away[1,8] 1.41 0.00 0.17 1.11 1.40 1.76 8007 1
## lambda_away[1,9] 1.18 0.00 0.19 0.86 1.16 1.57 8935 1
## lambda_away[2,1] 1.26 0.00 0.15 0.99 1.25 1.57 7729 1
## lambda_away[2,2] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[2,3] 2.74 0.00 0.25 2.27 2.73 3.25 8629 1
## lambda_away[2,4] 1.45 0.00 0.15 1.17 1.45 1.78 7816 1
## lambda_away[2,5] 1.50 0.00 0.15 1.22 1.49 1.81 8064 1
## lambda_away[2,6] 2.22 0.00 0.21 1.84 2.22 2.66 8540 1
## lambda_away[2,7] 1.52 0.00 0.15 1.24 1.52 1.84 8283 1
## lambda_away[2,8] 1.17 0.00 0.13 0.94 1.17 1.43 8315 1
## lambda_away[2,9] 0.98 0.00 0.15 0.71 0.97 1.30 8561 1
## lambda_away[3,1] 0.69 0.00 0.09 0.54 0.69 0.88 8333 1
## lambda_away[3,2] 0.83 0.00 0.09 0.67 0.83 1.02 8549 1
## lambda_away[3,3] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[3,4] 0.80 0.00 0.09 0.64 0.80 0.99 7939 1
## lambda_away[3,5] 0.83 0.00 0.09 0.67 0.82 1.01 8229 1
## lambda_away[3,6] 1.23 0.00 0.12 1.00 1.22 1.48 8234 1
## lambda_away[3,7] 0.84 0.00 0.09 0.68 0.84 1.03 8577 1
## lambda_away[3,8] 0.65 0.00 0.08 0.51 0.64 0.80 8257 1
## lambda_away[3,9] 0.54 0.00 0.09 0.39 0.53 0.73 8718 1
## lambda_away[4,1] 1.31 0.00 0.16 1.03 1.30 1.64 8295 1
## lambda_away[4,2] 1.57 0.00 0.16 1.27 1.56 1.92 8387 1
## lambda_away[4,3] 2.84 0.00 0.27 2.35 2.83 3.40 8182 1
## lambda_away[4,4] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[4,5] 1.56 0.00 0.16 1.26 1.55 1.89 8353 1
## lambda_away[4,6] 2.31 0.00 0.23 1.89 2.30 2.79 8186 1
## lambda_away[4,7] 1.58 0.00 0.16 1.28 1.58 1.93 8119 1
## lambda_away[4,8] 1.22 0.00 0.14 0.97 1.21 1.50 8233 1
## lambda_away[4,9] 1.02 0.00 0.16 0.74 1.00 1.36 8989 1
## lambda_away[5,1] 1.27 0.00 0.15 1.00 1.26 1.57 8277 1
## lambda_away[5,2] 1.52 0.00 0.15 1.24 1.51 1.84 8523 1
## lambda_away[5,3] 2.76 0.00 0.25 2.30 2.74 3.26 8575 1
## lambda_away[5,4] 1.46 0.00 0.15 1.19 1.46 1.79 8321 1
## lambda_away[5,5] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[5,6] 2.24 0.00 0.21 1.84 2.23 2.68 8406 1
## lambda_away[5,7] 1.54 0.00 0.15 1.25 1.53 1.86 8376 1
## lambda_away[5,8] 1.18 0.00 0.13 0.95 1.18 1.45 8243 1
## lambda_away[5,9] 0.99 0.00 0.15 0.73 0.97 1.31 8821 1
## lambda_away[6,1] 0.85 0.00 0.10 0.66 0.85 1.08 8165 1
## lambda_away[6,2] 1.02 0.00 0.11 0.83 1.02 1.25 8647 1
## lambda_away[6,3] 1.86 0.00 0.18 1.53 1.85 2.22 8763 1
## lambda_away[6,4] 0.99 0.00 0.11 0.78 0.98 1.22 8175 1
## lambda_away[6,5] 1.02 0.00 0.11 0.82 1.01 1.24 8234 1
## lambda_away[6,6] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[6,7] 1.03 0.00 0.11 0.83 1.03 1.26 8579 1
## lambda_away[6,8] 0.80 0.00 0.09 0.63 0.79 0.99 8469 1
## lambda_away[6,9] 0.66 0.00 0.11 0.48 0.66 0.90 8730 1
## lambda_away[7,1] 1.25 0.00 0.14 0.98 1.24 1.55 7796 1
## lambda_away[7,2] 1.49 0.00 0.15 1.21 1.49 1.81 8457 1
## lambda_away[7,3] 2.71 0.00 0.24 2.25 2.70 3.22 8565 1
## lambda_away[7,4] 1.44 0.00 0.15 1.16 1.43 1.75 7888 1
## lambda_away[7,5] 1.48 0.00 0.15 1.21 1.48 1.79 8169 1
## lambda_away[7,6] 2.20 0.00 0.21 1.82 2.19 2.62 8292 1
## lambda_away[7,7] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[7,8] 1.16 0.00 0.13 0.93 1.16 1.42 7986 1
## lambda_away[7,9] 0.97 0.00 0.15 0.71 0.96 1.28 8477 1
## lambda_away[8,1] 1.62 0.00 0.19 1.28 1.61 2.02 8039 1
## lambda_away[8,2] 1.94 0.00 0.20 1.58 1.93 2.36 8532 1
## lambda_away[8,3] 3.52 0.00 0.33 2.92 3.51 4.19 8486 1
## lambda_away[8,4] 1.87 0.00 0.20 1.51 1.86 2.29 8061 1
## lambda_away[8,5] 1.93 0.00 0.19 1.57 1.92 2.33 8109 1
## lambda_away[8,6] 2.86 0.00 0.28 2.35 2.86 3.43 8548 1
## lambda_away[8,7] 1.96 0.00 0.20 1.60 1.96 2.38 8101 1
## lambda_away[8,8] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lambda_away[8,9] 1.26 0.00 0.19 0.93 1.25 1.67 8943 1
## lambda_away[9,1] 1.96 0.00 0.30 1.44 1.94 2.60 8307 1
## lambda_away[9,2] 2.36 0.00 0.34 1.76 2.33 3.09 8113 1
## lambda_away[9,3] 4.27 0.01 0.60 3.21 4.23 5.56 8294 1
## lambda_away[9,4] 2.27 0.00 0.33 1.68 2.25 2.99 7982 1
## lambda_away[9,5] 2.34 0.00 0.34 1.74 2.32 3.07 8035 1
## lambda_away[9,6] 3.47 0.01 0.49 2.59 3.44 4.52 8030 1
## lambda_away[9,7] 2.38 0.00 0.34 1.77 2.35 3.11 7902 1
## lambda_away[9,8] 1.83 0.00 0.27 1.36 1.81 2.41 8459 1
## lambda_away[9,9] 1.50 0.00 0.06 1.39 1.50 1.62 8646 1
## lp__ -315.57 0.03 2.40 -321.02 -315.26 -311.84 5760 1
##
## Samples were drawn using NUTS(diag_e) at Tue Jan 23 16:05:50 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).
MCMCサンプルを使って、FCセビージャとFCバレンシアのトレース・プロットとスキル・パラメータの分布を可視化しておく。
plt_a <- stan_trace(fit1,
pars = col_name("skill", which(teams == "FC Sevilla")),
separate_chains = TRUE) + ggtitle('Traceplot of FC Sevilla')
plt_b <- stan_dens(fit1,
pars = col_name("skill", which(teams == "FC Sevilla")),
separate_chains = TRUE)
plt_c <- stan_trace(fit1,
pars = col_name("skill", which(teams == "FC Valencia")),
separate_chains = TRUE) + ggtitle('Traceplot of FC Valencia')
plt_d <- stan_dens(fit1,
pars = col_name("skill", which(teams == "FC Valencia")),
separate_chains = TRUE)
(plt_a | plt_b) / (plt_c | plt_d)
セビージャとバレンシアは似たようなスキルを持っているものの、バレンシアの方がわずかに強いと推定されている。MCMCサンプルを使ってチーム間の試合をシミュレートして、ゴール数の分布や、ホームチームの勝利、アウェイチームの勝利、引き分けの確率を見ることもできる。
# ホームチーム名、アウェイチーム名、MCMCサンプルが格納されたmsオブジェクトを受け取る
# MCMCサンプルからモデルのパラメータを取り出し、ホームチームとアウェイチームのスキルを取り出し
# ポアソン分布からランダムサンプリングを行い、予測されたホームチームとアウェイチームの得点を生成
# plot_goals関数を使用してこれらの予測された得点の分布を可視化
ms1 <- rstan::extract(fit1)
plot_pred_comp1 <- function(home_team, away_team, ms) {
old_par <- par(mfrow = c(2, 4))
baseline <- ms$baseline
home_skill <- ms$skill[, which(teams == home_team)]
away_skill <- ms$skill[, which(teams == away_team)]
home_goals_sim <- rpois(length(ms$lp__), exp(baseline + home_skill - away_skill))
away_goals_sim <- rpois(length(ms$lp__), exp(baseline + away_skill - home_skill))
plot_goals(home_goals_sim, away_goals_sim)
home_goals_obs <- d$HomeGoals[ d$HomeTeam == home_team & d$AwayTeam == away_team]
away_goals_obs <- d$AwayGoals[ d$HomeTeam == home_team & d$AwayTeam == away_team]
plot_goals(home_goals_obs, away_goals_obs)
par(old_par)
}
plot_goals <- function(home_goals, away_goals) {
n_matches <- length(home_goals)
goal_diff <- home_goals - away_goals
match_result <- ifelse(goal_diff < 0, "away_win", ifelse(goal_diff > 0, "home_win", "equal"))
hist(home_goals, xlim=c(-0.5, 10), breaks=(0:100) - 0.5)
hist(away_goals, xlim=c(-0.5, 10), breaks=(0:100) - 0.5)
hist(goal_diff, xlim=c(-6, 6), breaks=(-100:100) - 0.5 )
barplot(table(match_result) / n_matches , ylim=c(0, 1))
}
バレンシア(Home) vs
セビージャ(Away)の結果をシミュレートし、laliga
データの任意の試合の実際の結果と一緒に予測結果をプロットする。
plot_pred_comp1(home_team = 'FC Valencia', away_team = 'FC Sevilla', ms = ms1)
過去データはサンプリングしているので、あまりレコードがない。そのため、シュミレーションと似ているかの判断ができないが、シミュレーションではバレンシアの方がセビージャよりわずかに高い確率で勝つことを示している。ホームとアウェイを入れ替えて、セビージャ(Home) vs バレンシア(Away)とする。
plot_pred_comp1(home_team = 'FC Sevilla', away_team = 'FC Valencia', ms = ms1)
ここで現在のモデルに問題があることがわかる。シミュレーションデータは、ホームチームとアウェイチームが入れ替えても結果は変わらない。過去のデータを見ると、セビージャがバレンシアに勝つのはホームチームである場合が多い。ホームチームであることの利点を考慮しているモデルが必要である。
ホームアドバンテージを考慮するためのモデルの変更は、ベースラインをホームベースラインとアウェイベースラインの2つに分けること。
parameters {
real home_baseline; // ホームチームのベースラインのパラメータ
real away_baseline; // アウェイチームのベースラインのパラメータ
real skill[n_teams]; // 各チームのスキル
real<lower=0> group_sigma; // グループの標準偏差
// real group_skill; // グループのスキル
}
transformed parameters {
matrix[n_teams, n_teams] lambda_home; // ホームチームの得点率行列
matrix[n_teams, n_teams] lambda_away; // アウェイチームの得点率行列
for (home_i in 1:n_teams) {
for (away_i in 1:n_teams) {
lambda_home[home_i, away_i] = exp(home_baseline + skill[home_i] - skill[away_i]);
lambda_away[home_i, away_i] = exp(away_baseline + skill[away_i] - skill[home_i]);
}
}
}
変更したモデルでサンプリングする。
model2 <- stan_model('model2.stan')
sampling()
関数でサンプリングする。
fit2 <- sampling(object = model2, data = data_list, seed = 1989, chains = 4, iter = 10000, thin = 2)
推定結果を確認する。
ms2 <- rstan::extract(fit2)
print(fit2, prob = c(0.025, 0.5, 0.975), digits_summary = 2)
## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=5000; thin=2;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## home_baseline 0.56 0.00 0.05 0.47 0.57 0.66 9955 1
## away_baseline 0.22 0.00 0.06 0.10 0.22 0.33 9344 1
## skill[1] -0.20 0.00 0.16 -0.53 -0.20 0.11 4271 1
## skill[2] -0.01 0.00 0.15 -0.32 -0.01 0.30 4216 1
## skill[3] 0.59 0.00 0.15 0.28 0.58 0.90 4293 1
## skill[4] -0.05 0.00 0.15 -0.36 -0.05 0.26 4284 1
## skill[5] -0.02 0.00 0.15 -0.33 -0.02 0.29 4233 1
## skill[6] 0.38 0.00 0.15 0.08 0.38 0.70 4216 1
## skill[7] 0.00 0.00 0.15 -0.32 0.00 0.30 4100 1
## skill[8] -0.26 0.00 0.16 -0.58 -0.26 0.05 4264 1
## skill[9] -0.45 0.00 0.18 -0.82 -0.44 -0.11 4808 1
## group_sigma 0.39 0.00 0.13 0.22 0.36 0.71 6228 1
## lambda_home[1,1] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[1,2] 1.46 0.00 0.18 1.15 1.45 1.83 8555 1
## lambda_home[1,3] 0.81 0.00 0.11 0.62 0.80 1.03 8452 1
## lambda_home[1,4] 1.52 0.00 0.19 1.18 1.51 1.93 8677 1
## lambda_home[1,5] 1.47 0.00 0.18 1.15 1.47 1.85 8239 1
## lambda_home[1,6] 0.99 0.00 0.13 0.77 0.98 1.26 8507 1
## lambda_home[1,7] 1.45 0.00 0.18 1.14 1.44 1.83 8570 1
## lambda_home[1,8] 1.88 0.00 0.23 1.47 1.87 2.35 8212 1
## lambda_home[1,9] 2.28 0.00 0.36 1.66 2.25 3.05 8548 1
## lambda_home[2,1] 2.14 0.00 0.24 1.71 2.13 2.66 8522 1
## lambda_home[2,2] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[2,3] 0.98 0.00 0.11 0.78 0.97 1.20 9032 1
## lambda_home[2,4] 1.84 0.00 0.19 1.49 1.83 2.25 8812 1
## lambda_home[2,5] 1.78 0.00 0.19 1.45 1.77 2.16 8957 1
## lambda_home[2,6] 1.20 0.00 0.13 0.96 1.19 1.47 8955 1
## lambda_home[2,7] 1.76 0.00 0.18 1.43 1.75 2.15 8885 1
## lambda_home[2,8] 2.27 0.00 0.24 1.83 2.26 2.77 8809 1
## lambda_home[2,9] 2.75 0.00 0.40 2.05 2.72 3.63 8641 1
## lambda_home[3,1] 3.88 0.00 0.44 3.10 3.87 4.80 7878 1
## lambda_home[3,2] 3.20 0.00 0.31 2.64 3.19 3.86 9034 1
## lambda_home[3,3] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[3,4] 3.34 0.00 0.34 2.72 3.32 4.03 8825 1
## lambda_home[3,5] 3.23 0.00 0.31 2.65 3.22 3.86 9122 1
## lambda_home[3,6] 2.17 0.00 0.21 1.78 2.16 2.60 8972 1
## lambda_home[3,7] 3.18 0.00 0.30 2.61 3.17 3.80 9082 1
## lambda_home[3,8] 4.11 0.00 0.41 3.35 4.09 4.97 9172 1
## lambda_home[3,9] 4.99 0.01 0.74 3.72 4.94 6.55 9198 1
## lambda_home[4,1] 2.06 0.00 0.25 1.61 2.04 2.58 8397 1
## lambda_home[4,2] 1.70 0.00 0.18 1.37 1.69 2.08 8671 1
## lambda_home[4,3] 0.94 0.00 0.11 0.74 0.93 1.17 8529 1
## lambda_home[4,4] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[4,5] 1.71 0.00 0.19 1.37 1.70 2.10 8379 1
## lambda_home[4,6] 1.15 0.00 0.13 0.91 1.14 1.43 8416 1
## lambda_home[4,7] 1.69 0.00 0.18 1.35 1.68 2.08 8546 1
## lambda_home[4,8] 2.18 0.00 0.24 1.74 2.17 2.69 8171 1
## lambda_home[4,9] 2.65 0.00 0.40 1.95 2.62 3.50 8822 1
## lambda_home[5,1] 2.13 0.00 0.25 1.68 2.11 2.65 7982 1
## lambda_home[5,2] 1.76 0.00 0.18 1.42 1.75 2.14 8979 1
## lambda_home[5,3] 0.97 0.00 0.11 0.77 0.96 1.20 8980 1
## lambda_home[5,4] 1.83 0.00 0.20 1.47 1.82 2.24 8700 1
## lambda_home[5,5] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[5,6] 1.19 0.00 0.13 0.95 1.18 1.46 8870 1
## lambda_home[5,7] 1.74 0.00 0.18 1.42 1.74 2.13 8991 1
## lambda_home[5,8] 2.25 0.00 0.24 1.81 2.24 2.74 8838 1
## lambda_home[5,9] 2.73 0.00 0.41 2.03 2.70 3.61 8675 1
## lambda_home[6,1] 3.17 0.00 0.36 2.52 3.14 3.92 7710 1
## lambda_home[6,2] 2.61 0.00 0.26 2.13 2.60 3.16 8658 1
## lambda_home[6,3] 1.44 0.00 0.15 1.17 1.43 1.75 8715 1
## lambda_home[6,4] 2.72 0.00 0.28 2.21 2.71 3.31 8498 1
## lambda_home[6,5] 2.63 0.00 0.26 2.16 2.62 3.17 8749 1
## lambda_home[6,6] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[6,7] 2.60 0.00 0.25 2.13 2.58 3.14 8678 1
## lambda_home[6,8] 3.35 0.00 0.34 2.72 3.34 4.05 8888 1
## lambda_home[6,9] 4.07 0.01 0.60 3.02 4.04 5.37 8736 1
## lambda_home[7,1] 2.16 0.00 0.25 1.72 2.14 2.68 8309 1
## lambda_home[7,2] 1.78 0.00 0.18 1.44 1.77 2.16 8896 1
## lambda_home[7,3] 0.98 0.00 0.11 0.79 0.98 1.21 8703 1
## lambda_home[7,4] 1.85 0.00 0.20 1.49 1.84 2.27 8634 1
## lambda_home[7,5] 1.79 0.00 0.19 1.45 1.79 2.18 8944 1
## lambda_home[7,6] 1.20 0.00 0.13 0.97 1.20 1.47 8798 1
## lambda_home[7,7] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[7,8] 2.28 0.00 0.24 1.85 2.27 2.79 8867 1
## lambda_home[7,9] 2.77 0.00 0.41 2.05 2.74 3.65 8977 1
## lambda_home[8,1] 1.67 0.00 0.21 1.31 1.66 2.10 7736 1
## lambda_home[8,2] 1.38 0.00 0.16 1.09 1.37 1.71 8870 1
## lambda_home[8,3] 0.76 0.00 0.09 0.60 0.76 0.96 9045 1
## lambda_home[8,4] 1.44 0.00 0.17 1.14 1.43 1.78 8188 1
## lambda_home[8,5] 1.39 0.00 0.16 1.10 1.38 1.72 8721 1
## lambda_home[8,6] 0.93 0.00 0.11 0.74 0.93 1.17 9011 1
## lambda_home[8,7] 1.37 0.00 0.15 1.09 1.36 1.70 8974 1
## lambda_home[8,8] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_home[8,9] 2.15 0.00 0.33 1.58 2.12 2.85 8863 1
## lambda_home[9,1] 1.39 0.00 0.22 1.00 1.37 1.87 8619 1
## lambda_home[9,2] 1.15 0.00 0.18 0.84 1.13 1.53 8947 1
## lambda_home[9,3] 0.63 0.00 0.11 0.45 0.63 0.86 9401 1
## lambda_home[9,4] 1.20 0.00 0.19 0.87 1.18 1.60 9279 1
## lambda_home[9,5] 1.16 0.00 0.18 0.84 1.15 1.55 9165 1
## lambda_home[9,6] 0.78 0.00 0.13 0.56 0.77 1.05 9167 1
## lambda_home[9,7] 1.14 0.00 0.18 0.83 1.13 1.53 9275 1
## lambda_home[9,8] 1.47 0.00 0.23 1.08 1.46 1.97 9067 1
## lambda_home[9,9] 1.76 0.00 0.09 1.59 1.76 1.94 9953 1
## lambda_away[1,1] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[1,2] 1.51 0.00 0.18 1.20 1.51 1.88 8589 1
## lambda_away[1,3] 2.74 0.00 0.32 2.18 2.73 3.41 8396 1
## lambda_away[1,4] 1.46 0.00 0.18 1.13 1.44 1.83 8612 1
## lambda_away[1,5] 1.50 0.00 0.18 1.18 1.49 1.88 8102 1
## lambda_away[1,6] 2.24 0.00 0.26 1.77 2.23 2.77 8166 1
## lambda_away[1,7] 1.52 0.00 0.18 1.20 1.52 1.90 8607 1
## lambda_away[1,8] 1.18 0.00 0.15 0.92 1.17 1.49 8421 1
## lambda_away[1,9] 0.98 0.00 0.16 0.70 0.97 1.33 8717 1
## lambda_away[2,1] 1.03 0.00 0.13 0.80 1.03 1.30 8150 1
## lambda_away[2,2] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[2,3] 2.26 0.00 0.23 1.85 2.25 2.75 8919 1
## lambda_away[2,4] 1.20 0.00 0.13 0.96 1.19 1.48 8679 1
## lambda_away[2,5] 1.24 0.00 0.14 0.99 1.24 1.52 8866 1
## lambda_away[2,6] 1.85 0.00 0.19 1.50 1.84 2.24 8444 1
## lambda_away[2,7] 1.26 0.00 0.14 1.01 1.25 1.54 8778 1
## lambda_away[2,8] 0.98 0.00 0.11 0.77 0.97 1.22 8835 1
## lambda_away[2,9] 0.81 0.00 0.13 0.59 0.80 1.09 8867 1
## lambda_away[3,1] 0.57 0.00 0.08 0.44 0.57 0.73 7770 1
## lambda_away[3,2] 0.69 0.00 0.08 0.55 0.69 0.86 8914 1
## lambda_away[3,3] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[3,4] 0.66 0.00 0.08 0.52 0.66 0.83 8165 1
## lambda_away[3,5] 0.68 0.00 0.08 0.54 0.68 0.85 8568 1
## lambda_away[3,6] 1.02 0.00 0.11 0.82 1.02 1.24 8501 1
## lambda_away[3,7] 0.69 0.00 0.08 0.55 0.69 0.86 8505 1
## lambda_away[3,8] 0.54 0.00 0.07 0.42 0.53 0.68 9154 1
## lambda_away[3,9] 0.45 0.00 0.08 0.32 0.44 0.61 9313 1
## lambda_away[4,1] 1.08 0.00 0.14 0.83 1.07 1.38 8466 1
## lambda_away[4,2] 1.30 0.00 0.14 1.04 1.29 1.60 8662 1
## lambda_away[4,3] 2.36 0.00 0.25 1.91 2.35 2.88 8541 1
## lambda_away[4,4] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[4,5] 1.29 0.00 0.14 1.03 1.28 1.59 8475 1
## lambda_away[4,6] 1.92 0.00 0.21 1.55 1.91 2.36 8216 1
## lambda_away[4,7] 1.31 0.00 0.15 1.04 1.30 1.62 8501 1
## lambda_away[4,8] 1.02 0.00 0.12 0.80 1.01 1.27 8310 1
## lambda_away[4,9] 0.85 0.00 0.14 0.61 0.83 1.14 9118 1
## lambda_away[5,1] 1.04 0.00 0.13 0.81 1.04 1.32 7843 1
## lambda_away[5,2] 1.26 0.00 0.14 1.01 1.25 1.54 8916 1
## lambda_away[5,3] 2.28 0.00 0.23 1.86 2.28 2.75 8963 1
## lambda_away[5,4] 1.21 0.00 0.14 0.97 1.20 1.50 8450 1
## lambda_away[5,5] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[5,6] 1.86 0.00 0.19 1.51 1.85 2.26 8595 1
## lambda_away[5,7] 1.27 0.00 0.14 1.01 1.26 1.55 8855 1
## lambda_away[5,8] 0.98 0.00 0.11 0.78 0.98 1.22 8957 1
## lambda_away[5,9] 0.82 0.00 0.13 0.59 0.81 1.10 9099 1
## lambda_away[6,1] 0.70 0.00 0.09 0.54 0.69 0.90 8193 1
## lambda_away[6,2] 0.85 0.00 0.10 0.67 0.84 1.05 8896 1
## lambda_away[6,3] 1.53 0.00 0.16 1.24 1.52 1.86 8899 1
## lambda_away[6,4] 0.81 0.00 0.10 0.64 0.81 1.02 8587 1
## lambda_away[6,5] 0.84 0.00 0.10 0.66 0.84 1.04 8895 1
## lambda_away[6,6] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[6,7] 0.85 0.00 0.10 0.68 0.85 1.05 8807 1
## lambda_away[6,8] 0.66 0.00 0.08 0.51 0.66 0.83 9038 1
## lambda_away[6,9] 0.55 0.00 0.09 0.39 0.54 0.75 9164 1
## lambda_away[7,1] 1.03 0.00 0.13 0.80 1.02 1.30 8201 1
## lambda_away[7,2] 1.24 0.00 0.13 1.00 1.24 1.52 8875 1
## lambda_away[7,3] 2.25 0.00 0.22 1.84 2.24 2.70 8892 1
## lambda_away[7,4] 1.19 0.00 0.13 0.95 1.19 1.47 8556 1
## lambda_away[7,5] 1.23 0.00 0.13 0.99 1.23 1.51 8935 1
## lambda_away[7,6] 1.83 0.00 0.18 1.50 1.83 2.23 8427 1
## lambda_away[7,7] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[7,8] 0.97 0.00 0.11 0.77 0.96 1.21 8944 1
## lambda_away[7,9] 0.81 0.00 0.13 0.58 0.80 1.09 9203 1
## lambda_away[8,1] 1.33 0.00 0.17 1.03 1.32 1.68 7584 1
## lambda_away[8,2] 1.60 0.00 0.18 1.27 1.60 1.97 8755 1
## lambda_away[8,3] 2.90 0.00 0.30 2.35 2.89 3.55 9024 1
## lambda_away[8,4] 1.54 0.00 0.18 1.22 1.53 1.91 8090 1
## lambda_away[8,5] 1.59 0.00 0.18 1.27 1.58 1.96 8475 1
## lambda_away[8,6] 2.37 0.00 0.25 1.91 2.36 2.88 8688 1
## lambda_away[8,7] 1.61 0.00 0.18 1.29 1.61 1.99 8481 1
## lambda_away[8,8] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lambda_away[8,9] 1.04 0.00 0.17 0.75 1.03 1.40 8863 1
## lambda_away[9,1] 1.61 0.00 0.26 1.17 1.59 2.17 8375 1
## lambda_away[9,2] 1.95 0.00 0.29 1.44 1.93 2.57 8564 1
## lambda_away[9,3] 3.53 0.01 0.53 2.61 3.49 4.66 9121 1
## lambda_away[9,4] 1.87 0.00 0.29 1.37 1.85 2.48 8834 1
## lambda_away[9,5] 1.93 0.00 0.29 1.41 1.91 2.56 8637 1
## lambda_away[9,6] 2.88 0.00 0.43 2.12 2.85 3.78 8663 1
## lambda_away[9,7] 1.96 0.00 0.29 1.45 1.93 2.58 8915 1
## lambda_away[9,8] 1.52 0.00 0.23 1.11 1.50 2.02 8881 1
## lambda_away[9,9] 1.24 0.00 0.07 1.11 1.24 1.39 9343 1
## lp__ -304.47 0.03 2.49 -310.21 -304.12 -300.61 6134 1
##
## Samples were drawn using NUTS(diag_e) at Tue Jan 23 16:05:56 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).
home_baseline
とaway_baseline
のトレースプロットと分布を見ると、かなりのホームアドバンテージがあることがわかる。ホームベースラインとは0.564473で、アウェイベースラインは0.2169645である。
plt_e <- stan_trace(fit2,
pars = 'home_baseline',
separate_chains = TRUE) + ggtitle('Traceplot of home_baseline')
plt_f <- stan_dens(fit2,
pars = 'home_baseline',
separate_chains = TRUE)
plt_g <- stan_trace(fit2,
pars = 'away_baseline',
separate_chains = TRUE) + ggtitle('Traceplot of away_baseline')
plt_h <- stan_dens(fit2,
pars = 'away_baseline',
separate_chains = TRUE)
(plt_e | plt_f) / (plt_g | plt_h)
exp(home_baseline)
とexp(away_baseline)
の差を見ると、ホームアドバンテージはホームチームのゴール数でおよそ0.5162821点多いことがわかる。
plotPost(
paramSampleVec = exp(ms2$home_baseline) - exp(ms2$away_baseline),
compVal = 0,
xlab = 'Home advantage in number of goals'
)
最後に、バレンシア(ホームチーム)対セビージャ(アウェイチーム)のシミュレーション結果を見てみよう。
plot_pred_comp2 <- function(home_team, away_team, ms) {
old_par <- par(mfrow = c(2, 4))
home_baseline <- ms$home_baseline
away_baseline <- ms$away_baseline
home_skill <- ms$skill[, which(teams == home_team)]
away_skill <- ms$skill[, which(teams == away_team)]
home_goals_sim <- rpois(length(ms$lp__), exp(home_baseline + home_skill - away_skill))
away_goals_sim <- rpois(length(ms$lp__), exp(away_baseline + away_skill - home_skill))
plot_goals(home_goals_sim, away_goals_sim)
home_goals_obs <- d$HomeGoals[ d$HomeTeam == home_team & d$AwayTeam == away_team]
away_goals_obs <- d$AwayGoals[ d$HomeTeam == home_team & d$AwayTeam == away_team]
plot_goals(home_goals_obs, away_goals_obs)
par(old_par)
}
バレンシア(Home) vs セビージャ(Away)を可視化する。
plot_pred_comp2(home_team = 'FC Valencia', away_team = 'FC Sevilla', ms = ms2)
そして同様に、チームを入れ替えて、セビージャ(Home) vs バレンシア(Away)として可視化する。
plot_pred_comp2(home_team = 'FC Sevilla', away_team = 'FC Valencia', ms = ms2)
セビージャもバレンシアもホームチームとしてプレーした方が勝ちやすいことがわかる。
データセットlaliga
には5つの異なるシーズンのデータが含まれており、現在のモデルの仮定は、チームがすべてのシーズンで同じスキルを持つとしている。これはおそらく現実的な仮定ではなく、チームはおそらく年ごとにパフォーマンスが異なるはずである。1部リーグから脱落した場合、laliga
データのすべてのシーズンに参加していないチームもある。
Villarreal CF
:
2012/13
シーズンは1部リーグではないReal Zaragoza
:
2008/09
シーズンは1部リーグではないReal Valladolid
:
2010/11,2011/12
シーズンは1部リーグではないCD Tenerife
:
2009/10
シーズンのみ1部リーグplt_particicipation
チームのスキルの年ごとの変動を含めるように修正する。これは、各チームがシーズンごとに1つのスキル・パラメーターを持つようにし、シーズン\(t\)のチームのスキル・パラメタをシーズン\(t+1\)のチームのスキル・パラメタの事前分布に使用することによって、スキル・パラメタを引き継ぐことで行う。
\[ \text{skill}_{t+1} \sim \text{Normal}(\text{skill}_{t}, \sigma_{\text{season}}^2) \]
ここで\(\sigma^{2}_{season}\)はデータ全体を使って推定したパラメタである。
model3 <- stan_model('model3.stan')
モデルはこちら。
data {
int<lower=1> n_games;
int<lower=1> n_seasons;
int<lower=1> n_teams;
int<lower=1> Season[n_games];
int<lower=1> HomeTeam[n_games];
int<lower=1> AwayTeam[n_games];
int<lower=0> HomeGoals[n_games];
int<lower=0> AwayGoals[n_games];
}
parameters {
vector[n_seasons] home_baseline;
vector[n_seasons] away_baseline;
matrix[n_teams, n_seasons] skill;
real<lower=0> group_sigma;
// real<lower=0> group_skill;
real<lower=0> season_sigma;
}
transformed parameters {
vector[n_games] lambda_home;
vector[n_games] lambda_away;
for (i in 1:n_games) {
lambda_home[i] = exp(home_baseline[Season[i]] + skill[HomeTeam[i], Season[i]] - skill[AwayTeam[i], Season[i]]);
lambda_away[i] = exp(away_baseline[Season[i]] + skill[AwayTeam[i], Season[i]] - skill[HomeTeam[i], Season[i]]);
}
}
model {
group_sigma ~ uniform(0, 1000);
// group_skill ~ normal(0, 4);
home_baseline[1] ~ normal(0, 4);
away_baseline[1] ~ normal(0, 4);
for (n in 1:n_teams) {
skill[n, 1] ~ normal(0, group_sigma);
}
for (season_i in 2:n_seasons) {
for (n in 1:n_teams) {
skill[n, season_i] ~ normal(skill[n, season_i-1], season_sigma);
}
home_baseline[season_i] ~ normal(home_baseline[season_i-1], season_sigma);
away_baseline[season_i] ~ normal(away_baseline[season_i-1], season_sigma);
}
for (i in 1:n_games) {
HomeGoals[i] ~ poisson(lambda_home[i]);
AwayGoals[i] ~ poisson(lambda_away[i]);
}
}
sampling()
関数でサンプリングする。
fit3 <- sampling(object = model3, data = data_list, seed = 1989, chains = 4, iter = 10000, thin = 5)
推定結果を確認する。
print(fit3, prob = c(0.025, 0.5, 0.975), digits_summary = 3)
## Inference for Stan model: anon_model.
## 4 chains, each with iter=10000; warmup=5000; thin=5;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## home_baseline[1] 0.617 0.001 0.080 0.468 0.614 0.781 3338 1.002
## home_baseline[2] 0.615 0.001 0.066 0.490 0.614 0.748 2207 1.002
## home_baseline[3] 0.540 0.001 0.069 0.395 0.543 0.665 3079 1.000
## home_baseline[4] 0.486 0.002 0.083 0.306 0.492 0.634 1975 1.001
## home_baseline[5] 0.479 0.002 0.093 0.285 0.485 0.641 2391 1.000
## away_baseline[1] 0.256 0.001 0.087 0.090 0.254 0.431 3572 1.000
## away_baseline[2] 0.241 0.001 0.072 0.101 0.241 0.380 3472 1.000
## away_baseline[3] 0.197 0.001 0.077 0.034 0.198 0.340 3392 1.000
## away_baseline[4] 0.172 0.001 0.084 -0.007 0.177 0.326 3268 1.000
## away_baseline[5] 0.142 0.002 0.101 -0.078 0.148 0.320 2515 1.000
## skill[1,1] -0.160 0.004 0.162 -0.478 -0.162 0.162 2093 1.001
## skill[1,2] -0.224 0.004 0.163 -0.539 -0.225 0.088 2047 1.001
## skill[1,3] -0.212 0.004 0.176 -0.564 -0.210 0.132 2179 1.001
## skill[1,4] -0.197 0.004 0.187 -0.573 -0.197 0.168 2433 1.001
## skill[1,5] -0.186 0.004 0.187 -0.556 -0.186 0.172 2429 1.000
## skill[2,1] -0.023 0.003 0.158 -0.331 -0.024 0.298 2192 1.000
## skill[2,2] 0.000 0.003 0.156 -0.307 -0.003 0.314 2067 1.000
## skill[2,3] -0.016 0.004 0.163 -0.343 -0.018 0.302 2055 1.000
## skill[2,4] -0.009 0.004 0.168 -0.345 -0.009 0.315 2102 1.001
## skill[2,5] 0.025 0.004 0.180 -0.337 0.022 0.372 2325 1.000
## skill[3,1] 0.586 0.004 0.163 0.283 0.580 0.908 2077 1.001
## skill[3,2] 0.593 0.004 0.157 0.294 0.589 0.903 1895 1.001
## skill[3,3] 0.592 0.004 0.161 0.288 0.587 0.910 2071 1.000
## skill[3,4] 0.591 0.004 0.167 0.278 0.589 0.916 2017 1.001
## skill[3,5] 0.560 0.004 0.176 0.221 0.557 0.899 2215 1.001
## skill[4,1] -0.023 0.003 0.160 -0.329 -0.025 0.295 2216 1.001
## skill[4,2] -0.017 0.004 0.159 -0.322 -0.021 0.291 2001 1.001
## skill[4,3] -0.060 0.004 0.166 -0.391 -0.063 0.269 2048 1.001
## skill[4,4] -0.095 0.004 0.178 -0.454 -0.095 0.248 2035 1.002
## skill[4,5] -0.094 0.004 0.200 -0.496 -0.092 0.290 2217 1.001
## skill[5,1] 0.002 0.003 0.160 -0.313 -0.001 0.311 2138 1.000
## skill[5,2] -0.006 0.003 0.157 -0.315 -0.007 0.294 2016 1.000
## skill[5,3] -0.005 0.004 0.161 -0.315 -0.008 0.312 2052 1.000
## skill[5,4] -0.020 0.004 0.166 -0.347 -0.020 0.301 2057 1.001
## skill[5,5] -0.044 0.004 0.176 -0.399 -0.043 0.301 2107 1.001
## skill[6,1] 0.255 0.005 0.174 -0.087 0.255 0.616 1367 1.002
## skill[6,2] 0.356 0.004 0.156 0.059 0.354 0.673 1899 1.001
## skill[6,3] 0.411 0.004 0.162 0.097 0.406 0.731 2031 1.000
## skill[6,4] 0.474 0.004 0.171 0.155 0.469 0.825 1894 1.000
## skill[6,5] 0.482 0.004 0.180 0.144 0.475 0.854 2001 1.000
## skill[7,1] 0.009 0.003 0.158 -0.298 0.004 0.328 2075 1.000
## skill[7,2] -0.011 0.004 0.157 -0.320 -0.011 0.299 1973 1.000
## skill[7,3] -0.019 0.004 0.163 -0.347 -0.017 0.300 2031 1.001
## skill[7,4] -0.003 0.004 0.166 -0.332 -0.007 0.324 2159 1.001
## skill[7,5] 0.022 0.004 0.179 -0.315 0.017 0.389 2231 1.000
## skill[8,1] -0.217 0.004 0.174 -0.573 -0.216 0.135 2293 1.001
## skill[8,2] -0.235 0.003 0.159 -0.558 -0.235 0.067 2103 1.001
## skill[8,3] -0.236 0.003 0.163 -0.563 -0.236 0.090 2247 1.001
## skill[8,4] -0.282 0.004 0.169 -0.616 -0.280 0.044 2174 1.001
## skill[8,5] -0.304 0.004 0.180 -0.659 -0.305 0.045 2187 1.002
## skill[9,1] -0.394 0.004 0.195 -0.777 -0.391 -0.026 2238 1.002
## skill[9,2] -0.425 0.004 0.178 -0.777 -0.422 -0.084 2245 1.002
## skill[9,3] -0.424 0.004 0.199 -0.810 -0.421 -0.035 2533 1.002
## skill[9,4] -0.422 0.004 0.219 -0.867 -0.423 -0.002 2639 1.002
## skill[9,5] -0.420 0.004 0.237 -0.890 -0.422 0.048 2834 1.002
## group_sigma 0.362 0.003 0.133 0.185 0.336 0.682 1863 1.005
## season_sigma 0.083 0.001 0.036 0.026 0.080 0.161 775 1.005
## lambda_home[1] 1.636 0.005 0.268 1.181 1.614 2.231 3084 1.002
## lambda_home[2] 1.829 0.005 0.280 1.343 1.811 2.441 3242 1.000
## lambda_home[3] 3.441 0.009 0.506 2.559 3.397 4.558 3129 1.000
## lambda_home[4] 1.584 0.004 0.254 1.156 1.561 2.159 3585 1.001
## lambda_home[5] 1.425 0.007 0.252 1.016 1.393 2.010 1442 1.003
## lambda_home[6] 1.876 0.005 0.299 1.356 1.855 2.531 3396 1.000
## lambda_home[7] 1.596 0.005 0.256 1.166 1.569 2.157 3066 1.002
## lambda_home[8] 3.950 0.009 0.581 2.904 3.909 5.162 4098 1.001
## lambda_home[9] 1.246 0.006 0.244 0.867 1.216 1.813 1716 1.004
## lambda_home[10] 1.861 0.004 0.281 1.379 1.840 2.487 3915 1.000
## lambda_home[11] 2.152 0.005 0.333 1.568 2.120 2.871 3999 1.000
## lambda_home[12] 1.048 0.003 0.175 0.754 1.034 1.441 3133 1.004
## lambda_home[13] 3.332 0.007 0.468 2.501 3.304 4.350 4098 1.000
## lambda_home[14] 2.416 0.007 0.372 1.722 2.401 3.165 2676 1.000
## lambda_home[15] 2.613 0.009 0.427 1.910 2.565 3.595 2062 1.001
## lambda_home[16] 1.920 0.005 0.292 1.413 1.901 2.559 3303 1.002
## lambda_home[17] 2.398 0.007 0.363 1.727 2.385 3.140 2831 1.000
## lambda_home[18] 1.023 0.003 0.178 0.723 1.004 1.429 3123 1.003
## lambda_home[19] 1.934 0.005 0.298 1.417 1.910 2.596 3627 1.000
## lambda_home[20] 2.476 0.009 0.383 1.762 2.465 3.260 2028 1.000
## lambda_home[21] 1.933 0.005 0.289 1.431 1.913 2.562 3668 1.001
## lambda_home[22] 2.151 0.005 0.324 1.587 2.125 2.856 3676 1.001
## lambda_home[23] 1.922 0.005 0.302 1.408 1.901 2.597 3544 1.001
## lambda_home[24] 1.022 0.003 0.169 0.723 1.006 1.399 2550 1.003
## lambda_home[25] 2.219 0.005 0.333 1.624 2.204 2.931 3938 1.001
## lambda_home[26] 2.475 0.009 0.379 1.782 2.460 3.268 1796 1.000
## lambda_home[27] 1.874 0.005 0.281 1.377 1.855 2.463 3220 1.001
## lambda_home[28] 2.204 0.006 0.332 1.619 2.177 2.920 3606 1.002
## lambda_home[29] 0.892 0.003 0.160 0.628 0.874 1.265 2757 1.004
## lambda_home[30] 2.844 0.010 0.456 1.997 2.829 3.759 2183 1.001
## lambda_home[31] 1.887 0.005 0.288 1.388 1.860 2.509 3995 1.000
## lambda_home[32] 1.635 0.005 0.261 1.194 1.611 2.228 3207 1.002
## lambda_home[33] 3.357 0.007 0.487 2.517 3.314 4.412 4291 1.000
## lambda_home[34] 1.461 0.006 0.265 1.038 1.428 2.076 1743 1.004
## lambda_home[35] 1.054 0.003 0.174 0.763 1.038 1.446 3048 1.004
## lambda_home[36] 1.348 0.004 0.210 0.958 1.340 1.774 3100 0.999
## lambda_home[37] 1.830 0.005 0.288 1.330 1.806 2.470 3595 1.000
## lambda_home[38] 3.441 0.008 0.501 2.515 3.411 4.547 3929 1.000
## lambda_home[39] 1.471 0.006 0.266 1.053 1.438 2.063 1778 1.004
## lambda_home[40] 1.815 0.004 0.270 1.344 1.796 2.401 3707 1.000
## lambda_home[41] 1.427 0.006 0.267 1.002 1.392 2.037 1702 1.003
## lambda_home[42] 1.817 0.005 0.282 1.342 1.789 2.449 3860 1.000
## lambda_home[43] 1.856 0.004 0.246 1.430 1.834 2.381 3442 1.000
## lambda_home[44] 2.265 0.006 0.366 1.625 2.235 3.055 3368 1.000
## lambda_home[45] 2.346 0.005 0.310 1.792 2.332 3.027 3187 1.000
## lambda_home[46] 1.509 0.003 0.212 1.139 1.498 1.943 3984 1.000
## lambda_home[47] 3.373 0.009 0.424 2.615 3.341 4.317 2395 1.000
## lambda_home[48] 1.852 0.006 0.294 1.369 1.823 2.521 2643 1.001
## lambda_home[49] 1.287 0.004 0.190 0.958 1.269 1.699 2254 1.001
## lambda_home[50] 1.845 0.004 0.243 1.410 1.828 2.353 3165 1.000
## lambda_home[51] 4.085 0.017 0.642 2.974 4.035 5.440 1407 1.001
## lambda_home[52] 2.359 0.006 0.314 1.809 2.336 3.041 3120 1.001
## lambda_home[53] 1.301 0.004 0.183 0.996 1.284 1.707 2153 1.002
## lambda_home[54] 1.021 0.003 0.142 0.773 1.009 1.329 2810 1.002
## lambda_home[55] 3.333 0.008 0.453 2.550 3.297 4.305 3513 1.000
## lambda_home[56] 4.266 0.009 0.547 3.301 4.224 5.437 3778 1.000
## lambda_home[57] 2.815 0.009 0.454 2.041 2.771 3.840 2437 1.000
## lambda_home[58] 1.308 0.005 0.185 0.983 1.296 1.711 1536 1.002
## lambda_home[59] 1.887 0.004 0.246 1.447 1.872 2.408 3850 1.001
## lambda_home[60] 2.334 0.005 0.313 1.774 2.317 3.003 4081 1.000
## lambda_home[61] 1.234 0.004 0.210 0.866 1.219 1.682 2434 1.001
## lambda_home[62] 2.300 0.007 0.352 1.695 2.267 3.115 2408 1.000
## lambda_home[63] 2.364 0.006 0.308 1.842 2.335 3.047 2988 1.001
## lambda_home[64] 2.290 0.007 0.375 1.637 2.263 3.138 2832 1.000
## lambda_home[65] 2.324 0.006 0.336 1.761 2.290 3.079 2960 1.001
## lambda_home[66] 1.898 0.005 0.253 1.437 1.884 2.440 3126 1.001
## lambda_home[67] 1.294 0.004 0.181 0.974 1.281 1.689 2443 1.002
## lambda_home[68] 1.226 0.004 0.204 0.874 1.213 1.659 3000 1.001
## lambda_home[69] 3.366 0.007 0.432 2.608 3.331 4.299 3898 1.000
## lambda_home[70] 1.877 0.005 0.258 1.426 1.854 2.455 2294 1.001
## lambda_home[71] 3.431 0.007 0.428 2.680 3.406 4.352 3901 1.000
## lambda_home[72] 0.679 0.003 0.120 0.469 0.671 0.928 2258 1.002
## lambda_home[73] 1.493 0.004 0.210 1.119 1.478 1.938 3555 1.000
## lambda_home[74] 3.392 0.007 0.429 2.642 3.358 4.313 3808 1.000
## lambda_home[75] 1.878 0.004 0.247 1.441 1.864 2.409 3969 1.001
## lambda_home[76] 1.241 0.004 0.213 0.863 1.229 1.695 2351 1.000
## lambda_home[77] 0.826 0.002 0.122 0.611 0.817 1.088 3506 1.001
## lambda_home[78] 2.320 0.006 0.318 1.754 2.302 3.001 3156 1.000
## lambda_home[79] 1.550 0.004 0.254 1.086 1.530 2.091 3218 1.000
## lambda_home[80] 1.875 0.005 0.250 1.431 1.852 2.429 2910 1.001
## lambda_home[81] 2.312 0.006 0.328 1.728 2.281 3.033 3287 1.001
## lambda_home[82] 1.486 0.004 0.217 1.120 1.465 1.975 2563 1.001
## lambda_home[83] 1.032 0.003 0.146 0.775 1.023 1.354 1878 1.004
## lambda_home[84] 1.889 0.004 0.270 1.393 1.874 2.451 4046 0.999
## lambda_home[85] 2.707 0.006 0.340 2.092 2.684 3.441 3732 1.000
## lambda_home[86] 1.887 0.005 0.256 1.432 1.868 2.452 2176 1.001
## lambda_home[87] 0.861 0.004 0.152 0.594 0.851 1.186 1813 1.001
## lambda_home[88] 2.676 0.006 0.341 2.083 2.655 3.424 3675 1.000
## lambda_home[89] 1.477 0.004 0.212 1.123 1.463 1.942 3144 1.000
## lambda_home[90] 3.410 0.008 0.432 2.659 3.392 4.326 3094 1.000
## lambda_home[91] 1.047 0.002 0.153 0.782 1.035 1.374 3798 1.001
## lambda_home[92] 0.818 0.003 0.125 0.599 0.806 1.092 2288 1.003
## lambda_home[93] 1.248 0.004 0.207 0.888 1.232 1.696 2950 1.001
## lambda_home[94] 2.661 0.007 0.334 2.066 2.634 3.385 2475 0.999
## lambda_home[95] 1.494 0.004 0.221 1.121 1.476 1.984 2499 1.001
## lambda_home[96] 1.845 0.005 0.260 1.404 1.823 2.440 2644 1.001
## lambda_home[97] 2.847 0.008 0.449 2.065 2.805 3.817 3178 1.000
## lambda_home[98] 1.519 0.003 0.215 1.138 1.503 1.985 3930 1.000
## lambda_home[99] 1.537 0.006 0.271 1.061 1.518 2.119 2241 1.001
## lambda_home[100] 1.471 0.003 0.191 1.129 1.457 1.897 3476 1.002
## lambda_home[101] 1.501 0.003 0.211 1.127 1.490 1.958 3992 1.000
## lambda_home[102] 2.690 0.006 0.342 2.081 2.671 3.424 3442 1.000
## lambda_home[103] 1.835 0.005 0.256 1.381 1.819 2.368 2774 1.000
## lambda_home[104] 2.862 0.007 0.448 2.082 2.828 3.813 3594 1.001
## lambda_home[105] 1.036 0.004 0.160 0.761 1.022 1.402 1889 1.002
## lambda_home[106] 1.855 0.005 0.248 1.428 1.838 2.379 3011 1.000
## lambda_home[107] 1.015 0.003 0.148 0.757 1.002 1.342 2405 1.002
## lambda_home[108] 2.338 0.007 0.342 1.746 2.304 3.117 2089 1.002
## lambda_home[109] 5.179 0.020 0.816 3.774 5.108 6.954 1611 1.001
## lambda_home[110] 1.026 0.003 0.143 0.776 1.017 1.337 2954 1.003
## lambda_home[111] 1.855 0.005 0.264 1.409 1.835 2.435 2511 1.000
## lambda_home[112] 4.226 0.011 0.587 3.210 4.174 5.540 3034 1.000
## lambda_home[113] 2.833 0.008 0.453 2.045 2.792 3.846 2882 1.000
## lambda_home[114] 1.503 0.004 0.213 1.150 1.481 1.977 3533 1.001
## lambda_home[115] 0.944 0.003 0.139 0.681 0.936 1.239 2680 1.000
## lambda_home[116] 1.728 0.004 0.238 1.284 1.724 2.226 3659 1.000
## lambda_home[117] 2.157 0.005 0.296 1.608 2.145 2.763 3237 1.000
## lambda_home[118] 1.752 0.004 0.238 1.322 1.738 2.260 3908 1.000
## lambda_home[119] 3.185 0.006 0.402 2.454 3.171 4.006 4168 1.000
## lambda_home[120] 0.758 0.002 0.113 0.558 0.752 0.997 3093 1.000
## lambda_home[121] 1.658 0.004 0.233 1.236 1.648 2.150 3669 1.000
## lambda_home[122] 3.141 0.006 0.403 2.389 3.123 3.998 4017 0.999
## lambda_home[123] 2.153 0.005 0.297 1.587 2.148 2.758 3120 1.000
## lambda_home[124] 1.756 0.004 0.238 1.328 1.740 2.273 4072 1.000
## lambda_home[125] 2.652 0.006 0.348 2.026 2.626 3.403 3699 1.000
## lambda_home[126] 3.322 0.007 0.431 2.550 3.298 4.254 3871 1.000
## lambda_home[127] 1.375 0.003 0.196 1.031 1.364 1.790 3425 1.000
## lambda_home[128] 1.662 0.004 0.232 1.240 1.648 2.155 3907 0.999
## lambda_home[129] 2.074 0.005 0.273 1.578 2.061 2.652 3490 1.000
## lambda_home[130] 1.454 0.003 0.209 1.092 1.443 1.925 4019 1.000
## lambda_home[131] 2.658 0.006 0.348 2.044 2.631 3.434 3865 1.001
## lambda_home[132] 1.639 0.004 0.232 1.220 1.626 2.125 3648 0.999
## lambda_home[133] 0.908 0.002 0.134 0.669 0.901 1.196 3327 1.001
## lambda_home[134] 2.620 0.005 0.341 2.015 2.598 3.351 4042 1.000
## lambda_home[135] 2.772 0.006 0.370 2.135 2.741 3.572 3710 1.001
## lambda_home[136] 3.179 0.008 0.414 2.418 3.157 4.073 2985 1.000
## lambda_home[137] 1.736 0.004 0.241 1.304 1.723 2.235 3445 1.000
## lambda_home[138] 1.391 0.003 0.198 1.041 1.382 1.827 4004 1.000
## lambda_home[139] 1.712 0.004 0.239 1.272 1.701 2.202 3443 1.000
## lambda_home[140] 0.941 0.002 0.135 0.697 0.936 1.228 3134 1.000
## lambda_home[141] 1.810 0.004 0.252 1.345 1.800 2.339 3330 1.001
## lambda_home[142] 3.958 0.009 0.505 3.011 3.943 4.989 3374 1.000
## lambda_home[143] 0.954 0.002 0.136 0.715 0.945 1.248 3460 1.000
## lambda_home[144] 1.395 0.003 0.201 1.045 1.379 1.826 3486 1.000
## lambda_home[145] 1.131 0.003 0.166 0.828 1.125 1.468 2793 1.000
## lambda_home[146] 1.708 0.004 0.236 1.269 1.702 2.206 3875 1.000
## lambda_home[147] 2.182 0.005 0.292 1.646 2.165 2.792 3562 1.000
## lambda_home[148] 0.903 0.002 0.133 0.653 0.896 1.188 3053 1.000
## lambda_home[149] 1.806 0.004 0.248 1.348 1.792 2.336 3776 1.000
## lambda_home[150] 1.444 0.003 0.192 1.101 1.429 1.856 3754 1.000
## lambda_home[151] 2.065 0.005 0.289 1.520 2.056 2.659 3463 1.000
## lambda_home[152] 1.831 0.004 0.248 1.377 1.817 2.346 3885 1.001
## lambda_home[153] 1.129 0.003 0.163 0.827 1.121 1.469 2953 1.001
## lambda_home[154] 3.302 0.007 0.420 2.537 3.272 4.168 3845 1.001
## lambda_home[155] 1.144 0.003 0.162 0.859 1.136 1.481 3372 1.001
## lambda_home[156] 1.083 0.003 0.162 0.787 1.074 1.426 2863 1.000
## lambda_home[157] 3.263 0.008 0.472 2.407 3.226 4.251 3734 1.000
## lambda_home[158] 0.776 0.003 0.140 0.518 0.771 1.070 1752 1.002
## lambda_home[159] 1.654 0.005 0.248 1.193 1.641 2.144 2765 1.001
## lambda_home[160] 1.529 0.005 0.246 1.069 1.521 2.036 2338 1.001
## lambda_home[161] 0.908 0.003 0.141 0.646 0.903 1.197 2674 1.001
## lambda_home[162] 2.995 0.009 0.430 2.210 2.972 3.910 2086 1.002
## lambda_home[163] 1.617 0.004 0.239 1.173 1.609 2.107 3103 1.000
## lambda_home[164] 1.661 0.005 0.245 1.213 1.655 2.176 2795 1.000
## lambda_home[165] 1.985 0.006 0.311 1.408 1.972 2.632 2591 1.000
## lambda_home[166] 3.025 0.008 0.417 2.280 3.000 3.902 3063 0.999
## lambda_home[167] 2.902 0.008 0.441 2.151 2.859 3.882 2871 1.001
## lambda_home[168] 1.245 0.004 0.193 0.886 1.235 1.639 2785 1.000
## lambda_home[169] 2.157 0.005 0.304 1.602 2.142 2.790 3441 1.000
## lambda_home[170] 3.930 0.009 0.532 2.967 3.905 5.058 3215 1.000
## lambda_home[171] 1.024 0.004 0.173 0.699 1.018 1.373 1806 1.002
## lambda_home[172] 2.661 0.007 0.380 1.997 2.632 3.465 3329 1.000
## lambda_home[173] 1.266 0.004 0.193 0.907 1.260 1.665 2878 1.000
## lambda_home[174] 1.459 0.003 0.207 1.095 1.446 1.898 3777 1.000
## lambda_home[175] 1.007 0.004 0.170 0.689 1.005 1.353 1709 1.002
## lambda_home[176] 1.505 0.005 0.250 1.055 1.496 2.029 2322 1.001
## lambda_home[177] 1.792 0.005 0.276 1.298 1.774 2.379 3589 1.000
## lambda_home[178] 3.493 0.008 0.482 2.690 3.452 4.526 3430 1.000
## lambda_home[179] 0.831 0.003 0.144 0.567 0.826 1.126 1979 1.001
## lambda_home[180] 1.635 0.005 0.245 1.175 1.623 2.147 2894 1.001
## lambda_home[181] 1.772 0.004 0.269 1.289 1.761 2.337 3887 1.000
## lambda_home[182] 2.976 0.007 0.412 2.212 2.960 3.836 3245 1.000
## lambda_home[183] 0.903 0.003 0.143 0.635 0.897 1.190 2293 1.001
## lambda_home[184] 1.670 0.004 0.242 1.230 1.662 2.159 3313 1.000
## lambda_home[185] 1.627 0.005 0.245 1.175 1.620 2.114 2552 1.002
## lambda_home[186] 1.365 0.004 0.216 0.978 1.348 1.822 3788 1.000
## lambda_home[187] 0.893 0.003 0.140 0.633 0.887 1.180 2478 1.001
## lambda_home[188] 2.169 0.005 0.304 1.612 2.161 2.815 3340 1.000
## lambda_home[189] 0.938 0.005 0.178 0.611 0.932 1.298 1532 1.003
## lambda_home[190] 1.253 0.004 0.194 0.890 1.249 1.655 2508 1.001
## lambda_home[191] 2.644 0.006 0.369 2.025 2.610 3.470 3468 1.000
## lambda_home[192] 0.688 0.002 0.114 0.479 0.685 0.923 2437 1.001
## lambda_home[193] 1.018 0.004 0.175 0.685 1.011 1.372 1826 1.001
## lambda_home[194] 2.134 0.005 0.301 1.582 2.122 2.763 3507 1.000
## lambda_home[195] 1.852 0.007 0.290 1.315 1.847 2.429 1938 1.002
## lambda_home[196] 2.687 0.006 0.365 2.054 2.658 3.489 3321 1.000
## lambda_home[197] 1.802 0.004 0.277 1.314 1.778 2.399 3945 1.000
## lambda_home[198] 1.514 0.006 0.253 1.041 1.504 2.022 2068 1.002
## lambda_home[199] 2.592 0.007 0.421 1.830 2.556 3.509 3634 0.999
## lambda_home[200] 1.462 0.005 0.287 0.949 1.447 2.062 2845 1.001
## lambda_home[201] 2.805 0.009 0.465 1.919 2.791 3.721 2980 1.000
## lambda_home[202] 0.974 0.004 0.190 0.627 0.967 1.354 1823 1.001
## lambda_home[203] 2.027 0.006 0.372 1.378 2.003 2.824 3848 1.000
## lambda_home[204] 0.898 0.003 0.160 0.602 0.893 1.230 3301 1.000
## lambda_home[205] 2.266 0.006 0.372 1.607 2.244 3.054 3960 0.999
## lambda_home[206] 1.775 0.007 0.320 1.180 1.768 2.414 1866 1.000
## lambda_home[207] 1.263 0.004 0.217 0.874 1.255 1.725 3141 1.000
## lambda_home[208] 1.634 0.005 0.282 1.117 1.622 2.221 3151 1.000
## lambda_home[209] 3.589 0.010 0.582 2.606 3.524 4.882 3529 1.000
## lambda_home[210] 1.335 0.004 0.247 0.893 1.325 1.849 3717 1.000
## lambda_home[211] 3.876 0.010 0.598 2.789 3.855 5.094 3846 1.000
## lambda_home[212] 1.752 0.005 0.288 1.248 1.732 2.398 3888 1.000
## lambda_home[213] 1.895 0.007 0.356 1.249 1.879 2.633 2977 1.001
## lambda_home[214] 2.584 0.007 0.421 1.835 2.558 3.475 3333 0.999
## lambda_home[215] 0.846 0.003 0.171 0.535 0.838 1.208 2994 1.000
## lambda_home[216] 2.799 0.011 0.489 1.903 2.782 3.823 1857 1.000
## lambda_home[217] 0.781 0.002 0.151 0.517 0.770 1.113 3930 0.999
## lambda_home[218] 2.274 0.006 0.375 1.630 2.244 3.111 3989 1.000
## lambda_home[219] 1.747 0.005 0.289 1.238 1.731 2.372 3975 1.000
## lambda_home[220] 1.037 0.004 0.191 0.690 1.031 1.427 2547 1.000
## lambda_home[221] 1.847 0.005 0.330 1.266 1.820 2.575 3947 1.000
## lambda_home[222] 0.957 0.003 0.167 0.668 0.945 1.324 4085 0.999
## lambda_home[223] 2.764 0.008 0.435 2.014 2.723 3.731 3289 1.000
## lambda_home[224] 1.332 0.004 0.253 0.882 1.321 1.878 3469 1.000
## lambda_home[225] 2.987 0.008 0.463 2.164 2.966 3.982 3585 0.999
## lambda_home[226] 1.187 0.004 0.220 0.780 1.178 1.635 2615 1.001
## lambda_home[227] 1.514 0.004 0.253 1.087 1.490 2.100 3966 0.999
## lambda_home[228] 2.124 0.006 0.349 1.478 2.112 2.849 3748 1.000
## lambda_home[229] 1.643 0.005 0.278 1.136 1.627 2.228 3621 1.000
## lambda_home[230] 0.752 0.004 0.152 0.468 0.749 1.063 1863 1.001
## lambda_home[231] 2.019 0.006 0.368 1.378 1.995 2.803 3978 1.000
## lambda_home[232] 0.693 0.002 0.129 0.453 0.689 0.954 3110 1.000
## lambda_home[233] 1.534 0.006 0.281 1.003 1.528 2.080 2322 1.000
## lambda_home[234] 1.040 0.004 0.189 0.697 1.031 1.432 2036 1.000
## lambda_home[235] 1.423 0.004 0.253 0.983 1.400 1.990 3953 1.000
## lambda_away[1] 1.501 0.004 0.232 1.089 1.490 2.016 3516 1.001
## lambda_away[2] 1.341 0.003 0.216 0.973 1.322 1.819 3944 1.000
## lambda_away[3] 0.713 0.002 0.121 0.501 0.702 0.979 2393 1.001
## lambda_away[4] 1.548 0.004 0.235 1.132 1.531 2.058 3893 1.001
## lambda_away[5] 1.727 0.006 0.268 1.239 1.714 2.290 2037 1.000
## lambda_away[6] 1.308 0.004 0.204 0.944 1.294 1.750 2952 1.000
## lambda_away[7] 1.538 0.004 0.238 1.122 1.518 2.056 3919 1.001
## lambda_away[8] 0.623 0.002 0.114 0.435 0.610 0.884 2746 1.003
## lambda_away[9] 1.984 0.007 0.322 1.387 1.979 2.637 2288 1.002
## lambda_away[10] 1.317 0.003 0.203 0.968 1.299 1.763 3463 1.000
## lambda_away[11] 1.141 0.003 0.187 0.814 1.126 1.554 3350 1.001
## lambda_away[12] 2.343 0.006 0.349 1.738 2.318 3.129 3861 1.000
## lambda_away[13] 0.736 0.002 0.123 0.532 0.722 1.020 3033 1.002
## lambda_away[14] 1.020 0.004 0.190 0.711 1.000 1.449 1839 1.003
## lambda_away[15] 0.940 0.003 0.149 0.657 0.934 1.244 3140 0.999
## lambda_away[16] 1.277 0.004 0.205 0.929 1.257 1.727 3164 1.000
## lambda_away[17] 1.026 0.004 0.187 0.721 1.009 1.444 1845 1.002
## lambda_away[18] 2.401 0.006 0.362 1.767 2.381 3.173 3884 1.000
## lambda_away[19] 1.268 0.004 0.198 0.909 1.250 1.696 3190 0.999
## lambda_away[20] 0.996 0.005 0.190 0.692 0.972 1.425 1738 1.002
## lambda_away[21] 1.269 0.003 0.205 0.919 1.251 1.733 3712 1.000
## lambda_away[22] 1.142 0.003 0.191 0.813 1.126 1.567 3232 1.001
## lambda_away[23] 1.277 0.004 0.200 0.924 1.261 1.717 2628 0.999
## lambda_away[24] 2.401 0.006 0.362 1.778 2.366 3.163 3622 1.000
## lambda_away[25] 1.106 0.003 0.185 0.799 1.088 1.529 3602 1.001
## lambda_away[26] 0.995 0.004 0.180 0.697 0.976 1.411 1637 1.002
## lambda_away[27] 1.309 0.004 0.213 0.947 1.293 1.774 3684 1.000
## lambda_away[28] 1.114 0.003 0.183 0.795 1.096 1.519 2929 1.001
## lambda_away[29] 2.756 0.007 0.414 2.017 2.741 3.643 3846 1.001
## lambda_away[30] 0.870 0.004 0.174 0.588 0.850 1.271 1719 1.003
## lambda_away[31] 1.299 0.003 0.206 0.949 1.283 1.755 3871 1.000
## lambda_away[32] 1.501 0.004 0.236 1.090 1.480 2.034 3770 1.000
## lambda_away[33] 0.731 0.002 0.126 0.521 0.719 1.011 3476 1.002
## lambda_away[34] 1.686 0.005 0.264 1.191 1.678 2.224 2579 1.001
## lambda_away[35] 2.326 0.006 0.343 1.723 2.301 3.084 3884 1.000
## lambda_away[36] 1.824 0.007 0.306 1.324 1.788 2.498 2129 1.000
## lambda_away[37] 1.340 0.004 0.211 0.981 1.325 1.796 3605 1.001
## lambda_away[38] 0.714 0.002 0.127 0.507 0.699 0.991 3291 1.001
## lambda_away[39] 1.674 0.005 0.261 1.193 1.663 2.230 2809 1.001
## lambda_away[40] 1.349 0.003 0.210 0.988 1.333 1.805 3923 1.000
## lambda_away[41] 1.728 0.006 0.273 1.215 1.717 2.293 2287 1.001
## lambda_away[42] 1.349 0.003 0.205 0.992 1.332 1.797 3784 1.000
## lambda_away[43] 1.291 0.003 0.176 0.981 1.276 1.680 3907 1.000
## lambda_away[44] 1.068 0.003 0.182 0.752 1.055 1.457 3413 0.999
## lambda_away[45] 1.023 0.003 0.149 0.767 1.012 1.350 2653 1.000
## lambda_away[46] 1.591 0.004 0.229 1.192 1.572 2.083 3995 1.000
## lambda_away[47] 0.711 0.002 0.102 0.530 0.704 0.935 2058 1.002
## lambda_away[48] 1.301 0.003 0.193 0.948 1.290 1.706 3888 1.000
## lambda_away[49] 1.864 0.004 0.245 1.428 1.853 2.385 3645 1.000
## lambda_away[50] 1.299 0.004 0.179 0.968 1.287 1.676 2307 1.000
## lambda_away[51] 0.593 0.002 0.106 0.408 0.585 0.822 2046 1.000
## lambda_away[52] 1.017 0.002 0.146 0.763 1.007 1.337 3988 1.000
## lambda_away[53] 1.842 0.004 0.240 1.411 1.829 2.352 3389 1.001
## lambda_away[54] 2.347 0.005 0.302 1.804 2.329 2.990 3534 1.000
## lambda_away[55] 0.720 0.002 0.107 0.529 0.714 0.948 3877 1.000
## lambda_away[56] 0.563 0.002 0.087 0.406 0.555 0.755 2666 1.001
## lambda_away[57] 0.859 0.003 0.147 0.603 0.850 1.173 2897 1.000
## lambda_away[58] 1.832 0.004 0.237 1.409 1.816 2.346 2889 1.000
## lambda_away[59] 1.270 0.003 0.180 0.964 1.256 1.665 3090 1.001
## lambda_away[60] 1.028 0.003 0.152 0.766 1.014 1.358 2887 1.000
## lambda_away[61] 1.959 0.005 0.312 1.410 1.936 2.635 3339 1.000
## lambda_away[62] 1.046 0.003 0.153 0.775 1.037 1.366 3640 1.000
## lambda_away[63] 1.013 0.002 0.136 0.770 1.005 1.308 3737 1.000
## lambda_away[64] 1.058 0.004 0.189 0.731 1.046 1.472 2531 1.000
## lambda_away[65] 1.034 0.002 0.148 0.766 1.028 1.341 3638 1.001
## lambda_away[66] 1.263 0.003 0.177 0.944 1.251 1.640 3830 1.000
## lambda_away[67] 1.852 0.004 0.242 1.416 1.834 2.375 3627 1.000
## lambda_away[68] 1.970 0.005 0.311 1.419 1.947 2.624 3528 1.001
## lambda_away[69] 0.713 0.002 0.109 0.520 0.705 0.954 2468 1.001
## lambda_away[70] 1.277 0.003 0.174 0.963 1.267 1.640 4104 0.999
## lambda_away[71] 0.699 0.002 0.103 0.515 0.693 0.923 3128 1.000
## lambda_away[72] 3.564 0.014 0.567 2.606 3.513 4.809 1669 1.002
## lambda_away[73] 1.609 0.005 0.238 1.182 1.590 2.144 2447 1.001
## lambda_away[74] 0.707 0.002 0.101 0.531 0.701 0.925 3773 1.001
## lambda_away[75] 1.277 0.003 0.182 0.961 1.261 1.686 3329 1.000
## lambda_away[76] 1.950 0.006 0.316 1.402 1.923 2.657 2534 1.000
## lambda_away[77] 2.909 0.007 0.407 2.195 2.880 3.812 3694 1.000
## lambda_away[78] 1.034 0.002 0.148 0.774 1.022 1.351 3526 1.000
## lambda_away[79] 1.559 0.004 0.251 1.112 1.533 2.099 3706 1.000
## lambda_away[80] 1.278 0.003 0.173 0.975 1.262 1.654 3617 1.000
## lambda_away[81] 1.039 0.002 0.149 0.778 1.031 1.353 3958 1.000
## lambda_away[82] 1.615 0.004 0.222 1.210 1.600 2.095 3665 1.000
## lambda_away[83] 2.322 0.005 0.297 1.793 2.300 2.952 3009 1.001
## lambda_away[84] 1.274 0.004 0.201 0.938 1.255 1.727 3012 1.000
## lambda_away[85] 0.886 0.002 0.131 0.658 0.878 1.171 2945 1.000
## lambda_away[86] 1.270 0.003 0.172 0.960 1.262 1.634 4177 0.999
## lambda_away[87] 2.812 0.012 0.448 2.048 2.778 3.799 1417 1.002
## lambda_away[88] 0.896 0.002 0.128 0.669 0.885 1.164 2990 1.000
## lambda_away[89] 1.625 0.004 0.224 1.223 1.610 2.094 3315 1.000
## lambda_away[90] 0.703 0.002 0.100 0.526 0.697 0.916 3388 1.000
## lambda_away[91] 2.294 0.005 0.318 1.739 2.271 2.999 3873 1.000
## lambda_away[92] 2.938 0.006 0.390 2.246 2.909 3.751 3816 1.001
## lambda_away[93] 1.937 0.006 0.312 1.387 1.912 2.605 2930 1.000
## lambda_away[94] 0.901 0.003 0.128 0.668 0.894 1.178 1797 1.001
## lambda_away[95] 1.608 0.004 0.224 1.203 1.599 2.083 4050 0.999
## lambda_away[96] 1.300 0.003 0.176 0.991 1.290 1.679 3951 1.000
## lambda_away[97] 0.850 0.003 0.147 0.592 0.839 1.167 2432 1.000
## lambda_away[98] 1.583 0.004 0.241 1.163 1.561 2.122 3217 1.000
## lambda_away[99] 1.577 0.005 0.262 1.115 1.555 2.160 2614 1.001
## lambda_away[100] 1.627 0.004 0.213 1.252 1.612 2.096 3459 1.000
## lambda_away[101] 1.600 0.004 0.233 1.200 1.580 2.118 3810 1.000
## lambda_away[102] 0.891 0.002 0.126 0.668 0.883 1.167 3292 1.000
## lambda_away[103] 1.307 0.003 0.179 0.980 1.297 1.685 2725 1.000
## lambda_away[104] 0.844 0.002 0.144 0.593 0.835 1.151 3904 1.000
## lambda_away[105] 2.318 0.005 0.312 1.764 2.299 2.981 3793 1.001
## lambda_away[106] 1.292 0.004 0.179 0.970 1.282 1.669 2253 1.000
## lambda_away[107] 2.362 0.005 0.304 1.816 2.347 3.009 3915 1.000
## lambda_away[108] 1.028 0.002 0.148 0.764 1.019 1.337 3837 1.000
## lambda_away[109] 0.468 0.002 0.084 0.322 0.461 0.649 2529 1.000
## lambda_away[110] 2.335 0.005 0.300 1.803 2.318 2.974 3860 1.001
## lambda_away[111] 1.293 0.003 0.176 0.976 1.281 1.668 3898 1.000
## lambda_away[112] 0.569 0.001 0.086 0.415 0.563 0.758 3611 1.000
## lambda_away[113] 0.854 0.003 0.149 0.589 0.843 1.179 2705 1.000
## lambda_away[114] 1.597 0.004 0.224 1.192 1.582 2.073 3802 1.000
## lambda_away[115] 2.257 0.006 0.303 1.696 2.240 2.913 2944 1.001
## lambda_away[116] 1.233 0.003 0.178 0.913 1.222 1.615 3486 1.000
## lambda_away[117] 0.988 0.002 0.143 0.736 0.981 1.293 3956 1.000
## lambda_away[118] 1.216 0.003 0.177 0.904 1.205 1.584 3647 1.000
## lambda_away[119] 0.668 0.002 0.098 0.485 0.665 0.879 3282 1.000
## lambda_away[120] 2.811 0.007 0.378 2.121 2.790 3.604 3133 1.001
## lambda_away[121] 1.286 0.003 0.187 0.955 1.276 1.684 3371 1.000
## lambda_away[122] 0.678 0.002 0.099 0.500 0.670 0.888 3735 0.999
## lambda_away[123] 0.990 0.002 0.146 0.742 0.980 1.301 3572 0.999
## lambda_away[124] 1.212 0.003 0.172 0.887 1.208 1.573 3929 1.000
## lambda_away[125] 0.803 0.002 0.122 0.585 0.796 1.055 2985 1.000
## lambda_away[126] 0.641 0.002 0.097 0.464 0.637 0.843 3383 1.000
## lambda_away[127] 1.549 0.004 0.217 1.151 1.541 2.019 3200 1.000
## lambda_away[128] 1.282 0.003 0.182 0.948 1.276 1.668 3822 1.000
## lambda_away[129] 1.025 0.002 0.141 0.772 1.015 1.318 3866 1.000
## lambda_away[130] 1.466 0.004 0.213 1.069 1.460 1.909 3434 1.001
## lambda_away[131] 0.801 0.002 0.119 0.574 0.798 1.046 3133 1.000
## lambda_away[132] 1.300 0.003 0.183 0.968 1.288 1.684 3749 1.000
## lambda_away[133] 2.345 0.005 0.313 1.786 2.323 2.997 3866 1.002
## lambda_away[134] 0.812 0.002 0.118 0.597 0.805 1.054 3512 1.000
## lambda_away[135] 0.769 0.002 0.118 0.557 0.763 1.016 3055 1.001
## lambda_away[136] 0.670 0.002 0.102 0.486 0.664 0.883 2417 0.999
## lambda_away[137] 1.227 0.003 0.173 0.901 1.221 1.583 3807 1.000
## lambda_away[138] 1.532 0.004 0.221 1.130 1.524 1.988 3181 1.001
## lambda_away[139] 1.244 0.003 0.175 0.928 1.236 1.618 3817 1.000
## lambda_away[140] 2.261 0.005 0.297 1.721 2.246 2.914 3795 1.001
## lambda_away[141] 1.177 0.003 0.169 0.865 1.168 1.527 3689 1.001
## lambda_away[142] 0.538 0.001 0.082 0.392 0.531 0.715 3317 0.999
## lambda_away[143] 2.230 0.005 0.298 1.688 2.218 2.869 3867 1.000
## lambda_away[144] 1.529 0.004 0.219 1.103 1.522 1.969 3069 1.000
## lambda_away[145] 1.883 0.004 0.255 1.428 1.867 2.437 3542 1.001
## lambda_away[146] 1.247 0.003 0.174 0.937 1.235 1.616 4014 1.000
## lambda_away[147] 0.976 0.002 0.141 0.725 0.970 1.275 3834 0.999
## lambda_away[148] 2.359 0.005 0.318 1.796 2.345 3.053 3611 1.000
## lambda_away[149] 1.180 0.003 0.169 0.874 1.171 1.532 3916 1.000
## lambda_away[150] 1.473 0.003 0.200 1.100 1.465 1.894 3566 1.000
## lambda_away[151] 1.032 0.002 0.152 0.771 1.017 1.365 3982 0.999
## lambda_away[152] 1.163 0.003 0.169 0.860 1.156 1.535 3774 1.000
## lambda_away[153] 1.887 0.004 0.256 1.444 1.871 2.447 3649 1.001
## lambda_away[154] 0.645 0.002 0.097 0.472 0.639 0.852 3688 1.000
## lambda_away[155] 1.860 0.004 0.250 1.412 1.846 2.382 4120 1.001
## lambda_away[156] 1.968 0.005 0.273 1.505 1.948 2.567 3478 1.001
## lambda_away[157] 0.607 0.002 0.104 0.415 0.604 0.824 2320 1.001
## lambda_away[158] 2.555 0.007 0.374 1.941 2.515 3.405 3031 1.001
## lambda_away[159] 1.195 0.003 0.183 0.860 1.190 1.584 3488 1.001
## lambda_away[160] 1.295 0.003 0.201 0.947 1.277 1.732 3768 0.999
## lambda_away[161] 2.175 0.005 0.304 1.606 2.157 2.823 4041 1.000
## lambda_away[162] 0.660 0.002 0.107 0.460 0.654 0.885 2355 1.001
## lambda_away[163] 1.221 0.003 0.181 0.888 1.210 1.592 3781 1.000
## lambda_away[164] 1.188 0.003 0.176 0.867 1.182 1.564 3273 1.001
## lambda_away[165] 0.997 0.003 0.160 0.714 0.984 1.350 3877 0.999
## lambda_away[166] 0.652 0.002 0.101 0.462 0.650 0.860 2935 1.000
## lambda_away[167] 0.685 0.003 0.127 0.449 0.681 0.947 1775 1.002
## lambda_away[168] 1.586 0.004 0.226 1.173 1.575 2.068 3411 1.000
## lambda_away[169] 0.915 0.002 0.140 0.660 0.911 1.211 3432 1.000
## lambda_away[170] 0.503 0.002 0.083 0.350 0.499 0.676 3013 1.001
## lambda_away[171] 1.934 0.005 0.282 1.452 1.912 2.532 3383 1.000
## lambda_away[172] 0.744 0.003 0.127 0.506 0.739 1.004 1996 1.001
## lambda_away[173] 1.559 0.004 0.222 1.160 1.547 2.021 3595 0.999
## lambda_away[174] 1.352 0.004 0.207 0.965 1.346 1.780 2485 1.001
## lambda_away[175] 1.966 0.005 0.284 1.480 1.942 2.596 3115 1.001
## lambda_away[176] 1.317 0.003 0.207 0.961 1.302 1.773 3907 1.000
## lambda_away[177] 1.106 0.004 0.183 0.774 1.100 1.484 2520 1.001
## lambda_away[178] 0.567 0.002 0.100 0.383 0.562 0.773 2196 1.001
## lambda_away[179] 2.385 0.006 0.357 1.769 2.353 3.191 3739 1.000
## lambda_away[180] 1.208 0.003 0.180 0.882 1.200 1.570 3381 1.001
## lambda_away[181] 1.117 0.003 0.181 0.782 1.108 1.498 2978 1.001
## lambda_away[182] 0.663 0.002 0.103 0.472 0.660 0.878 2967 1.000
## lambda_away[183] 2.188 0.006 0.313 1.631 2.166 2.849 3057 1.001
## lambda_away[184] 1.181 0.003 0.175 0.864 1.171 1.544 3772 1.000
## lambda_away[185] 1.215 0.003 0.185 0.888 1.207 1.611 2908 1.001
## lambda_away[186] 1.451 0.004 0.228 1.027 1.443 1.918 3079 1.000
## lambda_away[187] 2.211 0.005 0.313 1.648 2.196 2.864 3975 1.000
## lambda_away[188] 0.910 0.002 0.141 0.654 0.899 1.204 3446 1.000
## lambda_away[189] 2.123 0.007 0.341 1.559 2.085 2.890 2557 1.001
## lambda_away[190] 1.577 0.004 0.230 1.161 1.567 2.062 3206 1.000
## lambda_away[191] 0.748 0.003 0.124 0.517 0.743 1.005 2101 1.001
## lambda_away[192] 2.873 0.007 0.399 2.167 2.845 3.768 3556 1.000
## lambda_away[193] 1.946 0.005 0.286 1.450 1.921 2.591 3387 1.000
## lambda_away[194] 0.925 0.002 0.142 0.664 0.916 1.226 3676 1.000
## lambda_away[195] 1.067 0.003 0.157 0.789 1.057 1.398 3448 1.000
## lambda_away[196] 0.735 0.003 0.121 0.508 0.729 0.978 1998 1.001
## lambda_away[197] 1.099 0.003 0.181 0.765 1.093 1.481 2976 1.001
## lambda_away[198] 1.310 0.004 0.210 0.945 1.292 1.771 3391 1.000
## lambda_away[199] 0.741 0.003 0.140 0.480 0.735 1.040 2404 1.000
## lambda_away[200] 1.321 0.004 0.251 0.878 1.302 1.870 4100 1.000
## lambda_away[201] 0.684 0.002 0.123 0.467 0.675 0.960 4170 0.999
## lambda_away[202] 1.976 0.005 0.327 1.417 1.949 2.711 3770 1.000
## lambda_away[203] 0.952 0.003 0.188 0.611 0.944 1.350 3791 0.999
## lambda_away[204] 2.135 0.006 0.349 1.490 2.118 2.872 3923 1.000
## lambda_away[205] 0.848 0.003 0.160 0.551 0.842 1.175 2758 1.001
## lambda_away[206] 1.081 0.003 0.184 0.776 1.065 1.501 4074 1.000
## lambda_away[207] 1.517 0.004 0.254 1.055 1.507 2.062 3512 1.000
## lambda_away[208] 1.174 0.003 0.206 0.804 1.165 1.605 3695 1.001
## lambda_away[209] 0.537 0.003 0.110 0.332 0.533 0.762 1868 1.001
## lambda_away[210] 1.441 0.004 0.261 0.978 1.422 1.982 3553 1.000
## lambda_away[211] 0.495 0.002 0.093 0.326 0.490 0.690 3112 1.000
## lambda_away[212] 1.095 0.004 0.200 0.727 1.092 1.508 2477 1.000
## lambda_away[213] 1.018 0.003 0.194 0.682 1.001 1.464 4184 1.000
## lambda_away[214] 0.743 0.003 0.140 0.488 0.738 1.031 2181 1.000
## lambda_away[215] 2.282 0.006 0.403 1.556 2.250 3.152 4042 1.000
## lambda_away[216] 0.687 0.002 0.126 0.457 0.678 0.968 3895 1.000
## lambda_away[217] 2.469 0.008 0.444 1.655 2.459 3.389 3467 1.000
## lambda_away[218] 0.845 0.003 0.160 0.550 0.838 1.174 2450 1.000
## lambda_away[219] 1.098 0.004 0.200 0.729 1.091 1.503 2903 1.001
## lambda_away[220] 1.851 0.005 0.306 1.302 1.828 2.500 3890 1.000
## lambda_away[221] 1.043 0.004 0.201 0.669 1.034 1.461 2759 1.001
## lambda_away[222] 2.004 0.006 0.341 1.343 1.998 2.705 3525 1.000
## lambda_away[223] 0.695 0.003 0.137 0.443 0.690 0.970 1795 1.001
## lambda_away[224] 1.447 0.004 0.266 0.974 1.428 2.016 3576 1.000
## lambda_away[225] 0.641 0.002 0.114 0.430 0.638 0.880 3116 1.000
## lambda_away[226] 1.620 0.005 0.279 1.120 1.603 2.222 3563 0.999
## lambda_away[227] 1.268 0.005 0.235 0.824 1.264 1.739 2209 1.000
## lambda_away[228] 0.903 0.003 0.161 0.609 0.892 1.250 3321 1.000
## lambda_away[229] 1.167 0.004 0.205 0.794 1.155 1.584 3131 0.999
## lambda_away[230] 2.564 0.007 0.430 1.839 2.519 3.563 3673 1.001
## lambda_away[231] 0.955 0.003 0.185 0.613 0.941 1.343 3962 1.000
## lambda_away[232] 2.770 0.007 0.448 1.951 2.751 3.717 3790 1.000
## lambda_away[233] 1.253 0.004 0.220 0.869 1.237 1.734 3491 1.000
## lambda_away[234] 1.845 0.005 0.303 1.291 1.824 2.497 3987 0.999
## lambda_away[235] 1.351 0.005 0.248 0.895 1.337 1.871 2843 1.001
## lp__ -209.853 0.763 18.612 -240.026 -212.424 -168.641 595 1.007
##
## Samples were drawn using NUTS(diag_e) at Tue Jan 23 16:06:26 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).
season_sigma
パラメータのトレースプロットと分布を示している。
plt_g <- stan_trace(fit3,
pars = 'season_sigma',
separate_chains = TRUE) + ggtitle('Traceplot of season_sigma')
plt_h <- stan_dens(fit3,
pars = 'season_sigma',
separate_chains = TRUE)
(plt_g | plt_h)
この部分は参考先のブログにないが、各チームのスキルの時系列の推移を可視化しておく。
ms3 <- rstan::extract(fit3)
d_est <- data.frame()
for (n in 1:length(t)) {
qua <- apply(ms3$skill[,n,], 2, quantile, prob = c(0.25, 0.5, 0.75))
d_est <- rbind(
d_est,
data.frame(team = t[n], year = 2008:2012, t(qua), check.names = FALSE)
)
}
d_est
## team year 25% 50% 75%
## 1 Real Valladolid 2008 -0.26249083 -0.1615685632 -0.05934306
## 2 Real Valladolid 2009 -0.32337130 -0.2249560375 -0.12417559
## 3 Real Valladolid 2010 -0.32043656 -0.2101265006 -0.10693089
## 4 Real Valladolid 2011 -0.31286318 -0.1969365122 -0.07746688
## 5 Real Valladolid 2012 -0.30868639 -0.1857577812 -0.06633374
## 6 Atlético Madrid 2008 -0.12379327 -0.0241530309 0.07641963
## 7 Atlético Madrid 2009 -0.09880099 -0.0031711579 0.09873116
## 8 Atlético Madrid 2010 -0.11832397 -0.0175747215 0.08387438
## 9 Atlético Madrid 2011 -0.11358870 -0.0091338154 0.09300948
## 10 Atlético Madrid 2012 -0.08970033 0.0223445652 0.13654988
## 11 FC Barcelona 2008 0.48162524 0.5799015219 0.68305575
## 12 FC Barcelona 2009 0.49190267 0.5889136526 0.68709019
## 13 FC Barcelona 2010 0.49156226 0.5871710590 0.69009676
## 14 FC Barcelona 2011 0.48233123 0.5894673868 0.69308181
## 15 FC Barcelona 2012 0.45090095 0.5567282833 0.66815718
## 16 Villarreal CF 2008 -0.12644633 -0.0254321826 0.07703012
## 17 Villarreal CF 2009 -0.11916071 -0.0205209857 0.08242666
## 18 Villarreal CF 2010 -0.16390142 -0.0629600117 0.03975689
## 19 Villarreal CF 2011 -0.20675863 -0.0954075940 0.01733701
## 20 Villarreal CF 2012 -0.21868405 -0.0924018560 0.03199746
## 21 FC Sevilla 2008 -0.09488935 -0.0006952173 0.10136193
## 22 FC Sevilla 2009 -0.10318872 -0.0071623203 0.09127598
## 23 FC Sevilla 2010 -0.10370095 -0.0077911016 0.09668257
## 24 FC Sevilla 2011 -0.12354466 -0.0195693763 0.08418438
## 25 FC Sevilla 2012 -0.15230075 -0.0429657293 0.06702181
## 26 Real Madrid CF 2008 0.14547865 0.2551388318 0.36205160
## 27 Real Madrid CF 2009 0.25647275 0.3535247693 0.44608168
## 28 Real Madrid CF 2010 0.30931737 0.4060032426 0.50742830
## 29 Real Madrid CF 2011 0.36050070 0.4685139643 0.58044966
## 30 Real Madrid CF 2012 0.36204427 0.4752129620 0.59292867
## 31 FC Valencia 2008 -0.09090930 0.0044240489 0.10325470
## 32 FC Valencia 2009 -0.10585777 -0.0109662292 0.08337908
## 33 FC Valencia 2010 -0.11749729 -0.0174772690 0.08069139
## 34 FC Valencia 2011 -0.10741789 -0.0071713985 0.09902800
## 35 FC Valencia 2012 -0.09384637 0.0172415133 0.12990448
## 36 Real Zaragoza 2008 -0.32193538 -0.2158215783 -0.11037751
## 37 Real Zaragoza 2009 -0.32775026 -0.2345322083 -0.13833023
## 38 Real Zaragoza 2010 -0.33571247 -0.2363916555 -0.13576214
## 39 Real Zaragoza 2011 -0.38773211 -0.2799226997 -0.17807933
## 40 Real Zaragoza 2012 -0.41428386 -0.3046215089 -0.19089745
## 41 CD Tenerife 2008 -0.52071400 -0.3909378123 -0.26386697
## 42 CD Tenerife 2009 -0.53980583 -0.4216834732 -0.30852706
## 43 CD Tenerife 2010 -0.55187257 -0.4206950800 -0.29589022
## 44 CD Tenerife 2011 -0.55980247 -0.4228324084 -0.28183474
## 45 CD Tenerife 2012 -0.57440452 -0.4223631739 -0.26829503
バルセロナ、レアル・マドリード、アトレチコ・マドリードの3チームに色付けをして可視化する。緑はアトレチコ・マドリード、赤はレアル・マドリード、青がバルセロナである。今回のデータセットの時期における実績は下記の通りで、圧倒的にバルセロナ優位であり、可視化した結果とも一致する。レアルマドリードについては、2011-12に優勝するまでに徐々に力をつけていたことがわかる。
col_team <- c('darkgreen', 'gray', 'cornflowerblue', 'gray', 'gray', 'coral', 'gray', 'gray', 'gray')
ggplot(data = d_est, aes(x = year, y = `50%`, group = team)) +
theme_bw(base_size = 12) +
geom_ribbon(aes(ymin = `25%`, ymax = `75%`, fill = team), alpha = 0.2) +
geom_line(aes(col = team), linewidth = 0.5) +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(x = 'Year', y = 'Strength', title = 'Time Series of Estimated Latent Strength') +
scale_x_continuous(breaks = 2008:2012) +
scale_y_continuous(breaks = seq(-1,1,1), limit = c(-1, 1)) +
scale_color_manual(values = col_team) +
scale_fill_manual(values = col_team)
## Ranking the Teams of La Liga
2012/2013シーズンから推定されたスキル・パラメタを用いて、リーガ・エスパニョーラのチームをランキングする。ブログには下記の記載があるが、Stanで実行しているモデルではスキルを0に固定せず、特定の分布から生成されるゆるい仮定のもとで推定しているので、そのまま推定されたスキルを利用する。
スキル・パラメタの値は、スキル・パラメタをゼロに「固定」したチームのスキルとの相対値なので、解釈が難しい。より解釈しやすい尺度に置き換えるために、まず全チームの平均スキルを差し引くことでスキルパラメータをゼロセンターにし、次にホームベースラインを加え、その結果を指数化。この再スケーリングされたスキルパラメーターは、ホームチームと対戦したときの予想ゴール数のスケールとなる。
tmp <- paste0('skill[', 1:length(teams), ',5]')
season1213_nm <- names(sort(get_posterior_mean(fit3, par = tmp)[, 'mean-all chains'], decreasing = TRUE))
stan_plot(
fit3,
pars = season1213_nm,
show_density = TRUE,
ci_level = 0.8,
outer_level = 1,
)
このままでは見ずらいが、チーム名は下記の通りで並んでいる。
teams[as.numeric(sub(".*\\[([0-9]+).*", "\\1", season1213_nm))]
## [1] "FC Barcelona" "Real Madrid CF" "Atlético Madrid" "FC Valencia"
## [5] "FC Sevilla" "Villarreal CF" "Real Valladolid" "Real Zaragoza"
## [9] "CD Tenerife"
FCバルセロナは、レアルマドリードに対して、71.05%の確率でより良いチームである。
plotPost(
paramSampleVec = ms3$skill[,3,5] - ms3$skill[,6,5],
compVal = 0,
xlab = "← Real Madrid vs Barcelona →"
)
laliga
のデータセットでは、2012/2013シーズンの最後の50試合の結果が欠落していた。モデルを用いて、これらの50試合の結果を予測し、シミュレーション。以下のRコードは、各試合(結果がわかっている試合とわかっていない試合の両方)について、いくつかの値を計算する。
dd <- laliga %>% filter(HomeTeam %in% t & AwayTeam %in% t)
n <- nrow(ms3$lp__)
m3_pred <- sapply(1:nrow(dd), function(i) {
home_team <- which(teams == dd$HomeTeam[1])
away_team <- which(teams == dd$AwayTeam[1])
season <- which(seasons == dd$Season[1])
home_skill <- ms3$skill[, home_team, season]
away_skill <- ms3$skill[, away_team, season]
home_baseline <- ms3$home_baseline[, season]
away_baseline <- ms3$away_baseline[, season]
home_goals <- rpois(n, exp(home_baseline + home_skill - away_skill))
away_goals <- rpois(n, exp(away_baseline + away_skill - home_skill))
home_goals_table <- table(home_goals)
away_goals_table <- table(away_goals)
match_results <- sign(home_goals - away_goals)
match_results_table <- table(match_results)
mode_home_goal <- as.numeric(names(home_goals_table)[which.max(home_goals_table)])
mode_away_goal <- as.numeric(names(away_goals_table)[which.max(away_goals_table)])
match_result <- as.numeric(names(match_results_table)[which.max(match_results_table)])
rand_i <- sample(seq_along(home_goals), 1)
c(
mode_home_goal = mode_home_goal,
mode_away_goal = mode_away_goal,
match_result = match_result,
mean_home_goal = mean(home_goals),
mean_away_goal = mean(away_goals),
rand_home_goal = home_goals[rand_i],
rand_away_goal = away_goals[rand_i],
rand_match_result = match_results[rand_i]
)
})
m3_pred <- t(m3_pred)
head(m3_pred, 30)
## mode_home_goal mode_away_goal match_result mean_home_goal mean_away_goal
## [1,] 1 1 1 1.64725 1.49800
## [2,] 1 1 1 1.69300 1.48000
## [3,] 1 1 1 1.63300 1.47375
## [4,] 1 1 1 1.60175 1.48050
## [5,] 1 1 1 1.66300 1.47350
## [6,] 1 1 1 1.64075 1.50950
## [7,] 1 1 1 1.62525 1.50275
## [8,] 1 1 1 1.59325 1.49575
## [9,] 1 1 1 1.64275 1.52100
## [10,] 1 1 1 1.65325 1.48200
## [11,] 1 1 1 1.65225 1.52575
## [12,] 1 1 1 1.62375 1.48800
## [13,] 1 1 1 1.64625 1.52400
## [14,] 1 1 1 1.66550 1.50800
## [15,] 1 1 1 1.66000 1.51700
## [16,] 1 1 1 1.63375 1.52525
## [17,] 1 1 1 1.63550 1.50850
## [18,] 1 1 1 1.62850 1.51600
## [19,] 1 1 1 1.63700 1.52525
## [20,] 1 1 1 1.63050 1.51000
## [21,] 1 1 1 1.63975 1.51475
## [22,] 1 1 1 1.65000 1.50325
## [23,] 1 1 1 1.65350 1.49975
## [24,] 1 1 1 1.62175 1.49150
## [25,] 1 1 1 1.64250 1.51050
## [26,] 1 1 1 1.62575 1.52600
## [27,] 1 1 1 1.60825 1.48425
## [28,] 1 1 1 1.61200 1.47000
## [29,] 1 1 1 1.61700 1.51800
## [30,] 1 1 1 1.59550 1.47800
## rand_home_goal rand_away_goal rand_match_result
## [1,] 1 3 -1
## [2,] 1 1 0
## [3,] 2 0 1
## [4,] 1 2 -1
## [5,] 0 1 -1
## [6,] 3 0 1
## [7,] 5 0 1
## [8,] 2 3 -1
## [9,] 0 0 0
## [10,] 0 1 -1
## [11,] 2 2 0
## [12,] 1 3 -1
## [13,] 2 0 1
## [14,] 2 2 0
## [15,] 2 3 -1
## [16,] 2 2 0
## [17,] 1 2 -1
## [18,] 4 1 1
## [19,] 0 2 -1
## [20,] 3 2 1
## [21,] 0 2 -1
## [22,] 1 1 0
## [23,] 3 1 1
## [24,] 0 6 -1
## [25,] 1 1 0
## [26,] 1 0 1
## [27,] 2 2 0
## [28,] 2 1 1
## [29,] 0 1 -1
## [30,] 0 4 -1
データ中のゴール数の分布を、全試合の予測最頻値、平均値、ランダムゴール数を比較してみる。 まず、ホームチームの観測された実際のゴール数分布。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = dd, aes(HomeGoals), fill = 'lightsteelblue', alpha = 1/2) +
labs(title = 'Distribution of home goal observations')
ほとんどすべての試合で、最も可能性の高いゴール数は1つ。実際、リーガ・エスパニョーラの試合について何も知らない場合、ホームチームの1ゴールに賭けることは、100 %の確率でベストベットです。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = data.frame(m3_pred), aes(mode_home_goal), fill = 'darkgreen', alpha = 1/2) +
scale_x_continuous(breaks = seq(0, 10, 1)) +
labs(title = 'Distribution of simlation Mode home goal')
ほとんどの試合で、予想ゴール数は1.6前後。つまり、最も安全なベットが1ゴールだったとしても、2ゴール前後が予想される。
ggplot() +
theme_bw(base_size = 15) +
geom_histogram(data = data.frame(m3_pred), aes(mean_home_goal), fill = 'coral', alpha = 1/2, bins = 25) +
labs(title = 'Distribution of simlation Mean home goal')
最頻値と平均ゴール数の分布は、実際のゴール数とは似ても似つかない。しかし、ランダム化されたゴール数(各試合のゴール数がその試合のホームゴール数分布からランダムに抽出されたもの)の分布は、実際のホームゴール数に似ていると予想される。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = data.frame(m3_pred), aes(rand_home_goal), fill = 'cornflowerblue', alpha = 1/2) +
scale_x_continuous(breaks = seq(0,10,1)) +
labs(title = 'Distribution of simlation Random home goal')
また、モデルがどの程度データを予測するかを見ることもできる。これはおそらくクロスバリデーションを用いて行うべきですが、有効なパラメータの数はデータポイントの数よりはるかに少ないので、直接比較することで少なくとも適切な範囲の予測精度を推定することができる。
list(
mean(laliga$HomeGoals == m3_pred[ , "mode_home_goal"], na.rm = TRUE),
mean((laliga$HomeGoals - m3_pred[ , "mean_home_goal"])^2, na.rm = TRUE)
)
## [[1]]
## [1] 0.3259459
##
## [[2]]
## [1] 1.816843
つまり平均して、モデルは正しいホームゴール数を 33 %の確率で予測し、平均二乗誤差 1.82 の平均ゴール数を推測する。
次に、実際の試合結果と予測された試合結果を見てみる。下のグラフは、1がホームチームの勝利、0が引き分け、-1がアウェイチームの勝利で、データの試合結果を示している。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = dd, aes(MatchResult), fill = 'darkgray', alpha = 1/2) +
scale_x_continuous(breaks = seq(-1, 1, 1)) +
labs(title = 'Distribution of MatchResult')
ほとんどすべての試合において、最も安全な賭けはホームチームに賭けることである。引き分けは珍しくないが、最も安全な賭けであることはない。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = data.frame(m3_pred), aes(match_result), fill = 'darkgray', alpha = 1/2) +
scale_x_continuous(breaks = seq(-1, 1, 1)) +
labs(title = 'Distribution of MatchResult')
ホームゴール数の場合と同様に、ランダム化された試合結果は、実際の試合結果と似た分布を持つ(データを限定しているので、あまり似ていないが)。
ggplot() +
theme_bw(base_size = 15) +
geom_bar(data = data.frame(m3_pred), aes(rand_match_result), fill = 'darkgray', alpha = 1/2) +
scale_x_continuous(breaks = seq(-1, 1, 1)) +
labs(title = 'Distribution of Simulation MatchResult')
このモデルは、正しい試合結果を 0.501%の確率で予測している。
mean(laliga$MatchResult == m3_pred[ , "match_result"], na.rm = TRUE)
## [1] 0.5005405
さて、モデルがリーガ・エスパニョーラの歴史を合理的に予測することをチェックしたので、リーガ・エスパニョーラの終盤戦を予測する。以下のコードは、終盤戦の予測されたモードと平均ゴール数、そして各ゲームの予測された勝者を表示する。
laliga_forecast <- dd[is.na(dd$HomeGoals), c("Season", "Week", "HomeTeam", "AwayTeam")]
m3_forecast <- m3_pred[is.na(dd$HomeGoals),]
laliga_forecast$mean_home_goals <- round(m3_forecast[,"mean_home_goal"], 1)
laliga_forecast$mean_away_goals <- round(m3_forecast[,"mean_away_goal"], 1)
laliga_forecast$mode_home_goals <- m3_forecast[,"mode_home_goal"]
laliga_forecast$mode_away_goals <- m3_forecast[,"mode_away_goal"]
laliga_forecast$predicted_winner <- ifelse(m3_forecast[ , "match_result"] == 1, laliga_forecast$HomeTeam,
ifelse(m3_forecast[ , "match_result"] == -1, laliga_forecast$AwayTeam, "Draw"))
rownames(laliga_forecast) <- NULL
print(xtable::xtable(laliga_forecast, align="cccccccccc"), type="html")
Season | Week | HomeTeam | AwayTeam | mean_home_goals | mean_away_goals | mode_home_goals | mode_away_goals | predicted_winner | |
---|---|---|---|---|---|---|---|---|---|
1 | 2012/13 | 34 | Real Madrid CF | Real Valladolid | 1.60 | 1.50 | 1.00 | 1.00 | Real Madrid CF |
2 | 2012/13 | 35 | Atlético Madrid | FC Barcelona | 1.60 | 1.50 | 1.00 | 1.00 | Atlético Madrid |
3 | 2012/13 | 36 | FC Barcelona | Real Valladolid | 1.60 | 1.50 | 1.00 | 1.00 | FC Barcelona |
4 | 2012/13 | 38 | FC Sevilla | FC Valencia | 1.70 | 1.50 | 1.00 | 1.00 | FC Sevilla |
5 | 2012/13 | 38 | Real Zaragoza | Atlético Madrid | 1.60 | 1.50 | 1.00 | 1.00 | Real Zaragoza |
これらの予測は、勝者になる可能性の高いチームに賭けるには良いが、実際の終盤戦がどのように展開するかは反映していない。
m3_pred
のmatch_result
はhome_goals - away_goals
の差をsign()
関数で変換し、ホームゴールが多ければホームチームの勝利(=1)、アウェイチームのゴールが多ければアウェイチームの勝利(=-1)として、MCMCのサンプル分シュミレーションした結果を集計し、最も頻度が多いものがmatch_result
となる。
最後に、先に計算したシミュレーション結果を表示して、リーガ・エスパニョーラ終盤戦の可能なバージョンを見る。つまり、各チームのスキルをもとにしてポアソン分布で発生させたゴール数の試合結果を見る。
laliga_sim <- dd[is.na(dd$HomeGoals), c("Season", "Week", "HomeTeam", "AwayTeam")]
laliga_sim$home_goals <- m3_forecast[,"rand_home_goal"]
laliga_sim$away_goals <- m3_forecast[,"rand_away_goal"]
laliga_sim$winner <- ifelse(m3_forecast[ , "rand_match_result"] == 1, laliga_forecast$HomeTeam,
ifelse(m3_forecast[ , "rand_match_result"] == -1, laliga_forecast$AwayTeam, "Draw"))
rownames(laliga_sim) <- NULL
print(xtable::xtable(laliga_sim, align="cccccccc"), type="html")
Season | Week | HomeTeam | AwayTeam | home_goals | away_goals | winner | |
---|---|---|---|---|---|---|---|
1 | 2012/13 | 34 | Real Madrid CF | Real Valladolid | 3.00 | 4.00 | Real Valladolid |
2 | 2012/13 | 35 | Atlético Madrid | FC Barcelona | 1.00 | 3.00 | FC Barcelona |
3 | 2012/13 | 36 | FC Barcelona | Real Valladolid | 4.00 | 0.00 | FC Barcelona |
4 | 2012/13 | 38 | FC Sevilla | FC Valencia | 4.00 | 0.00 | FC Sevilla |
5 | 2012/13 | 38 | Real Zaragoza | Atlético Madrid | 2.00 | 2.00 | Draw |
ベイズ・モデリングとMCMCサンプリングを使用することのパワーの1つは、いったんパラメータのMCMCサンプルが得られれば、パラメタ推定値の不確実性を保持したまま、これらの推定値から得られる任意の量を計算するのが簡単だということ。それでは、セビージャ対バレンシアの試合のゴール数の予測分布を計算するためにMCMCサンプルを使うことから始める。
n <- nrow(ms3$lp__)
home_team <- which(teams == "FC Sevilla")
away_team <- which(teams == "FC Valencia")
season <- which(seasons == "2012/13")
home_skill <- ms3$skill[, home_team, season]
away_skill <- ms3$skill[, away_team, season]
home_baseline <- ms3$home_baseline[, season]
away_baseline <- ms3$away_baseline[, season]
home_goals <- rpois(n, exp(home_baseline + home_skill - away_skill))
away_goals <- rpois(n, exp(away_baseline + away_skill - home_skill))
この2つの分布の要約を見ると、接戦になるだろうが、ホームチームのセビージャがわずかに有利であることがわかる。
old_par <- par(mfrow = c(2, 2), mar=rep(2.2, 4))
plot_goals(home_goals, away_goals)
par(old_par)
この記事を書いている時点(2013-05-28)では、ベッティングサイトwww.betsson.comでこのゲームの結果に賭けた場合、以下のようなペイアウト(つまり、ベットが成功した場合、いくら戻ってくるか)が得られます:
Sevilla Draw Valencia
3.2 3.35 2.15
シミュレーションしたゴール数分布を使って、私のモデルの予測ペイアウトを計算することができる。
1 / c(
Sevilla = mean(home_goals > away_goals),
Draw = mean(home_goals == away_goals),
Valencia = mean(home_goals < away_goals)
)
## Sevilla Draw Valencia
## 2.287021 4.119464 3.125000
私のモデルは2.29のペイアウト(つまり、セビージャの勝利の可能性が高い)を予測しているのに対し、betsson.comは3.2というはるかに高いペイアウトを示しているので、私は明らかにセビージャに賭けるべきだ。言い換えると、自作のモデルはホームのセビージャが勝つ予測しているので、セビージャのペイアウトが小さいが、betsson.comはセビージャが負けると予想しているため、ペイアウトが高い。
最終的なゴールの結果にベットすることも可能なので、異なるゴールの結果に対して私のモデルが予測するペイアウトを計算する。
goals_payout <- plyr::laply(0:6, function(home_goal) {
plyr::laply(0:6, function(away_goal) {
1 / mean(home_goals == home_goal & away_goals == away_goal)
})
})
colnames(goals_payout) <- paste("Valencia", 0:6, sep=" - ")
rownames(goals_payout) <- paste("Sevilla", 0:6, sep=" - ")
goals_payout <- round(goals_payout, 1)
print(xtable::xtable(goals_payout, align="cccccccc"), type="html")
Valencia - 0 | Valencia - 1 | Valencia - 2 | Valencia - 3 | Valencia - 4 | Valencia - 5 | Valencia - 6 | |
---|---|---|---|---|---|---|---|
Sevilla - 0 | 16.90 | 12.70 | 20.90 | 42.60 | 142.90 | 666.70 | 4000.00 |
Sevilla - 1 | 11.00 | 8.70 | 14.20 | 33.30 | 105.30 | 250.00 | 1000.00 |
Sevilla - 2 | 14.00 | 11.80 | 18.20 | 38.80 | 102.60 | 666.70 | 1333.30 |
Sevilla - 3 | 25.00 | 22.50 | 35.70 | 93.00 | 166.70 | 2000.00 | 2000.00 |
Sevilla - 4 | 51.90 | 51.30 | 100.00 | 181.80 | 444.40 | 2000.00 | Inf |
Sevilla - 5 | 222.20 | 129.00 | 363.60 | 1000.00 | 4000.00 | 4000.00 | 4000.00 |
Sevilla - 6 | 666.70 | 500.00 | 666.70 | Inf | Inf | Inf | Inf |
最も可能性の高い結果は1 - 1で、予想配当は8.7。betsson.comもこれに同意し、このベットの最低配当である5.3を提供している。十分ではない!bettson.comのペイアウトを見ると、Sevilla - Valencia: 2 - 0のペイアウトは16.0と、私の予想ペイアウトである14よりはるかに良い。これは私の予想配当である 14 よりもはるかに良い配当である!
私のモデルには多くの利点があると思う。概念的に非常にシンプルで、理解しやすく、実装しやすく、拡張しやすい。データのパターンと分布をよく捉えている。これによって、リーガ・エスパニョーラのどのシーズンのどのチームとの試合でも、どのような結果になるかの確率を計算することができる。