---
title: "Modélisation du nombre de buts"
output: html_notebook
---

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.

```{r}
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}
## 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).
```{r}
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
```{r}
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 ```r acc_rate``` 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*
```{r}
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
```{r}
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
```{r}
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)))

```
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).
