Modélisation du nombre de buts

Ici nous considérons le modèle Bayésien suivant \[ N_i \mid \lambda \stackrel{iid}{\sim} Poisson(\lambda), \qquad i = 1, \ldots, n, \] avec la loi a priori \(\lambda \sim Gamma(\alpha, \beta)\). Attention \(\beta\) est le paramètre d’intensité (rate).

Un peu de calcul et nous trouvons la loi a posteriori qui est une \(Gamma(\tilde \alpha, \tilde \beta)\) avec \[ \tilde \alpha = \alpha + \sum_{i=1}^n N_i, \qquad \tilde \beta = n + \beta, \]

Afin de vérifier que nous savons coder un Metropolis–Hastings sur un vrai modèle statistique on va utiliser la loi cible \[ \pi(\lambda \mid N_1, \ldots, N_n) \propto \lambda^{\alpha + \sum_{i=1}^n N_i} - 1 - (n + \beta) \log \lambda. \] Evidemment cela ne sert à rien normalement mais cela a un intérêt pédagogique.

football <- function(nchain, data, hyper_param,
                     prop_sd, init = 1){
  chain <- rep(NA, nchain + 1)
  chain[1] <- init
  nobs <- length(data)
  
  acc_rate <- 0
  
  for (i in 1:nchain){
    ## Générons un candidat (marche alétoire à l'échelle log)
    prop <- rlnorm(1, log(chain[i]), prop_sd)
    
    ## Calcul de la probabilité d'acceptation
    top <- -(nobs + hyper_param[2]) * prop +
      (hyper_param[1] + sum(data) - 1) * log(prop)
    bottom <- -(nobs + hyper_param[2]) * chain[i] +
      (hyper_param[1] + sum(data) - 1) * log(chain[i])
    
    ## Attention ici le ratio des noyaux de proposition n'est pas 1
    acc_prob <- exp(top - bottom) * prop / chain[i]
    
    ## Mise a jour
    if (runif(1) < acc_prob){
      chain[i+1] <- prop
      acc_rate <- acc_rate + 1
    } else
      chain[i+1] <- chain[i]
  }
  
  acc_rate <- acc_rate / nchain
  return(list(chain = chain, acc_rate = acc_rate))
}

Essai sur les données.

## Récupération des données depuis le web
url <- "https://www.football-data.co.uk/mmz4281/2324/F1.csv"
data <- read.csv(url)

## Data pre-processing
team <- "Lyon"
idx_home <- which(data$HomeTeam == team)
idx_away <- which(data$AwayTeam == team)

lyon_goals <- c(data[idx_home,"FTHG"],
                data[idx_away,"FTAG"])

Choisissons des hyper-paramètres raisonnables (prior elicitation).

alpha <- 1
beta <- 0.1
prior <- function(x)
  dgamma(x, alpha, beta)

plot(prior, from = 0, to = 10,
     xlab = expression(lambda),
     ylab = "Density", ylim = c(0, 1))

Ici on a choisi des hyper-paramètres afin que la loi a priori soit pas trop informative.

Testons notre function

nchain <- 50000
hyper_param <- c(alpha, beta)
prop_sd <- 0.35
ans <- football(nchain, lyon_goals, hyper_param,
                prop_sd)

chain <- ans$chain
acc_rate <- ans$acc_rate

Le taux d’acceptation de notre chaîne est 0.55166 ce qui est une bonne valeur pour une chaine en dimension 1. Confirmons cela en visualisant notre chaine. Pour cela on peut utiliser le package coda

library(coda)
plot(mcmc(chain))

acf(chain)

Ici pas besoin d’enlever la période de chauffe car l’on semble être déjà dans le régime stationnaire. De plus l’ACF montre une dépendance assez limitée. On est donc plus que content.

Regardons avec impatience si ça colle bien la théorie

nobs <- length(lyon_goals)
true_posterior <- function(x)
  dgamma(x, alpha + sum(lyon_goals),
         nobs + beta)

col <- RColorBrewer::brewer.pal(3, "Set2")
par(mar = c(4, 5, 0.5, 0.5))
plot(density(chain), col = col[1], main = "",
     xlab = expression(lambda), lwd = 2)
curve(true_posterior, from = 0, to = 3,
      add = TRUE, col = col[2], lwd = 2 )

legend("topright", c("Estimated", "True"),
       lty = 1, bty = "n", col = col)

Conclusion on est trop fort !

Passons maintenant à la prédiction. Pour cela nous allons prendre comme estimateur la moyenne de la loi prédictive. On se rappelle que la loi prédictive a posteriori est donnée par \[ \pi(x_* \mid N_1, \ldots, N_n) = \int f(x_* \mid N_1, \ldots, N_n, \lambda) \pi(\lambda \mid N_1, \ldots, N_n) d\lambda, \] or ici puisque les matchs sont supposés indépendants, cela se simplifie en \[ \pi(x_* \mid N_1, \ldots, N_n) = \int f(x_* \mid \lambda) \pi(\lambda \mid N_1, \ldots, N_n) d\lambda. \] L’espérance \(\mu = \mathbb{E}(X_*)\) sous cette loi est donc \[ \mu = \sum_{k \geq 0} k \int f(k \mid \lambda) \pi(\lambda \mid N_1, \ldots, N_n) d\lambda = \int \sum_{k \geq 0} k f(k \mid \lambda) \pi(\lambda \mid N_1, \ldots, N_n) d\lambda, \] qui s’écrit comme \[ \mu_N = \mathbb{E}\left\{ \mathbb{E}\left[X_* \mid \lambda \right]\right\}, \] que l’on peut facilement estimer selon la loi des grands nombres. Concrètement, en utilisant encore notre chaine de Markov, on fait

sim <- rpois(length(chain), chain)
sim_factor <- factor(sim, levels = 0:10)
freq <- table(sim_factor) / length(chain)
plot(freq, xlab = "Nombre de buts",
     ylab = "Fonction de masse",
     col = col[1])

## Estimation et IC
c(mean(sim), quantile(sim, prob = c(0.025, 0.975)))

             2.5%    97.5% 
1.082778 0.000000 4.000000 

Ici on devrait plutôt utiliser comme estimateur le mode de la loi prédictive qui ici est 1 (mais 0 est très proche aussi). Donc un pari raisonnable serait de dire que Lyon marque moins de 1.5 buts. Bon il y a aussi le pari super sur “marque moins de 2.5 buts” mais ça ne devrait rien rapporter.