Exercice 3 : Algorithme(s) de Metropolis-Hastings pour l'exemple historique (modèle Beta-Bernoulli)

En reprenant l’exemple historique, programmer un algorithme de Metropolis-Hastings indépendant pour estimer la loi a posteriori du paramètre \(\theta\) (la probabilité d’avoir une fille à la naissance). On prendra comme loi instrumentale la loi a priori de \(\theta\), et on utilisera dans un premier temps un a priori uniforme. On considérera les \(493\,472\) naissances observées à Paris entre 1745 et 1770, dont \(241\,945\) furent des filles.

1. Programmer une fonction calculant le numérateur du posterior, dont on rappelle l’expression : \(p(\theta|n,S)\propto\theta^S(1-\theta)^{n-S}\) où \(S = 241\,945\) et \(n = 493\,472\) (prévoir un argument permettant de renvoyer - ou non - le logarithme de ce numérateur).

   post_num_hist <- function(theta, log = FALSE) {
   
       # data
       n <- 493472
       S <- 241945
   
       if (log) {
           num <- ...
       } else {
           num <- ...
       }
       return(num)
   }

   post_num_hist(0.2, log = FALSE)
   post_num_hist(0.6, log = FALSE)
   post_num_hist(0.2, log = TRUE)
   post_num_hist(0.6, log = TRUE)

2. Programmer l’algorithme de Metropolis-Hastings correspondant qui renvoie un vecteur de taille \(N\) échantillonné selon la loi a posteriori (on utilisera l’a priori — la loi Uniforme — comme loi instrumentale). La fonction doit également renvoyer le vecteur des probabilités d’acceptation. Que se passe-t-il si l’on ne calcule pas cette probabilité d’acceptation dans l’échelle logarithmique ?

   monMH <- function(niter, post_num) {
   
       x_save <- numeric(length = niter)  # sauver les valeurs échantillonées - initialisation d'un vecteur de 0s de taille niter
       alpha <- numeric(length = niter)  # sauver les probabilités d'acceptation - initialisation d'un vecteur de 0s de taille niter
   
       # initialisation de x0
       x <- ...
   
       for (t in 1:niter) {
   
           # proposer une valeur a partir de la loi instumentale (a priori
           # uniforme)
           y <- ...
   
           # calculer la probabilité d'acceptation-rejet
           alpha[t] <- ...
           # message('alpha :', alpha) # permet d'afficher alpha
   
           # etape d'acceptation-rejet
           u <- ...
           if (...) {
               x_save[t] <- ...
           } else {
               x_save[t] <- ...
           }
   
           # mise-à-jour de la valeur courante x
           x <- x_save[t]
       }
       return(list(theta = x_save, alpha = alpha))
   }

3. Comparer la densité a posteriori obtenue avec l’algorithme M-H sur 2000 itérations à la densité théorique. On prendra soin d’enlever les 500 premières itérations afin d’obtenir la convergence de la chaîne de Markov. Commenter ces résultats, notamment au vu des différentes probabilités d’acceptation calculées au cours de l’algorithme et des différentes valeurs échantillonnées pour \(\theta\).

4. Afin de mieux illustrer le fonctionnement de l’algorithme de Metropolis-Hastings, nous allons modifier la loi stationnaire cible de la chaîne de Markov en ne considérant que \(100\) observations, parmi lesquelles \(49\) naissances sont des filles, et en changeant l’a priori par une loi \(\text{Beta}(\alpha = 3, \beta=3)\).
   Reprogrammer une fonction post_num_beta() calculant le numérateur de la loi a posteriori correspondante (d’après la formule vue dans le cours), puis faire tourner l’algorithme monMH() avec \(10\,000\) itérations.
   NB : On conserve donc une loi instrumentale uniforme (celle-ci n’est pas obligatoirement égale à l’a priori) — on reste donc dans le cas d’un algorithme de Metropolis-Hastings indépendant. Si le temps le permet, utiliser également l’a priori \(\text{Beta}\) pour la loi instrumentale (ce qui nécessite de modifier le calcul de la probabilité \(\alpha_t\).

   post_num_beta <- function(theta, log = TRUE, a = 3, b = 3) {
       ...
   }
   
   sampleMH_beta <- monMH(10000, post_num = post_num_beta)
   
   # evaluation graphique
   ...

5. Reprenons maintenant la loi stationnaire initiale utilisant les données de l’exemple historique avec un a priori uniforme.
   Programmer un algorithme de Metropolis-Hastings à marche aléatoire (on pourra commencer par prendre un pas aléatoire uniforme d’amplitude \(0,02\), puis (si le temps le permet) un pas gaussien d’écart-type \(0,02\)). Étudier à nouveau les résultats ainsi obtenus (on pourra faire varier la taille du pas aléatoire).

   monMH_marchalea <- function(niter, post_num, a = 3, b = 3) {
       ...
       y <- x + runif(n = 1, min = -0.01, max = 0.01)
       ...
   }
   sampleMH_marche <- monMH_marchalea(20000, post_num = post_num_beta_hist)
   ...