Introduction à la statistique bayésienne pour les ingénieurs et les médecins

• La seule différence avec l'algorithme de Metropolis est que la distribution proposition n'est plus symétrique, ceci signifie que d'aller de \(\boldsymbol{\theta}_1\) vers \(\boldsymbol{\theta}_2\), n'est pas aussi facile (ou difficile) que d'aller de \(\boldsymbol{\theta}_2\) vers \(\boldsymbol{\theta}_1\)

• On calcule \(\alpha = \dfrac{p(\tilde{\boldsymbol{\theta}}|y) q(\boldsymbol{\theta}^{(k)} |\tilde{\boldsymbol{\theta}})}{p(\boldsymbol{\theta}^{(k)}|y) q(\tilde{\boldsymbol{\theta}}|\boldsymbol{\theta}^{(k)}) }\) et l'algorithme se déroule comme pour l'algorithme de Metropolis

• Utilise comme distribution de proposition une distribution qui satisfait à la propriété \(q(\tilde{\boldsymbol{\theta}}|\boldsymbol{\theta}) = q(\tilde{\boldsymbol{\theta}} - \boldsymbol{\theta})\)

• Souvent la distribution de proposition est une loi multinormale, avec un taux d'acceptation attendu de l'ordre de 25%

  • Moyenne \(\boldsymbol{\theta}^k\) et de matrice de variance-covariance \(\Sigma\), dont le choix est laissé à l'utilisateur (ou au logiciel)

  • Ce choix est important car la matrice intervient dans le temps de convergence : plus les valeurs dans la diagonale de \(\Sigma\) seront élevées et plus la convergence sera rapide mais le taux d'acceptation bas (donc les estimations médiocres)

• On calibre la matrice \(\Sigma\) pour obtenir un taux d'acceptation de l'ordre de 25%

  • Étude pilote avec un matrice diagonale pour obtenir ce taux

  • On modifie ensuite la matrice en introduisant des valeurs hors de la diagonale (corrélations supposées entre les paramètres à estimer)

  • Finalement on retient 5,8/d (nombre de paramètres) fois cette matrice (Gelman et al 1995)

Plutôt que modéliser directement un pourcentage, on peut préférer modéliser le logit, transformant ainsi une quantité entre 0 et 1 en une quantité dont le domaine est \(\mathbb{R}\). On peut transformer la formule de la vraisemblance binomiale pour faire apparaître \(\theta = \log\left(\dfrac{\pi}{1-\pi}\right)\) plutôt que \(\pi\) : \(p(y|\pi) = \left(\begin{array}{c} N \\ y \end{array}\right) \pi^{y} (1-\pi){N-y}\) devient \(p(y|\theta) = \left(\begin{array}{c} N \\ y \end{array}\right) \left( \dfrac{e^{\theta}}{1+ e^{\theta}} \right)^{y} \left( \dfrac{1}{1+ e^{\theta}} \right)^{N-y} = \left(\begin{array}{c} N \\ y \end{array}\right) \dfrac{ e^{y \theta}}{(1+ e^{\theta})^N}\)

• Vraisemblance: \(p(y|\theta) = \left(\begin{array}{c} N \\ y \end{array}\right) \dfrac{ e^{y \theta}}{(1+ e^{\theta})^N}\)

• Prior sur \(\theta\): \(\theta \sim \mathcal{N}(\mu_0, \sigma^2_0)\), où \(\mu_0\) et \(\sigma^2_0\) sont des quantités fixes: 0 et 10.000

• A posteriori: \(p(\theta|y) \propto \dfrac{e^{y \theta}}{(1 + e^{\theta})^N} \exp{\left[-\dfrac{1}{2}\left( \dfrac{\theta - \mu_0}{\sigma_0} \right)^2\right]}\)

• Pour estimer \(\theta\), on va utiliser un algorithme « Metropolis à marches aléatoires » ("random-walk Metropolis")

• On ajoute un petit truc: pour éviter les dépassements de capacités de calcul (à cause des exponentielles), on va travailler sur l'échelle logarithmique

lorsque la chaîne est en \(\theta^{(k)}\), l'algorithme est ainsi :

1. échantillonner \(\tilde{\theta}\) de la distribution \(\mathcal{N}(\theta^{(k)}, s^2)\)

2. calculer le \(r\), minimum entre 0 et \(\alpha = \log{\dfrac{p(\tilde{\theta}|y)}{ p(\theta^{(k)}|y)}}\)

3. \(\alpha = (\tilde{\theta} - \theta^{(k)}) \left( y - \dfrac{\tilde{\theta} + \theta^{(k)} - 2\mu_0}{2 \sigma^2_0} \right) - N \log{\dfrac{1+e^{\tilde{\theta}}}{1+e^{\theta^{(k)}}}}\)

4. si \(r \geqslant 1\), \(\tilde{\theta}\) est accepté et \(\theta^{(k+1)} = \tilde{\theta}\)

5. si \(r < 1\), \(\tilde{\theta}\) est accepté comme nouvelle valeur avec une probabilité \(r\) : on tire \(u\) d'une uniforme [0,1], si \(u \leqslant r\), on accepte la nouvelle valeur et \(\theta^{(k+1)} = \tilde{\theta}\) et si \(u > r\), on rejette la nouvelle valeur et \(\theta^{(k+1)} = \theta^{(k)}\)

\(\rightarrow\) décider de la valeur de \(s^2\) \(\Leftarrow\) acceptation de l'ordre de 25%