Exercice 3
On observe les réalisations \(\boldsymbol{x} = \{x_1,\ldots,x_n\}\) d’une suite de variables aléatoires \(iid\) \(\{X_i\}_{i=1, \dots, n}\) réelles et supérieures à 1, dont la loi \(P_\theta\) est supposée connue à un paramètre \(\theta > 0\) près. Cette loi \(P_\theta\) est une loi continue, appelée loi de Pareto de paramètres \((\theta+1,1)\) dont la densité est définie, pour \(x>1\), par : \[\begin{align*} f_\theta(x) &=\frac{\theta+1}{x^{\theta+2}} \end{align*}\]
- La loi a priori utilisée pour \(\theta\) est une loi exponentielle de paramètre 1, dont la fonction de densité s’écrit : \(g(\theta) = e^{-\theta}\). Écrire alors le modèle bayésien associé.
Question d’intérêt
Elle porte ici sur l’estimation de \(\theta\), le paramètre de la loi de Pareto supposée pour les observations.Modèle d’échantillonnage \[X_i | \theta \overset{iid}{\sim} \text{Pareto}(\theta +1, 1)\] et donc \(\displaystyle f(x_i|\theta) =\frac{\theta+1}{x_i^{\theta+2}}\)
loi(s) a priori \[\theta \sim \mathcal{E}(1)\]
- Montrer que la densité de la loi a posteriori de \(\theta | \boldsymbol{x}\) — notée \(p(\theta|\boldsymbol{x})\) — est proportionnelle à :
\[e^{-\theta}\left(\theta+1\right)^{n}\prod_{i=1}^{n}x_i^{-\theta} \quad ; \quad \theta > 0\]
D’après le théorème de Bayes, on peut écrire que le numérateur de la loi a posteriori est proportionnel au produit de la densité de la loi a priori avec et la vraisemblance des données :
\[\begin{align*} p(\theta | \boldsymbol{x}) &\quad \propto \quad \pi(\theta)f(x_1, ..., x_n | \theta) \\ &\quad \propto \quad e^{-\theta}\prod_{i=1}^n f(x_i| \theta) \\ &\quad \propto \quad e^{-\theta} (\theta+1)^n \prod_{i=1}^n \frac{1}{x_i^{\theta+2}}\\ &\quad \propto \quad e^{-\theta} (\theta+1)^n \left(\prod_{i=1}^n x_i^{-\theta}\right) \left(\prod_{i=1}^n x_i^{-2}\right) \end{align*}\]
Or \((\prod_{i=1}^n x_i^{-2})\) est indépendant de \(\theta\), donc : \[\begin{align*} p(\theta | \boldsymbol{x})\quad \propto \quad e^{-\theta} (\theta+1)^n \prod_{i=1}^n x_i^{-\theta} \qquad \theta \geq 0 \end{align*}\]
- Proposer un algorithme de Metropolis-Hastings pour estimer la loi a posteriori de \(\theta | X_1,...,X_n\). On prendra comme loi instrumentale la loi a priori de \(\theta\) (il s’agira donc d’un algorithme de Metropolis-Hastings indépendant). Expliciter l’estimateur Bayésien de \(\theta\) construit pour le coût quadratique.
NB: Ne pas oublier de faire apparaître les calculs et la formule de la probabilité d’acceptation.
- Initialiser \(\theta^{(0)}\)
- Pour \(t=1, \dots, n+N\) faire :
- Générer aléatoirement \(y^{(t)}\) à partir de la loi \(\mathcal{E}(1)\)
- Calcul de la probabilité d’acceptation \(\alpha^{(t)}\) (voir après)
- Étape d’acceptation-rejet de \(y^{(t)}\)
- Générer aléatoirement \(u^{(t)}\) à partir de la loi \(U[0;1]\)
- \(\theta^{(t)} := \left\{ \begin{array}{l} y^{(t)} \text{ si } u^{(t)} \leq \alpha^{(t)} \\ \theta^{(t-1)} \text{ sinon} \end{array} \right.\)
On détaille maintenant le calcul de la probabilité d’acceptation : \[\alpha^{(t)} = min \left( 1, \frac{p(y^{(t)} | \boldsymbol{x})}{p(\theta^{(t-1)} | \boldsymbol{x})}\frac{\pi(\theta^{(t-1)})}{\pi(y^{(t)})} \right)\]
\[\begin{align*} \require{cancel} \dfrac{p(y^{(t)} | \boldsymbol{x})}{p(\theta^{(t-1)} | \boldsymbol{x})}\dfrac{\pi(\theta^{(t-1)})}{\pi(y^{(t)})} & = \dfrac{p(y^{(t)} | \boldsymbol{x})}{p(\theta^{(t-1)} | \boldsymbol{x})}\dfrac{e^{-\theta^{(t-1)}}}{e^{-y^{(t)}}} \\ & = \dfrac{(y^{(t)} +1)^n(\prod_{i=1}^n x_i^{-y^{(t)}})\cancel{e^{-y^{(t)}}}}{(\theta^{(t-1)}+1)^n\left(\prod_{k=1}^n x_i^{-\theta^{(t-1)} }\right)\cancel{e^{-\theta^{(t-1)}}}}\dfrac{\cancel{e^{-\theta^{(t-1)}}}}{\cancel{e^{-y^{(t)}}}} \end{align*}\]
Finalement \(\alpha = min \left( 1, \dfrac{(y +1)^n \prod_{i=1}^n x_i^{-y}}{(\theta^{(t-1)} +1)^n \prod_{k=1}^n x_i^{-\theta^{(t-1)} }} \right)\)
Cet algorithme échantillonne alors selon la loi a posteriori du paramètre \(\theta\) :
Après une phase de chauffe de longueure \(n\), nous recueillons un \(N\)-échantillon de \(\theta\) généré par cet algorithme de Metropolis-Hasting \(\left(\theta^{(n+1)}, ... , \theta^{(n+N)}\right)\). Appliquer la méthode de Monte-Carlo permet d’obtenir un estimateur de la moyenne a posteriori de \(\theta\).
L’estimateur bayésien pour une fonction de coût quadratique est la moyenne a posteriori, calculée par : \[\widehat{E}(\theta|\boldsymbol{x}) = \frac{1}{N} \sum_{k=1}^{N}\theta^{(n + k)}\]
- Quel résultat théorique garantit sa convergence ? Expliquer brièvement.
La convergence de l’algorithme de Metropolis-Hastings est garantie par le théorème ergodique. L’algorithme de Metropolis-Hastings échantillonne les réalisations de \(\theta\) à l’aide d’une chaine de Markov dont la distribution stationnaire est la loi a posteriori de \(\theta\). Le théorème ergodique permet donc d’appliquer la loi des grands nombres aux réalisations de cette chaîne.