UPDATE: 2024-01-23 16:05:44.967038

はじめに

このノートは「ベイズ統計」に関する何らかの内容をまとめ、ベイズ統計への理解を深めていくために作成している。今回はRasmus Bååth先生のブログに記載されているサッカーチームの強さを推定する内容を参考にさせていただき、写経しながら、ところどころ自分用の補足をメモすることで、自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。

こちらのブログではJAGSを利用しており、再現はStanで行うため、結果は完全に一致しない。また、JAGSとStanでは正規分布のパラメタの与え方が異なるので、そのあたりは修正を行っているが、ブログの数式はそのまま表示しているので、数式とStanモデルでは差異がある点は注意が必要である。

Modeling Match Results in La Liga Using a Hierarchical Bayesian Poisson Model.

この分析は、

これが本分析の目的とのこと。まずは必要なデータやパッケージを準備しておく。ブログで使用されている通り、使用するデータはラ・リーガの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

特徴に関して、下記の通り説明を記載しておく。

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

Modeling Match Results: Iteration 1

サッカーの試合はポアソン分布に従うと一般的に知られている。サッカーの試合はアディショナルタイムがあれど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)

ここで現在のモデルに問題があることがわかる。シミュレーションデータは、ホームチームとアウェイチームが入れ替えても結果は変わらない。過去のデータを見ると、セビージャがバレンシアに勝つのはホームチームである場合が多い。ホームチームであることの利点を考慮しているモデルが必要である。

Modeling Match Results: Iteration 2

ホームアドバンテージを考慮するためのモデルの変更は、ベースラインをホームベースラインとアウェイベースラインの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_baselineaway_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)

セビージャもバレンシアもホームチームとしてプレーした方が勝ちやすいことがわかる。

Modeling Match Results: Iteration 3

データセット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に優勝するまでに徐々に力をつけていたことがわかる。

  • 2008-09: バルセロナ / レアル・マドリード
  • 2009-10: バルセロナ / レアル・マドリード
  • 2010-11: バルセロナ / レアル・マドリード
  • 2011-12: レアル・マドリード / バルセロナ
  • 2012-13: バルセロナ / レアル・マドリード
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 →"
  )

Predicting the End Game of La Liga 2012/2013

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_predmatch_resulthome_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

Calculating the Predicted Payout for Sevilla vs Valencia, 2013-06-01

ベイズ・モデリングと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 よりもはるかに良い配当である!

Wrap-up

私のモデルには多くの利点があると思う。概念的に非常にシンプルで、理解しやすく、実装しやすく、拡張しやすい。データのパターンと分布をよく捉えている。これによって、リーガ・エスパニョーラのどのシーズンのどのチームとの試合でも、どのような結果になるかの確率を計算することができる。