class: center, middle, inverse, title-slide .title[ # Statistiques computationnelles ] .author[ ###
Charlotte Baey
] .date[ ###
M1 MA - 2025/2026
] --- <style> .remark-slide-content { background-color: #FFFFFF; border-top: 80px solid #16A085; font-size: 20px; line-height: 1.5; padding: 1em 2em 1em 2em } .my-one-page-font { font-size: 20px; } .remark-slide-content > h1 { font-size: 38px; margin-top: -85px; } .inverse { background-color: #16A085; border-top: 80px solid #16A085; text-shadow: none; background-position: 50% 75%; background-size: 150px; font-size: 40px } .title-slide { background-color: #16A085; border-top: 80px solid #16A085; background-image: none; } .remark-slide-number { position: absolute; } .remark-slide-number .progress-bar-container { position: absolute; bottom: 0; height: 4px; display: block; left: 0; right: 0; } .remark-slide-number .progress-bar { height: 100%; background-color: grey; } .left-column { width: 20%; height: 92%; float: left; } .left-column h2:last-of-type, .left-column h3:last-child { color: #000; } .right-column { width: 75%; float: right; padding-top: 1em; } .left-column2 { width: 60%; height: 92%; float: left; } .right-column2 { width: 35%; height: 92%; float: left; } </style> <script type="text/x-mathjax-config"> MathJax.Hub.Config({ TeX: { Macros: { yellow: ["{\\color{yellow}{#1}}", 1], orange: ["{\\color{orange}{#1}}", 1], green: ["{\\color{green}{#1}}", 1] }, loader: {load: ['[tex]/color']}, tex: {packages: {'[+]': ['color']}} } }); </script> <style type="text/css"> .left-code { width: 40%; height: 92%; float: left; } .right-plot { width: 59%; float: right; padding-left: 1%; } </style> <style type="text/css"> .left-plot { width: 59%; float: left; } .right-code { width: 40%; float: right; padding-left: 1%; } </style> # Quelques informations pratiques ### Plan du cours 1. Méthodes de ré-échantillonnage 2. Méthodes de Monte-Carlo 3. Introduction aux statistiques bayésiennes 4. Algorithme EM (s'il reste du temps) ### Organisation - 2 séances de cours d'1h30 par semaine (`\(\times\)` 11 semaines) - 2 séances de TD/TP de 2h par semaine (`\(\times\)` 12 semaines) ### Evaluation - 1 DS intermédiaire d'une durée de 2h - 1 Projet **à effectuer en binôme** - 1 DS final d'une durée de 3h .red[**Aucun document autorisé lors des examens.**] --- # Sommaire .pull-left[ **1. Méthodes de ré-échantillonnage** - [Cours 1](#c1) (12/01/2026) - [Cours 2](#c2) (12/01/2026) - [Cours 3](#c3) (19/01/2026) - [Cours 4](#c4) (20/01/2026) - [Cours 5](#c5) (26/01/2026) **2. Méthodes de Monte-Carlo** - [Cours 6](#c6) (27/01/2026) - [Cours 7](#c7) (03/02/2026) - [Cours 8](#c8) (06/02/2026) - [Cours 9](#c9) (09/02/2026) - [Cours 10](#c10) (10/02/2026) - [Cours 11](#c11) (16/02/2026) - [Cours 12](#c12) (02/03/2026) - [Cours 13](#c13) (03/03/2026) - [Cours 14](#c14) (09/03/2026) ] --- name: c1 class: inverse, middle, center # Introduction --- class: my-one-page-font # C'est quoi les statistiques computationnelles ? <br> <br> - C'est le recours (plus ou moins) intensif à l'ordinateur pour répondre à des questions statistiques que l'on ne sait pas (ou difficilement) résoudre autrement. - On utilise/développe/étudie des algorithmes, des astuces numériques/statistiques/computationnelles - L'objectif est de faire de l'inférence, d'étudier la robustesse de méthodes statistiques, de traiter de grands jeux de données, ... --- class: inverse, middle, center # I. Méthodes de ré-échantillonnage --- # Notion d'échantillon Qu'est-ce qu'un échantillon ? - une suite de variables aléatoires `\(\mathcal{X} = (X_1, \dots, X_n)\)` - dont on observe une réalisation `\(\mathcal{X}(\omega) = (X_1(\omega), \dots, X_n(\omega))\)` --- # Notion d'échantillon Qu'est-ce qu'un échantillon ? - une suite de variables aléatoires `\(\mathcal{X} = (X_1, \dots, X_n)\)` - dont on observe une réalisation `\(\mathcal{X}(\omega) = (X_1(\omega), \dots, X_n(\omega))\)` Que représente `\(\omega\)` ? - l'aléa autour de l'expérience (ex. : `\(n\)` lancers d'une pièce de monnaie) - cet aléa `\(\omega \in \Omega\)` est transporté dans `\(\mathbb{R}\)` via `\(X_i\)` <!-- - on a seulement accès à la variabilité sur `\(\mathbb{R}\)` --> - en général, on ne dispose que d'une seule réalisation, pour un `\(\omega\)` donné --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-4-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-5-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-6-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-7-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-8-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-9-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-10-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-11-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-12-1.png" width="2000px" /> --- # Statistique paramétrique <img src="slides_25_files/figure-html/unnamed-chunk-13-1.png" width="2000px" /> --- # Un exemple simple Soit `\((X_1,\dots,X_n)\)` un échantillon gaussien i.i.d. de loi `\(\mathcal{N}(\theta,1)\)`, et `\(\hat{\theta}\)` l'EMV de `\(\theta\)`. Quelle est la loi de `\(\hat{\theta}\)` ? -- .pull-left[ ``` r theta_vrai <- 2; n <- 100; N <- 500 hat_theta <- rep(0,N) for (i in 1:N){ ech_i <- rnorm(n,theta_vrai,1) hat_theta[i] <- mean(ech_i) } hist(hat_theta,freq=F) curve(dnorm(x,theta_vrai,1/sqrt(n))) ``` ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-15-1.png" width="432" /> ] --- class: my-one-page-font # Jackknife Comment construire de nouveaux échantillons ? <img src="slides_25_files/figure-html/unnamed-chunk-16-1.png" width="720" /> --- class: my-one-page-font # Jackknife Comment construire de nouveaux échantillons ? <img src="slides_25_files/figure-html/unnamed-chunk-17-1.png" width="720" /> --- # Jackknife #### 1. Réduction du biais - Estimation du biais : `$$\hat{b}_{jack} = (n-1)(\hat{\theta}_{jack} - \hat{\theta}),$$` avec `\(\hat{\theta}_{jack} = \frac{1}{n}\sum_{i=1}^n \hat{\theta}_{(i)}\)`. - Pseudo-valeurs : `$$\tilde{\theta}_{(i)} = n\hat{\theta} - (n-1)\hat{\theta}_{(i)}$$` <!-- Sur l'exemple précédent, cela donne : --> <!-- \tilde{\theta}_{(1)} = \tilde{\theta}_{(2)} = \tilde{\theta}_{(3)} = \tilde{\theta}_{(7)} = \tilde{\theta}_{(8)} = \tilde{\theta}_{(10)} = \tilde{\theta}_{(12)} = 13\times0.5384 - 12\times0.50 = 0.9992 \tilde{\theta}_{(4)} = \tilde{\theta}_{(5)} = \tilde{\theta}_{(6)} = \tilde{\theta}_{(9)} = \tilde{\theta}_{(11)} = \tilde{\theta}_{(13)} = 13\times0.5384 - 12\times0.5833 = -0.0004 --> - Estimateur jackknife corrigé du biais : `$$\hat{\theta}^*_{jack} = \hat{\theta} - \hat{b}_{jack} = \frac{1}{n}\sum_{i=1}^n \tilde{\theta}_{(i)}$$` .red[**Réduire le biais n'implique pas nécessairement une amélioration de l'estimateur (au sens du risque quadratique)**] --- # Jackknife #### 2. Estimation de la variance `$$\hat{s}^2_{jack} = \frac{n-1}{n} \sum_{i=1}^n \big(\hat{\theta}_{(i)} - \hat{\theta}^*_{jack}\big)^2$$` avec les pseudo-valeurs : `\(\hat{s}^2_{jack} = \frac{1}{n(n-1)} \sum_{i=1}^n \big(\tilde{\theta}_{(i)} - \tilde{\theta}\big)^2\)` -- #### 3. Construction d'intervalles de confiance Si existence d'un TCL : - en utilisant l'estimateur jackknife de la variance : `$$\hat{I}_{jack} = \left[\hat{\theta} - q^{\mathcal{N(0,1)}}_{1-\alpha/2} \hat{s}_{jack} ; \hat{\theta} + q^{\mathcal{N(0,1)}}_{1-\alpha/2} \hat{s}_{jack} \right]$$` - en utilisant les deux estimateurs : `$$\hat{I}_{jack} = \left[\hat{\theta}_{jack} - q^{\mathcal{N(0,1)}}_{1-\alpha/2} \hat{s}_{jack} ; \hat{\theta}_{jack} + q^{\mathcal{N(0,1)}}_{1-\alpha/2} \hat{s}_{jack} \right]$$` --- name: c2 # Résumé et remarques <br> - le jackknife est une méthode **non-paramétrique** permettant d'estimer le biais et la variance d'un estimateur à l'aide de simulations - la consistance est garantie pour un grand nombre d'estimateurs **suffisamment réguliers** - les hypothèses de régularité sont toutefois plus strictes que pour un TCL par exemple : - la delta-méthode requiert la différentiabilité de `\(g\)` en `\(\mu=\mathbb{E}(X_1)\)` - la consistance de `\(\hat{s}_{jack}\)` requiert que `\(g'\)` soit continue en `\(\mu\)` - ex. d'estimateur pour lequel `\(\hat{s}_{jack}\)` n'est pas consistant : la médiane empirique (malgré existence TCL) - extensions possibles reposant sur des hypothèses moins fortes : - le *delete-d* jackknife - le jackknife infinitésimal <!-- - peut-être utilisé pour la médiane empirique --> <!-- - peut permettre d'estimer la distribution de `\(\hat{\theta}\)` --> <span style="color:#16A085">**fin du cours 1**</span>- --- # Du jackknife au bootstrap - `\(\hat{\theta} = T(\underbrace{X_1,\dots,X_n}_{\text{observations}},\underbrace{\omega_1,\dots,\omega_n}_{\text{poids des obs.}}) = T(X_1,\dots,X_n,\frac{1}{n},\dots,\frac{1}{n})\)` - réplication jackknife : `\(\hat{\theta}_{(i)} = T(X_1,\dots,X_{i-1},X_i,X_{i+1},X_n,\frac{1}{n-1},\dots,\frac{1}{n-1},0,\frac{1}{n-1},\dots,\frac{1}{n-1})\)` -- - Idée du bootstrap : mettre des poids **aléatoires** -- - Procédure de ré-échantillonnage : tirer uniformément et **avec remise** parmi les observations de `\(\mathcal{X}\)`, pour construire un échantillon *bootstrap* de taille `\(n\)` noté `\(\mathcal{X}^*\)` : - Puis sur chaque échantillon bootstrap, on construit la statistique bootstrapée `\(\hat{\theta}^* = T(\mathcal{X}^*)\)` --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-18-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-19-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-20-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-21-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-22-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-23-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-24-1.png" width="2000px" /> --- # Bootstrap <img src="slides_25_files/figure-html/unnamed-chunk-25-1.png" width="2000px" /> --- # Bootstrap .pull-left[ ### Monde réel .left[ <span style="color:white">**on raisonne conditionnellement à `\(F_n\)`**</span> - échantillon `\(\mathcal{X} = (X_1,\dots,X_n)\)` - `\(X_i\)` de loi inconnue `\(F\)` - paramètre `\(\theta(F)\)` - estimateur `\(\hat{\theta} = T(\mathcal{X})\)` - loi de `\(\hat{\theta}\)` : `\(G\)` inconnue ] ] .pull-right[ ### Monde Bootstrap .left[ <span style="color:red">**on raisonne conditionnellement à `\(F_n\)`**</span> - échantillon bootstrap `\(\mathcal{X}^* = (X_{1}^*,\dots,X_{n}^*)\)` - `\(X_{i}^*\)` de loi connue `\(F_n\)` - paramètre `\(\theta(F_n) = \hat{\theta}\)` - statistique bootstrapée `\(\hat{\theta}^* = T(\mathcal{X}^*)\)` - loi de `\(\hat{\theta}^*\)` : `\(G^*\)` connue ] ] </br> .center[ `\(G\)` inconnue `\(\longrightarrow\)` `\(G^*\)` connue `\(\longrightarrow\)` `\(\hat{G}^*_B\)` approximation bootstrap ] </br> <span style="color:#16A085">**fin du cours 2**</span> --- name: c3 # Fonction de répartition empirique ``` ## [1] -0.88 -0.75 1.38 0.24 0.11 1.20 -0.46 0.64 0.42 0.78 ``` -- <img src="slides_25_files/figure-html/unnamed-chunk-27-1.png" width="720" /> --- # Fonction de répartition empirique ``` ## [1] -0.88 -0.75 1.38 0.24 0.11 1.20 -0.46 0.64 0.42 0.78 ``` <img src="slides_25_files/figure-html/unnamed-chunk-29-1.png" width="720" /> --- # Fonction de répartition empirique - Echantillon exponentiel <img src="slides_25_files/figure-html/unnamed-chunk-30-1.png" width="1152" /> - Echantillon uniforme <img src="slides_25_files/figure-html/unnamed-chunk-31-1.png" width="1152" /> <span style="color:#16A085">**fin du cours 3**</span> --- name: c4 # Bootstrap : principe du plug-in - on remplace `\(F\)` par sa version empirique `\(F_n\)` <img src="slides_25_files/figure-html/unnamed-chunk-32-1.png" width="864" /> <!-- .pull-left[ --> <!-- .center[ --> <!-- espérance : `\(\theta = \int x dF(x)\)` ]]--> <!-- .pull-right[ --> <!-- espérance : `\(\int x dF_n(x) = \frac{1}{n} \sum_{i=1}^n X_i\)`] --> --- # Bootstrap : principe du plug-in #### Exemple `\(\theta = \mathbb{E}(X_1) = \int x dF(x)\)` -- Méthode du plug-in : on estime `\(\theta\)` par `\(\hat{\theta} = \int x dF_n(x) = \frac{1}{n}\sum_{i=1}^n X_i\)`. -- Dans l'échantillon bootstrap, les `\(X_{b,i}^*\)` sont distribuées selon la loi `\(F_n\)`, et ont pour espérance `\(\theta^* = \int x dF_n(x) = \hat{\theta}\)` -- La statistique bootstrapée `\(\hat{\theta}^*_b = \frac{1}{n}\sum_{i=1}^n X_{b,i}^*\)` est donc à `\(\theta^*\)` ce que `\(\hat{\theta}\)` est à `\(\theta\)`. --- # Bootstrap #### Exemple 2 : estimation du biais de `\(\hat{\theta}\)` Paramètre `\(\theta(F)\)` estimé par `\(\hat{\theta} = \theta(F_n)\)` -- Biais : `\(b(\hat{\theta}) = \mathbb{E}_F(\hat{\theta}) - \theta = \int xdG(x) - \theta(F)\)` -- Dans le monde bootstrap : `$$b^*(\hat{\theta}) = \int x dG^*(x) - \theta(F_n) = \int x dG^*(x) - \hat{\theta}$$` -- Estimateur du bootstrap : `$$\hat{b}^*(\hat{\theta}) = \int x d\hat{G}^*_B(x) - \hat{\theta} = \frac{1}{B} \sum_{i=1}^B \hat{\theta}^*_b - \hat{\theta}$$` --- # Bootstrap #### Exemple 3 : estimation de la variance de `\(\hat{\theta}\)` Variance : `\(\text{Var} (\hat{\theta}) = \mathbb{E}_F\big[(\hat{\theta} - \mathbb{E}(\hat{\theta}))^2\big] = \int (x - \int x dG(x))^2dG(x)\)` -- Dans le monde bootstrap : `$$\text{Var}^*(\hat{\theta}) = \int (x - \int x dG^*(x))^2 dG^*(x)$$` -- Estimateur du bootstrap : `$$\hat{\text{Var}}^*(\hat{\theta}) = \int (x - \int x d\hat{G}_B^*(x))^2 d\hat{G}_B^*(x) = \frac{1}{B} \sum_{b=1}^B \big(\hat{\theta}^*_b - \frac{1}{B} \sum_{b=1}^B \hat{\theta}^*_b \big)^2$$` `\(\rightarrow\)` c'est la variance empirique des statistiques bootstrapées `\(\hat{\theta}^*_1,\dots,\hat{\theta}^*_B\)`. </br> --- # Construction d'intervalles de confiance #### Méthode du bootstrap classique Point de départ : l'existence d'une quantité pivotale, par exemple `\(\hat{\theta} - \theta\)`, de loi `\(H\)`. -- <img src="slides_25_files/figure-html/unnamed-chunk-33-1.png" width="360" /> -- IC : `\(I(1-\alpha) = \big[\hat{\theta} - H^{-1}(1-\alpha/2); \hat{\theta} - H^{-1}(\alpha/2)\big]\)` --- # Construction d'intervalles de confiance Version bootstrap empirique : `\(\hat{I}^*(1-\alpha) = \big[\hat{\theta} - (\hat{H}^*_B)^{-1}(1-\alpha/2); \hat{\theta} - (\hat{H}^*_B)^{-1}(\alpha/2)\big]\)` #### Calcul des quantiles de `\(\hat{H}^*_B\)` - `\(H^*(x) = G^*(x + \hat{\theta})\)` -- <img src="slides_25_files/figure-html/unnamed-chunk-34-1.png" width="504" /> -- - `\((H^*)^{-1}(y) = (G^*)^{-1}(y) - \hat{\theta}\)` --- # Construction d'intervalles de confiance #### Méthode du bootstrap classique `$$\hat{I}_{boot}^*(1-\alpha) = \big[2\hat{\theta} - \hat{\theta}^*_{(\lceil B(1-\alpha/2)\rceil )} ; 2\hat{\theta} - \hat{\theta}^*_{(\lceil B \alpha/2\rceil )} \big]$$` --- # Construction d'intervalles de confiance #### Méthode percentile <span style="color:red">Hypothèse : il existe une transformation `\(h\)` **croissante** telle que la loi de `\(h(\hat{\theta})\)` soit symétrique autour de `\(\eta = h(\theta)\)`</span> IC pour `\(\eta\)` : $$IC_\eta(1-\alpha) = \left[U - H^{-1}_U\left( 1 - \frac{\alpha}{2}\right) ; U - H^{-1}_U \left(\frac{\alpha}{2}\right)\right]. $$ - `\(h\)` inconnue ? - `\(H_U\)` inconnue ? <span style="color:#16A085">**fin du cours 4**</span> --- name: c5 # Construction d'intervalles de confiance Démarche : 1. construire les échantillons bootstrap `\(\mathcal{X}_1^*, \dots, \mathcal{X}_B^*\)` 2. construire les statistiques bootstrapées `\(\hat{\theta}^*_1, \dots, \hat{\theta}^*_B\)` et leurs transformations par `\(h\)` : `\(U_b^* = h(\hat{\theta}^*_b), \forall b\)` 3. définir la fonction de répartition de la statistique bootstrap : `$$H^*_{U} (x) = \mathbb{P}(U_b^* - U \leq x \mid F_n)$$` 4. intervalle de confiance bootstrap pour `\(\eta = h(\theta)\)` : `$$IC_\eta^*(1-\alpha) = \left[U - H^{*-1}_{U}\left(1 -\frac{\alpha}{2} \right) ; U - H^{* \ -1}_{U}\left(\frac{\alpha}{2} \right) \right]$$` --- # Construction d'intervalles de confiance .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-35-1.png" width="360" style="display: block; margin: auto;" /> ] .pull-right[ </br> </br> </br> </br> `\(H^{* \ -1}_{U}\left(\frac{\alpha}{2} \right) = - H^{* \ -1}_{U}\left(1 - \frac{\alpha}{2} \right)\)` ] `$$IC_{\eta}^*(1-\alpha) = \left[U + H^{*-1}_{U}\left(\frac{\alpha}{2} \right) ; U + H^{* \ -1}_{U}\left(1-\frac{\alpha}{2} \right) \right]$$` -- `$$IC_{\eta}^*(1-\alpha) = \left[G^{*-1}_{U}\left(\frac{\alpha}{2} \right) ; G^{* \ -1}_{U}\left(1-\frac{\alpha}{2} \right) \right]$$` --- # Construction d'intervalles de confiance #### Méthode `\(t\)`-percentile Repose sur l'existence d'une quantité pivotale de la forme (de loi notée `\(J_n\)`): `$$S_n = \sqrt{n} \ \frac{\hat{\theta} - \theta}{\hat{\sigma}}$$` -- IC classique : `\(IC(1-\alpha) = \left[\hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} J_n^{-1}\left(1-\frac{\alpha}{2}\right) ; \hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} J_n^{-1}\left(\frac{\alpha}{2}\right) \right]\)` -- Démarche bootstrap : 1. construire les échantillons bootstrap `\(\mathcal{X}_1^*, \dots, \mathcal{X}_B^*\)` 2. construire les statistiques bootstrapées `\(S^*_b = \sqrt{n} \ \frac{\hat{\theta}^*_b - \hat{\theta}}{\hat{\sigma}^*_b}, \quad b=1,\dots,B\)` 3. approcher `\(J_n^{-1}\)` par les quantiles empiriques des statistiques bootstrapées `$$\hat{I}_{t-\text{boot}}^*(1-\alpha) = \left[\hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} S^*_{\big(\lceil B(1-\frac{\alpha}{2})\rceil\big)}; \hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} S^*_{\big(\lceil B(\frac{\alpha}{2})\rceil \big)} \right]$$` --- # Intervalles de confiance bootstrap - IC du bootstrap classique : `$$\hat{I}_{boot}^*(1-\alpha) = \big[2\hat{\theta} - \hat{\theta}^*_{(\lceil B(1-\alpha/2)\rceil )} ; 2\hat{\theta} - \hat{\theta}^*_{(\lceil B \alpha/2\rceil )} \big]$$` - IC du bootstrap percentile : `$$\hat{I}_{perc}^*(1-\alpha) = \big[\hat{\theta}^*_{(\lceil B\alpha/2\rceil )} ; \hat{\theta}^*_{(\lceil B(1- \alpha/2)\rceil )} \big]$$` - IC du bootstrap `\(t\)`-percentile : `$$\hat{I}_{t-\text{boot}}^*(1-\alpha) = \left[\hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} S^*_{(\lceil B(1-\frac{\alpha}{2})\rceil)}; \hat{\theta} - \frac{\hat{\sigma}}{\sqrt{n}} S^*_{(\lceil B(\frac{\alpha}{2})\rceil)} \right]$$` --- # Paramétrique ou non paramétrique ? - Jusqu'à présent, on a décrit des procédures de bootstrap dit **non-paramétrique** -- - En effet, on n'a fait aucune hypothèse sur la loi `\(F\)`, ni exploité sa forme paramétrique -- - Au contraire, on a échantillonné les `\(X_{b,i}^*\)` selon la loi `\(F_n\)`, qui est un estimateur **non-paramétrique** de `\(F\)` </br> -- - Il est possible de faire du bootstrap **paramétrique** -- - Dans ce cas, on va utiliser la forme paramétrique de la loi `\(F\)`, que l'on note `\(F_\theta\)`. -- - On échantillonne alors les `\(X_{b,i}^*\)` selon la loi `\(F_{\hat{\theta}}\)`, qui est un estimateur **paramétrique** de `\(F\)` <!-- <table> --> <!-- <tr style="font-weight:bold"> --> <!-- <td> </td> --> <!-- <td>Bootstrap non paramétrique</td> --> <!-- <td>Bootstrap paramétrique</td> --> <!-- </tr> --> <!-- <tr> --> <!-- <td>hypothèse sur `\(F\)`</td> --> <!-- <td>pas d'hypothèse particulière</td> --> <!-- <td>$F = F_\theta$</td> --> <!-- </tr> --> <!-- <tr> --> <!-- <td>échantillonnage</td> --> <!-- <td>selon loi `\(F_n\)`</td> --> <!-- <td>selon loi `\(F_{\hat{\theta}}\)`</td> --> <!-- </tr> --> <!-- </table> --> </br> -- .center[.red[**RQ : c'est parfois la seule approche possible**]] --- # Tests bootstrap Procédure de test usuelle pour un test de niveau `\(\alpha\)` : -- - choix des hypothèses `\(H_0\)` et `\(H_1\)` -- - choix d'une statistique de test `\(T(\mathcal{X})\)` -- - identification de la loi (éventuellement asymptotique) de `\(T(\mathcal{X})\)` -- - déterminer la région de rejet **ou** calculer la `\(p\)`-valeur du test -- - conclure </br> -- .center[.red[→** quelles sont les difficultés possibles ?**]] --- # Tests bootstrap Procédure de test usuelle pour un test de niveau `\(\alpha\)` : - choix des hypothèses `\(H_0\)` et `\(H_1\)` .green[**→ ok**] - choix d'une statistique de test `\(T(\mathcal{X})\)` .green[**→ ok**] - identification de la loi (éventuellement asymptotique) de `\(T(\mathcal{X})\)` .red[**→ loi inconnue**] - déterminer la région de rejet **ou** calculer la `\(p\)`-valeur du test .red[**→ non calculable**] - conclure --- # Tests bootstrap #### Exemple avec du bootstrap non paramétrique `\(X_1,\dots,X_n\)` i.i.d. de loi `\(F_1\)` et `\(Y_1,\dots,Y_m\)` i.i.d. de loi `\(F_2\)`. `$$H_0 : F_1 = F_2 \quad \text{vs.} \quad H_1 : F_1 \neq F_2$$` -- On suppose que : - l'égalité en loi se traduit par une égalité de certains paramètres -- - l'on dispose d'une statistique de test `\(T(\mathcal{X})\)` -- </br> **Objectif** : construire une version bootstrap de `\(T\)` <span style="color:red">**sous `\(H_0\)`**</span> --- # Tests bootstrap #### Exemple avec du boostrap paramétrique `\(X_1,\dots,X_n\)` i.i.d. de loi `\(F_{\theta,\eta}\)`, et `$$H_0 : \theta = \theta_0 \quad \text{vs.} \quad H_1 : \theta \neq \theta_0$$` → `\(\eta\)` est un *paramètre de nuisance*. </br> -- Exemple du TRV : $$ T(X_1,\dots,X_n) = \frac{L(X_1,\dots,X_n ; \hat{\theta}, \hat{\eta}_1)}{L(X_1,\dots,X_n ; \theta_0, \hat{\eta}_0)},$$ </br> -- <span style="color:red">**→ la loi de `\(T\)` sous `\(H_0\)` n'est pas connue en présence de paramètres de nuisnace**</span> </br> <span style="color:#16A085">**fin du cours 5**</span> --- name: c6 class: inverse, middle, center # II. Méthodes de Monte-Carlo --- # Un peu d'histoire ... - première expérience de type Monte Carlo dûe à Buffon au XVIIIème siècle → l'aiguille de Buffon .center[ <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Buffon_needle.svg/1920px-Buffon_needle.svg.png" width="200" align="center">] -- - les méthodes actuelles sont nées pendant la seconde guerre mondiale au laboratoire américain de Los Alamos → nom de code *Monte-Carlo* .center[ <img src="https://frenchriviera.travel/wp-content/uploads/2018/03/Monte-Carlo-Casino1.jpg" width="300" align="center"> ] --- # Objectif - Question principale : **l'approximation d'intégrales** $$\int h(x) dx $$ - Outil théorique à la base des méthodes de Monte-Carlo : **la loi forte des grands nombres** - Outil pratique nécessaire : **la simulation de variables aléatoires** --- class: inverse, middle, center # 1. Génération de variables aléatoires --- # Génération de variables aléatoires Lois usuelles → disponibles sous la plupart des langages ou logiciels ``` r # Sous R rexp(n=10, rate=1/5) ``` ``` ## [1] 2.3466261 1.9510857 5.9010095 0.4384056 0.2316887 0.4712571 0.5493962 ## [8] 3.1169243 1.4005475 1.2041220 ``` ``` python # Sous Python from scipy.stats import expon expon.rvs(scale=5, size=10) ``` .red[**attention aux conventions utilisées dans chaque langage !**] --- # Méthode de la fonction inverse Repose sur le résultat suivant : > Soit `\(F\)` une f.d.r. et `\(F^{-}\)` son inverse (généralisée). Soit `\(U \sim \mathcal{U}([0,1])\)`. Alors `\(X = F^{-}(U)\)` suit la loi `\(F\)`. -- <img src="slides_25_files/figure-html/unnamed-chunk-39-1.png" width="504" style="display: block; margin: auto;" /> --- # Méthode de la fonction inverse Repose sur le résultat suivant : > Soit `\(F\)` une f.d.r. et `\(F^{-}\)` son inverse (généralisée). Soit `\(U \sim \mathcal{U}([0,1])\)`. Alors `\(X = F^{-}(U)\)` suit la loi `\(F\)`. <img src="slides_25_files/figure-html/unnamed-chunk-40-1.png" width="504" style="display: block; margin: auto;" /> --- # Méthode de la fonction inverse Repose sur le résultat suivant : > Soit `\(F\)` une f.d.r. et `\(F^{-}\)` son inverse (généralisée). Soit `\(U \sim \mathcal{U}([0,1])\)`. Alors `\(X = F^{-}(U)\)` suit la loi `\(F\)`. <img src="slides_25_files/figure-html/unnamed-chunk-41-1.png" width="504" style="display: block; margin: auto;" /> --- # Méthode de la fonction inverse Repose sur le résultat suivant : > Soit `\(F\)` une f.d.r. et `\(F^{-}\)` son inverse (généralisée). Soit `\(U \sim \mathcal{U}([0,1])\)`. Alors `\(X = F^{-}(U)\)` suit la loi `\(F\)`. <img src="slides_25_files/figure-html/unnamed-chunk-42-1.png" width="504" style="display: block; margin: auto;" /> --- # Exemple Loi de Cauchy, de densité `\(f(x) = \frac{1}{\pi(1+x^2)}\)` -- ``` r u <- runif(10000) finv <- function(u){tan(pi*(u-0.5))} x <- finv(u) ``` -- <img src="slides_25_files/figure-html/unnamed-chunk-44-1.png" width="432" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet Proposition : > Soit `\(X\)` une v.a. de densité `\(f\)` et soit `\(g\)` une densité de probabilité et une constante `\(M \geq 1\)` t.q. `\(\forall x, f(x) \leq M g(x)\)`. Alors pour simuler selon la loi `\(f\)` il suffit de : 1. simuler `\(Y \sim g\)` 2. simuler `\(U | Y=y \sim \mathcal{U}([0,Mg(y)])\)` 3. si `\(0 < U < f(Y)\)`, poser `\(X=Y\)`, sinon reprendre l'étape 1. *Preuve* <div class="horizontalgap" style="width:10px"></div> --- # Méthode d'acceptation-rejet Illustration <img src="slides_25_files/figure-html/unnamed-chunk-45-1.png" width="720" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet - Tirage de `\(Y\)` selon la loi `\(g\)` <img src="slides_25_files/figure-html/unnamed-chunk-46-1.png" width="720" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet - Tirage de `\(U\)` conditionnellement à `\(Y=y\)` selon la loi `\(\mathcal{U}([0,Mg(y)])\)` → on rejette <img src="slides_25_files/figure-html/unnamed-chunk-47-1.png" width="720" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet - Tirage de `\(U\)` conditionnellement à `\(Y=y\)` selon la loi `\(\mathcal{U}([0,Mg(y)])\)` → on accepte <img src="slides_25_files/figure-html/unnamed-chunk-48-1.png" width="720" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet A la fin : <img src="slides_25_files/figure-html/unnamed-chunk-49-1.png" width="720" style="display: block; margin: auto;" /> --- # Méthode d'acceptation-rejet Proposition : > Soit `\(X\)` une v.a. de densité `\(f\)` et soit `\(g\)` une densité de probabilité et une constante `\(M \geq 1\)` t.q. `\(\forall x, f(x) \leq M g(x)\)`. Alors pour simuler selon la loi `\(f\)` il suffit de : 1. simuler `\(Y \sim g\)` 2. simuler `\(U | Y=y \sim \mathcal{U}([0,Mg(y)])\)` 3. si `\(0 < U < f(Y)\)`, poser `\(X=Y\)`, sinon reprendre l'étape 1. *Preuve* <div class="horizontalgap" style="width:10px"></div> -- <span style="color:#16A085">**fin du cours 6**</span> --- name: c7 # Retour sur l'acceptation-rejet Rappel de l'algorithme de simulation : > Soit `\(X\)` une v.a. de densité `\(f\)` et soit `\(g\)` une densité de probabilité et une constante `\(M \geq 1\)` t.q. `\(\forall x, f(x) \leq M g(x)\)`. Alors pour simuler selon la loi `\(f\)` il suffit de : 1. simuler `\(Y \sim g\)` 2. simuler `\(U | Y=y \sim \mathcal{U}([0,Mg(y)])\)` 3. si `\(0 < U < f(Y)\)`, poser `\(X=Y\)`, sinon reprendre l'étape 1. -- Remarques : - on a seulement besoin de connaître `\(f\)` à une constante multiplicative près - `\(\forall x, f(x) \leq M g(x) \Rightarrow \text{supp}(f) \subset \text{supp}(g)\)` - probabilité d'accepter un candidat : `\(\frac{1}{M}\)` (influence du choix de `\(g\)`) --- # Autre formulation (plus générale) > Soit `\(X\)` une v.a. de densité `\(f\)`, `\(g\)` une densité de probabilité et `\(M\)` t.q. `\(\forall x, f(x) \leq M g(x)\)`. Alors pour simuler selon la loi `\(f\)` il suffit de : 1. simuler `\(Y \sim g\)` 2. simuler `\(U \sim \mathcal{U}([0,1])\)` 3. si `\(0 < U < f(Y)/Mg(Y)\)`, poser `\(X=Y\)`, sinon reprendre l'étape 1. RQ. : cela marche aussi si on connaît `\(f\)` à une constante multiplicative près. --- # Exemple On veut simuler `\(X\)` selon la loi normale centrée réduite à l'aide d'une loi exponentielle de paramètre 1. 1. par symétrie de la loi normale, il suffit de savoir simuler selon la loi de `\(|X|\)` -- 2. on génère ensuite `\(Z\)`, une Bernoulli de paramètre 1/2 et on pose `\(X = |X|\)` si `\(Z=1\)` et `\(X=-|X|\)` si `\(Z=0\)`. -- 3. il faut trouver une constante `\(M\)` telle que `\(f(x)\leq M g(x)\)` pour tout `\(x\)` avec : - densité cible `\(f(x) = \frac{2}{\sqrt{2\pi}} e^{-x^2/2} \mathbb{1}_{x\geq 0}\)` - densité de la loi exponentielle `\(g(x) = e^{-x} \mathbb{1}_{x>0}\)` RQ : on peut chercher `\(\tilde{M}\)` telle que `\(e^{-x^2/2} \leq \tilde{M} e^{-x}\)` pour `\(x >0\)`. --- class: my-one-page-font # Exemple ``` r n <- 1000 f <- function(x){ifelse(x>0,sqrt(2/pi)*exp(-x^2/2),0)} g <- function(x){ifelse(x>0,exp(-x),0)} M <- sqrt(2*exp(1)/pi) # env. 1.31 -> 1/M = 0.76 U1 <- runif(n) Y <- -log(U1) U2 <- runif(n,min = 0, max = M*g(Y)) absX <- Y[U2<f(Y)] Z <- 2*rbinom(length(absX),size = 1, p=1/2) - 1 X <- Z*absX ``` -- <img src="slides_25_files/figure-html/unnamed-chunk-51-1.png" width="648" style="display: block; margin: auto;" /> --- # Un cas particulier Un cas particulier intéressant : le cas `\(f\)` bornée à support compact. -- - on prend pour `\(g\)` la loi uniforme sur le support de `\(f\)` - et on simule `\(U\)` uniforme sur `\([0,m]\)` où `\(m=\max_x f(x)\)` -- ex. : `\(\displaystyle f(x) = \frac{1}{8} |x^4 - 5x^2 + 4| \ \mathbb{1}_{[-2,2]}(x)\)` <img src="slides_25_files/figure-html/unnamed-chunk-52-1.png" width="576" style="display: block; margin: auto;" /> --- # Un cas particulier Un cas particulier intéressant : le cas `\(f\)` bornée à support compact. - on prend pour `\(g\)` la loi uniforme sur le support de `\(f\)` - et on simule `\(U\)` uniforme sur `\([0,m]\)` où `\(m=\max_x f(x)\)` ex. : `\(\displaystyle f(x) = \frac{1}{8} |x^4 - 5x^2 + 4| \ \mathbb{1}_{[-2,2]}(x)\)` <img src="slides_25_files/figure-html/unnamed-chunk-53-1.png" width="576" style="display: block; margin: auto;" /> --- # Un cas particulier `\(f(x) = \frac{1}{8} |x^4 - 5x^2 + 4| \ \mathbb{1}_{[-2,2]}(x)\)` .pull-left[ ``` r Y <- runif(10000,-2,2) U <- runif(10000,0,0.5) accept <- (U<fdens(Y)) X <- Y[accept] hist(X,breaks=50,freq=F) points(x,fx,type="l",col="red") ``` ``` r mean(accept) ``` ``` ## [1] 0.4954 ``` → la moitié des points simulés est 'perdue' ] .pull-right[ <!-- --> ] --- # Un cas particulier Visuellement : on peut représenter les points acceptés dans le rectangle `\([-2,2] \times [0,1/2]\)` </br> <img src="slides_25_files/figure-html/unnamed-chunk-58-1.png" style="display: block; margin: auto;" /> --- class: inverse, middle, center # 2. Méthodes de Monte-Carlo --- # Monte Carlo classique **Objectif** : calculer une intégrale de la forme `\(\mu = \int h(x) f(x) dx\)` où `\(f\)` est une densité de probabilité. -- > *Définition.* Soit `\(X_1,\dots,X_n\)` un échantillon i.i.d. de loi `\(f\)`. L'estimateur de Monte Carlo de `\(\mu\)` est : `$$\hat{\mu}_n = \frac{1}{n} \sum_{i=1}^n h(X_i)$$` Propriétés : - estimateur sans biais - fortement consistant - IC avec le TCL Remarques : - vitesse de convergence en `\(\sqrt{n}\)`, indépendante de la dimension - ne dépend pas de la régularité de `\(h\)` --- name: c8 # Monte Carlo classique **Exemple :** on veut calculer `\(\int_0^{2\pi} (\sin^2(x) + 2\cos(3x))^2 dx\)` -- <img src="slides_25_files/figure-html/unnamed-chunk-59-1.png" width="504" style="display: block; margin: auto;" /> On peut ré-écrire l'intégrale sous la forme $$ `\begin{eqnarray} I & = & 2\pi \int_{\mathbb{R}} (\sin^2(x) + 2\cos(3x))^2 \frac{1}{2\pi}1_{[0,2\pi]}(x)dx \\ & = & 2\pi \ \mathbb{E}((\sin^2(X) + 2\cos(3X))^2), \quad X \sim \mathcal{U}([0,2\pi]) \end{eqnarray}` $$ --- # Monte Carlo classique **Exemple :** on veut calculer `\(I=2\pi \ \mathbb{E}((\sin^2(X) + 2\cos(3X))^2), \quad X \sim \mathcal{U}([0,2\pi])\)` On peut utiliser l'estimateur de Monte-Carlo suivant : 1. on génère `\(X_1,\dots,X_n\)` i.i.d. de loi uniforme sur `\([0,2\pi]\)` <img src="slides_25_files/figure-html/unnamed-chunk-60-1.png" width="504" style="display: block; margin: auto;" /> --- # Monte Carlo classique **Exemple :** on veut calculer `\(I=2\pi \ \mathbb{E}((\sin^2(X) + 2\cos(3X))^2), \quad X \sim \mathcal{U}([0,2\pi])\)` On peut utiliser l'estimateur de Monte-Carlo suivant : 1. on génère `\(X_1,\dots,X_n\)` i.i.d. de loi uniforme sur `\([0,2\pi]\)` <img src="slides_25_files/figure-html/unnamed-chunk-61-1.png" width="504" style="display: block; margin: auto;" /> --- # Monte Carlo classique **Exemple :** on veut calculer `\(I=2\pi \ \mathbb{E}((\sin^2(X) + 2\cos(3X))^2), \quad X \sim \mathcal{U}([0,2\pi])\)` On peut utiliser l'estimateur de Monte-Carlo suivant : 1. on génère `\(X_1,\dots,X_n\)` i.i.d. de loi uniforme sur `\([0,2\pi]\)` 2. on approche `\(I\)` par `\(\hat{I} = 2 \pi \frac{1}{n}\sum_{i=1}^n (\sin^2(X_i) + 2\cos(3X_i))^2\)` <img src="slides_25_files/figure-html/unnamed-chunk-62-1.png" width="504" style="display: block; margin: auto;" /> --- # Monte Carlo classique <img src="slides_25_files/figure-html/unnamed-chunk-63-1.png" width="504" style="display: block; margin: auto;" /> On trouve : ``` ## [1] "Estimation :13.67753, variance :2.03864, IC = [10.879,16.476]" ``` ``` ## [1] "La vraie valeur est : 14.92257" ``` --- # Monte Carlo classique En augmentant la taille de l'échantillon (`\(n=100000\)`) : ``` r n <- 100000 u <- runif(n,0,2*pi) mi = mean((sin(u)^2+2*cos(3*u))^2)*2*pi vari = var((sin(u)^2+2*cos(3*u))^2)*4*pi*pi/n ``` ``` ## [1] "Estimation : 14.93206, variance :0.00202, IC = [14.844,15.02]" ``` ``` ## [1] "La vraie valeur est : 14.92257" ``` --- # Choix de la loi `\(f\)` **Exemple** : on veut calculer `\(\int_0^1 e^{-x}x^{-a} dx\)` pour `\(0 < a < 1\)`. -- Monte Carlo classique → `\(f\)` densité de la loi uniforme sur `\([0,1]\)` et `\(h(x) = e^{-x}x^{-a}\)`. -- ``` r set.seed(04022025) a <- 1/2 n <- 1000 X <- runif(n) hX <- exp(-X) * X^(-a) print(paste0("Estimation : ",mean(hX),", variance : ",var(hX)/n)) ``` ``` ## [1] "Estimation : 1.51278946224264, variance : 0.00492448440750204" ``` -- Peut-on faire mieux ? --- # Choix de la loi `\(f\)` ``` r set.seed(04022025) a <- 1/2 n <- 1000 Y <- rexp(n) hY <- ifelse(Y < 1,Y^(-a),0) print(paste0("Estimation : ",mean(hY),", variance : ",var(hY)/n)) ``` ``` ## [1] "Estimation : 1.4764602688381, variance : 0.00605753103598533" ``` En simulant selon une loi exponentielle à support dans `\(\mathbb{R}^+\)`, on a "perdu" certains points : ceux qui sont plus grands que 1. Cela correspondà une proportion : ``` r mean(Y < 1) ``` ``` ## [1] 0.622 ``` --- # Choix de la loi `\(f\)` A quoi ressemble la fonction `\(h(x) = e^{-x}x^{-a}\)` ? <img src="slides_25_files/figure-html/unnamed-chunk-70-1.png" width="360" style="display: block; margin: auto;" /> --- # Choix de la loi `\(f\)` A quoi ressemble la fonction `\(h(x) = e^{-x}x^{-a}\)` ? <img src="slides_25_files/figure-html/unnamed-chunk-71-1.png" width="504" style="display: block; margin: auto;" /> → en échantillonnant selon la loi uniforme, on obtient des points répartis de façon **uniforme** sur `\([0,1]\)`. → Peut-on trouver une loi d'échantillonnage qui produise plus de points proches de 0 ? --- # Choix de la loi `\(f\)` Prenons `\(f(x) = (1-a) \ x^{-a} \mathbb{1}_{[0,1]}(x)\)` et `\(\tilde{h}(x) = e^{-x}/ (1-a)\)` -- <img src="slides_25_files/figure-html/unnamed-chunk-72-1.png" width="360" style="display: block; margin: auto;" /> -- ``` r Z <- X^(1/(1-a)) hZ <- exp(-Z) / (1-a) print(paste0("Estimation : ",mean(hZ),", variance : ",var(hZ)/n)) ``` ``` ## [1] "Estimation : 1.51749207512399, variance : 0.000157397528531785" ``` <span style="color:#16A085">**fin du cours 8**</span> --- name: c9 # Echantillonnage préférentiel **Exemple 2** : calculer `\(\mathbb{P}(X>10)\)` pour `\(X \sim \mathcal{E}(1)\)`. -- ``` r 1 - pexp(10) # vraie valeur calculée à l'aide de la fonction de répartition mean(rexp(1000)>10) # estimation par Monte Carlo naïf ``` ``` ## [1] 4.539993e-05 ## [1] 0 ``` -- L'estimateur est nul ... on n'a peut-être pas eu de chance sur ce tirage de 1000 ? Et si on ré-essayait ? --- # Echantillonnage préférentiel On tire 100 fois un échantillon de taille `\(n=1000\)` et on regarde ce que vaut l'estimateur par Monte Carlo sur chacun de ces 100 tirages. ``` r repet_1000 <- sapply(1:100,FUN = function(i){mean(rexp(1000)>10)}) ``` -- <img src="slides_25_files/figure-html/unnamed-chunk-76-1.png" width="288" style="display: block; margin: auto;" /> → l'estimateur vaut presque tout le temps 0 (96 fois sur 100), le reste du temps il vaut 0.001 (dans ces 4 cas, cela veut dire qu'une seule observation sur les 1000 du tirage était plus grande que 10). --- # Echantillonnage préférentiel Que se passe t-il si on augmente la taille des échantillons ? ``` r repet_1000 <- sapply(1:100,FUN = function(i){mean(rexp(1000)>10)}) repet_10000 <- sapply(1:100,FUN = function(i){mean(rexp(10000)>10)}) repet_100000 <- sapply(1:100,FUN = function(i){mean(rexp(100000)>10)}) ``` <img src="slides_25_files/figure-html/unnamed-chunk-78-1.png" width="864" /> → on obtient des résultats "raisonnables" quand `\(n=100 000\)`, mais la variance reste très élevée. --- # Echantillonnage préférentiel En prenant comme loi instrumentale `\(\mathcal{E}(1/10)\)`, on peut ré-écrire : `$$\mathbb{P}(X>10) =\mathbb{E}\left(\mathbf{1}_{Z>10} \ 10 \ e^{-9Z/10}\right) \quad \text{où } Z \sim \mathcal{E}(1/10)$$` En utilisant Monte-Carlo, on approche cette espérance à l'aide de la moyenne empirique d'un échantillon `\(Z_1,\dots,Z_n\)` de taille `\(n\)` de v.a. iid de loi exponentielle de paramètre 1/10. `$$\hat{p} = \frac{1}{n} \sum_{i=1}^n 10 \ e^{-0.9 Z_i} \mathbf{1}_{Z_i > 10}$$` --- # Echantillonnage préférentiel Code correspondant (pour un tirage, donc une estimation) : ``` r Z <- rexp(n,1/10) w <- 10*exp(-9*Z/10) # ou plus généralement : dexp(Z,1)/dexp(Z,1/10) mean(w*(Z>10)) ``` ``` ## [1] 4.907244e-05 ``` -- On fait plusieurs tirages (plusieurs estimations de `\(p\)`) : <img src="slides_25_files/figure-html/unnamed-chunk-81-1.png" width="360" style="display: block; margin: auto;" /> C'est beaucoup mieux qu'avec la loi "naïve" ! --- # Echantillonnage préférentiel Avec la loi `\(\mathcal{E}(1)\)` : <img src="slides_25_files/figure-html/unnamed-chunk-83-1.png" width="864" /> Avec la loi `\(\mathcal{E}(1/10)\)` : <img src="slides_25_files/figure-html/unnamed-chunk-85-1.png" width="864" /> --- # Choix de la loi instrumentale - Calcul de la variance de `\(\tilde{h}_n\)` -- - Choix optimal pour `\(g\)` → celui qui minimise cette variance. -- - Variance finie si le ratio `\(f/g\)` est borné (i.e. queues de distribution de `\(g\)` plus lourdes) - La variance est optimale lorsque `\(g \propto |h|f\)` ... mais cette fonction est inconnue -- - En pratique : choisir `\(g\)` telle que `\(|h|f/g\)` soit (presque) constant, et de variance finie. --- class: my-one-page-font # Choix de la loi instrumentale Que se passe t-il lorsque le ratio `\(f/g\)` n'est pas borné ? -- **Exemple** : on veut approcher `\(\int x^2 e^{-x^2/2} dx = \sqrt{2\pi} \ \mathbb{E}(X^2)\)` avec `\(X \sim \mathcal{N}(0,1)\)`. -- 1. en utilisant comme loi instrumentale la loi `\(\mathcal{N}(0,0.5^2)\)`. <img src="slides_25_files/figure-html/unnamed-chunk-86-1.png" width="576" style="display: block; margin: auto;" /> --- # Choix de la loi instrumentale Code et approximation : ``` r set.seed(10022025) n <- 1000 Z1 <- rnorm(n,0,sd=0.5) w1 <- dnorm(Z1,0,1)/dnorm(Z1,0,0.5) print(paste0("Estimation : ",sqrt(2*pi)*mean(Z1^2*w1),", variance : ",2*pi*var(Z1^2*w1)/n)) print(paste0("IC : [",paste0(round(sqrt(2*pi)*mean(Z1^2*w1)+qnorm(c(0.025,0.975))*sqrt(2*pi*var(Z1^2*w1)/n),3),collapse=","),"]")) ``` ``` ## [1] "Estimation : 1.76571757560692, variance : 0.0666274508629393" ## [1] "IC : [1.26,2.272]" ``` Vraie valeur : `\(\sqrt{2 \pi} =\)` 2.5066283. --- # Choix de la loi instrumentale 2. avec la loi instrumentale `\(\mathcal{N}(0,1.2^2)\)` <img src="slides_25_files/figure-html/unnamed-chunk-88-1.png" width="576" style="display: block; margin: auto;" /> --- # Choix de la loi instrumentale Code et approximation : ``` r set.seed(10022025) n <- 1000 Z2 <- rnorm(n,0,sd=1.2) w2 <- dnorm(Z2,0,1)/dnorm(Z2,0,1.2) print(paste0("Estimation : ",sqrt(2*pi)*mean(Z2^2*w2),", variance : ",2*pi*var(Z2^2*w2)/n)) print(paste0("IC : [",paste0(round(sqrt(2*pi)*mean(Z2^2*w2)+qnorm(c(0.025,0.975))*sqrt(2*pi*var(Z2^2*w2)/n),3),collapse=","),"]")) ``` ``` ## [1] "Estimation : 2.58312378482309, variance : 0.00559442405726649" ## [1] "IC : [2.437,2.73]" ``` On a divisé la variance par 11. Vraie valeur : `\(\sqrt{2 \pi} =\)` 2.5066283. → calcul de la loi Gaussienne optimale. --- class: my-one-page-font # Choix de la loi instrumentale Si le ratio `\(f/g\)` n'est pas borné, certains poids d'importance peuvent être *trop élevés* ce qui donne trop de poids à certaines observations extrêmes. <img src="slides_25_files/figure-html/unnamed-chunk-90-1.png" width="576" style="display: block; margin: auto;" /> <img src="slides_25_files/figure-html/unnamed-chunk-91-1.png" width="576" style="display: block; margin: auto;" /> --- # Version auto-normalisée On peut définir une version *auto-normalisée* de l'estimateur (sous l'hypothèse que `\(g>0\)` dès que `\(f>0\)`): `$$\bar{\mu}_n = \frac{\sum_{i=1}^n h(Z_i) f(Z_i)/g(Z_i)}{\sum_{i=1}^n f(Z_i)/g(Z_i)} = \frac{\sum_{i=1}^n w_i h(Z_i)}{\sum_{i=1}^n w_i}$$` -- Sur notre exemple, cela donne : <img src="slides_25_files/figure-html/unnamed-chunk-92-1.png" width="720" style="display: block; margin: auto;" /> --- # Version auto-normalisée En utilisant la version auto-normalisée, on obtient : - un estimateur **biaisé** - fortement consistant Intérêts/avantages : - peut s'utiliser lorsque `\(f/g\)` est connue à une constante multiplicative près - une variance plus faible .red[**dans certains cas**] - permet également de *simuler* selon la loi `\(f\)` Inconvénients : - hypothèse plus forte sur `\(g\)` (`\(g>0\)` dès que `\(f>0\)`, alors qu'en EP classique l'hypothèse `\(g>0\)` dès que `\(hf>0\)` suffit) - la variance est plus compliquée à calculer et ne peut pas être arbitrairement faible --- # Version auto-normalisée Simulation selon la loi `\(f\)` : → en échantillonnant selon une loi **discrète** à support sur l'ensemble des `\(\{Z_1,\dots,Z_n\}\)`, et où la probabilité de tirer la valeur `\(Z_i\)` est `\(w_i /\sum_j w_j\)`. ``` r Z <- rnorm(10000,0,sd=1.5) w2 <- dnorm(Z,0,1)/dnorm(Z,0,1.5) X <- sample(Z, size = 2000, prob = w2/sum(w2)) ``` <img src="slides_25_files/figure-html/unnamed-chunk-94-1.png" width="720" style="display: block; margin: auto;" /> --- # Version auto-normalisée .red[**Attention au biais introduit par la méthode**] - Il faut une **taille d'échantillon suffisante** car la méthode est seulement **asymptotiquement sans biais** - Pour échantillonner à l'aide des poids `\(w_i\)`, il faut tirer un échantillon plus petit que l'échantillon des `\((Z_1,\dots,Z_n)\)` - On perd le caractère "i.i.d." --- # Version auto-normalisée **Exemple : ** on utilise les poids `\(w_i\)` obtenus après un EP de taille `\(n=1000\)` ou `\(n=10000\)` pour tirer un échantillon de taille 800 selon la loi `\(f\)` ``` r Z <- rnorm(1000,0,sd=1.5) w_1000 <- dnorm(Z,0,1)/dnorm(Z,0,1.5) X_1000 <- sample(Z, size = 800, prob = w_1000/sum(w_1000)) Z <- rnorm(10000,0,sd=1.5) w_10000 <- dnorm(Z,0,1)/dnorm(Z,0,1.5) X_10000 <- sample(Z, size = 800, prob = w_10000/sum(w_10000)) ``` <img src="slides_25_files/figure-html/unnamed-chunk-96-1.png" width="720" style="display: block; margin: auto;" /> <span style="color:#16A085">**fin du cours 9**</span> --- name: c10 # Taille effective de l'échantillon - *effective sampling size* : pour identifier un mauvais choix de `\(g\)` - par exemple, on a `\(\mathbb{E}(f(Z)/g(Z)) = 1\)` pour `\(Z \sim g\)` → on peut comparer la moyenne empirique des `\(w_i\)` à 1 - critère ESS : `$$ESS = \frac{\left( \sum_{i=1} w_i \right)^2}{\sum_{i=1}^n w_i^2} = \frac{1}{\sum_{i=1}^n \bar{w}_i^2},$$` avec `\(\bar{w}_i = \frac{w_i}{\sum_{i=1}^n w_i}\)` le poids d'importance normalisé Sur l'exemple précédent : | Loi | Moy. empirique des poids | ESS | |:---------:|:------------------------:|:---------:| | N(0,0.5²) | 0.9763350 | 91.65902 | | N(0,0.8²) | 0.9968692 | 859.86719 | | N(0,1.2²) | 1.0096987 | 956.92246 | --- class: inverse, middle, center # 4. Réduction de variance --- # Variables antithétiques - **Objectif : ** approcher `\(\mu = \int h(x)f(x)dx\)` - Permet de réduire la variance de l'estimateur en exploitant les propriétés de *symétrie* de la loi d'échantillonnage - Elle s'utilise lorsque la loi considérée est stable par transformation, i.e. s'il existe une transformation `\(\nu\)` telle que si `\(X \sim f\)` alors `\(\nu(X) \sim f\)`. - Quelques exemples de telles lois : uniformes, normales, student, (lois symétriques de façon générale), ... > *Définition.* Soit `\(X_1,\dots,X_n\)` un échantillon i.i.d. de loi `\(f\)`. L'estimateur de `\(\mu\)` par la méthode des variables antithétiques est : `$$\hat{\mu}_a = \frac{1}{n} \sum_{i=1}^n \frac{h(X_i) + h(\nu(X_i))}{2}$$` --- # Variables antithétiques Propriétés : - Sans biais - Fortement consistant - IC par TCL : `$$\sqrt{n} (\hat{\mu}_a - \mu) \longrightarrow \mathcal{N}(0,\sigma_a^2)$$` Variance de l'estimateur : -- `$$\text{Var}(\hat{\mu}_a) = \frac{1}{2n} \text{Var}(h(X))(1+\rho),$$` où `\(\rho = \text{Cov}(h(X),h(\nu(X)))/\text{Var}(h(X))\)` --- # Variables antithétiques Dans quel(s) cas cela fonctionne t-il ? > **Proposition :** si `\(h\)` est monotone et si `\(\nu\)` est décroissante, alors la méthode des variables antithétiques fournit un estimateur de variance plus faible qu'avec l'approche classique. *Preuve* --- # Variables antithétiques **Exemple : ** estimer `\(\int_0^1 \frac{1}{1+x^2}dx\)`. ``` r n <- 1000 # Monte Carlo classique u <- runif(n) h <- function(x){1/(1+x^2)} est1 <- mean(h(u)) var1 <- mean((h(u)-est1)^2)/n # Monte Carlo avec variables antithétiques u1 <- u[1:(n/2)]; u2 <- 1-u1 v <- c(u1,u2) est2 <- mean(h(v)) var2 <- mean((0.5*(h(u1)+h(u2)) - est2)^2)/(n/2) ``` ``` ## [1] "Estimateur 1 (MC classique) : 0.780915 variance : 2.67175e-05" ``` ``` ## [1] "Estimateur 2 (variables antithétiques): 0.784972 variance : 4.041e-07" ``` Vraie valeur : `\(\arctan(1) =\)` 0.7853982. <span style="color:#16A085">**fin du cours 10**</span> --- name: c11 class: inverse, middle, center # 5. Chaînes de Markov --- # Introduction - Quelques éléments théoriques pour comprendre les algos MCMC - Une chaîne de Markov est une séquence de variables aléatoires qui évolue avec le temps, et dont la probabilité de transition ne dépend que de la valeur de la chaîne à l'instant présent - chaîne à espace d'états discret ou continu, à temps discret ou continu → on va se concentrer sur le cas espace d'états continu et temps discret --- # Introduction Pourquoi s'intéresse t-on aux chaînes de Markov ? - pour générer des variables aléatoires *dépendantes* selon une chaîne de Markov convergente dont la loi stationnaire est la loi qui nous intéresse - le théorème ergodique nous permet alors d'utiliser cet échantillon dépendant pour approcher les mêmes quantités qu'avec le MC classique --- # Introduction Différence avec les méthodes de Monte Carlo : la façon dont on explore l'espace. Exemple avec l'approximation du nombre `\(\pi\)`. -- <img src="slides_25_files/figure-html/unnamed-chunk-100-1.png" width="750px" style="display: block; margin: auto;" /> --- class: my-one-page-font # Noyau de transition - chaîne de Markov `\(X_1, X_2, \dots, X_n\)` de noyau de transition `\(K\)` : `$$X_n | X_{n-1}, X_{n-2}, \dots, X_2, X_1 \sim K(X_{n-1},X_{n})$$` -- - chaîne homogène en temps si `\(\mathbb{P}(X_n \in A | X_{n-1} = x) = \mathbb{P}(X_1 \in A | X_0 = x)\)` -- - exemple : marche aléatoire simple `$$X_n = X_{n-1} + \varepsilon_n, \quad \varepsilon \sim \mathcal{N}(0,1)$$` le noyau de transition `\(K(X_{n-1},X_n)\)` est celui de la loi `\(\mathcal{N}(X_{n-1},1)\)`. <img src="slides_25_files/figure-html/unnamed-chunk-101-1.png" width="360" style="display: block; margin: auto;" /> --- # Noyau de transition - Propriété : `\(\forall n > 1, \forall x \in E, \forall A \in \mathcal{E}, \quad K^n(x,A) := \int_E K^{n-1}(y,A) K(x,dy)\)` → pour passer de l'état `\(x\)` à l'ensemble `\(A\)` en `\(n\)` étapes, on passe d'abord de `\(x\)` à `\(y\)` en une étape, puis de `\(y\)` à `\(A\)` en `\(n-1\)` étapes - Plus généralement, équation de Chapman-Kolmogorov : `\(K^{n+m}(x,A) := \int_E K^{m}(y,A) K^n(x,dy)\)` → pour passer de `\(x\)` à `\(A\)` en `\(n+m\)` étapes, il faut passer par un état intermédiaire `\(y\)` à l'étape `\(n\)` --- # Temps d'arrêt et nombre de passages - pour un ensemble `\(A\)`, on définit le *temps d'arrêt* en `\(A\)` par : `\(\tau_A = \inf \{ n | X_n \in A \}\)` → c'est le premier temps d'atteinte de l'ensemble `\(A\)` <img src="slides_25_files/figure-html/unnamed-chunk-102-1.png" width="360" style="display: block; margin: auto;" /> - de même on définit le *nombre de passages* en `\(A\)` par : `$$\eta_A = \sum_{n=1}^{+\infty} \mathbb{1}_{X_n \in A}$$` --- # Mesure invariante `\(\pi\)` est une mesure de probabilité **invariante** pour la chaîne de noyau `\(K\)` si `$$\forall B \in \mathcal{E}, \pi(B) = \int_E K(x,B) \pi(dx)$$` - exemple : `\(X_n = \rho X_{n-1} + \varepsilon_n\)`, avec `\(\varepsilon \sim \mathcal{N}(0,1)\)` -- ``` r rho <- 0.75; X <- rep(0,10000) for (i in 1:10000){ Xold <- ifelse(i==1,rnorm(1,0,1/(1-0.75^2)),X[i-1]) X[i] <- rho*Xold + rnorm(1)} ``` <img src="slides_25_files/figure-html/unnamed-chunk-104-1.png" width="360" style="display: block; margin: auto;" /> --- # Irréductibilité - mesure la sensibilité de la chaîne aux conditions initiales - très important en pratique pour les algorithmes MCMC car apporte une garantie de convergence quel que soit le point de départ - dans le cas espace d'états discret, la chaîne est dite irréductible si tous les états communiquent - définition équivalente : `\(\forall x, \forall A\)` t.q. `\(\mu(A)>0\)`, `\(\mathbb{E}_x(\eta_A) > 0\)` - <span style="color:red">condition suffisante : `\(K(x,\cdot) > 0\)`</span> --- # Apériodicité - Irréductible : tout ensemble de mesure non nulle peut être atteint à partir de n'importe quel point de `\(E\)` en un temps fini -- - Mais la chaîne peut rester bloquée dans un cycle -- → apériodicité -- `\(\Rightarrow\)` irréductibilité et apériodicité assurent que la chaîne peut aller de n'importe quel point `\(x\)` à n'importe quel point `\(y\)` **en une seule itération** --- # Récurrence - Déplacement possible de `\(x\)` à `\(y\)` en une étape, pour tous `\(x\)` et `\(y\)` -- - Mais ces déplacements doivent être possibles **une infinité de fois** → récurrence --- # Résumé - Si une chaîne de Markov est récurrente et admet une loi stationnaire, cette loi est également la loi limite de la chaîne - C'est ce qu'on appelle l'**ergodicité** - Alors, si on peut générer une chaîne de Markov dont la loi stationnaire est la loi cible, on obtient un échantillon .red[**(non i.i.d.)**] selon cette loi → théorème ergodique --- class: inverse, middle, center # 6. Algorithmes MCMC --- # Introduction **Objectif** : évaluer `\(\int h(x) f(x) dx\)` - Idée : construire une chaîne de Markov ergodique, admettant pour loi stationnaire la loi cible `\(f\)` - Utiliser les réalisations de la chaîne pour approcher l'intégrale par une moyenne empirique - Intérêts : - réduction de variance - plus facile à mettre en place - plus adapté en grande dimension - Inconvénients : - vitesse de convergence qui peut être lente - forte corrélation possible dans l'échantillon --- # Algorithme de Metropolis-Hastings - Article fondateur en 1953 par N Metropolis, A Rosenbluth, M Rosenbluth, A Teller, E Teller pour des distributions apparaissant en physique nucléaire <img src="slides_25_files/figure-html/unnamed-chunk-105-1.png" width="800px" style="display: block; margin: auto;" /> - Généralisé en 1970 par W Hastings pour des distributions quelconques <img src="slides_25_files/figure-html/unnamed-chunk-106-1.png" width="150px" style="display: block; margin: auto;" /> --- # Algorithme de Metropolis-Hastings Construction d'une chaîne de Markov selon l'algorithme suivant : 1. **initialisation** : `\(X_0\)` 2. **pour** `\(n=1, \dots, N\)` : - `\(Y_n \sim q(\cdot | X_{n-1})\)` où `\(q\)` est la *loi de proposition* - `\(\begin{equation}X_n = \begin{cases}Y_n & \text{avec probabilité } \alpha(X_{n-1},Y_n) \\ X_{n-1} & \text{avec probabilité } 1-\alpha(X_{n-1},Y_n) \end{cases}\end{equation}\)` où : $$\alpha(x,y) = \min \left(1, \frac{f(y) q(x | y)}{f(x) q(y | x)} \right) $$ RQ : fonctionne si on connaît `\(f\)` et/ou `\(q\)` à une constante près. <span style="color:#16A085">**fin du cours 11**</span> --- name: c12 # Algorithme de Metropolis-Hastings Construction d'une chaîne de Markov selon l'algorithme suivant : 1. **initialisation** : `\(X_0\)` 2. **pour** `\(n=1, \dots, N\)` : - `\(Y_n \sim q(\cdot | X_{n-1})\)` où `\(q\)` est la *loi de proposition* - `\(\begin{equation}X_n = \begin{cases}Y_n & \text{avec probabilité } \alpha(X_{n-1},Y_n) \\ X_{n-1} & \text{avec probabilité } 1-\alpha(X_{n-1},Y_n) \end{cases}\end{equation}\)` où : $$\alpha(x,y) = \min \left(1, \frac{f(y) q(x | y)}{f(x) q(y | x)} \right) $$ --- # Algorithme de Metropolis-Hastings **Exemple** : loi cible `\(f\)` loi `\(\mathrm{Beta}(2.5,8)\)` et loi de proposition `\(q\)` la loi `\(\mathcal{U}([0,1])\)` -- ``` r a=2.5; b=8 N <- 2000 X <- rep(0,N); X[1]=runif(1) # initialisation for (i in 2:N){ Y <- runif(1) # candidat généré selon la loi q alpha <- dbeta(Y,a,b)/dbeta(X[i-1],a,b) U <- runif(1) X[i] <- ifelse(U<alpha,Y,X[i-1])} ``` <img src="slides_25_files/figure-html/unnamed-chunk-108-1.png" width="576" style="display: block; margin: auto;" /> --- # Algorithme de Metropolis-Hastings <!-- --> --- # Algorithme de Metropolis-Hastings Construction d'une chaîne de Markov (**dans** `\(\mathbb{R}^p\)`) selon l'algorithme suivant : 1. **initialisation** : `\(X_0\)` 2. **pour** `\(n=1, \dots, N\)` : - `\(Y_n \sim q(\cdot | X_{n-1})\)` où `\(q\)` est la *loi de proposition* - `\(\begin{equation}X_n = \begin{cases}Y_n & \text{avec probabilité } \alpha(X_{n-1},Y_n) \\ X_{n-1} & \text{avec probabilité } 1-\alpha(X_{n-1},Y_n) \end{cases}\end{equation}\)` où : `\(\alpha(X_{n-1},Y_n) = \min \left(1, \frac{f(Y_n)}{f(X_{n-1})}\frac{ q(X_{n-1} | Y_n)}{q(Y_n | X_{n-1})} \right)\)` Toutes choses égales par ailleurs, `\(\alpha(X_{n-1},Y_n)\)` augmente : - avec le ratio `\(\frac{f(Y_n)}{f(X_{n-1})}\)`, - avec le ratio `\(\frac{q(X_{n-1} | Y_n)}{q(Y_n | X_{n-1})}\)` - <span style="color:red">**un candidat t.q. `\(\frac{f(Y_n)}{q(Y_n | X_{n-1})} > \frac{f(X_{n-1})}{q(X_{n-1} | Y_n)}\)` est toujours accepté**</span> --- # Algorithme de Metropolis-Hastings **Exemple** : loi cible `\(f\)` loi `\(\mathrm{Beta}(2.5,8)\)` et loi de proposition `\(q\)` la loi `\(\mathcal{U}([0.8X_{n-1},1.2X_{n-1}])\)` -- ``` r set.seed(19022024) a=2.5; b=8 N <- 5000 X <- rep(0,N); X[1] <- runif(1) for (i in 2:10){ Y <- runif(1,0.8*X[i-1],1.2*X[i-1]) # candidat généré selon la loi q alpha <- (dbeta(Y,a,b)/dbeta(X[i-1],a,b)) * (dunif(X[i-1],0.8*Y,1.2*Y)/dunif(Y,0.8*X[i-1],1.2*X[i-1])) U <- runif(1) X[i] <- ifelse(U<alpha,Y,X[i-1])} print(X[1:10]) ``` ``` ## [1] 0.7837793 0.7837793 0.8002693 0.7793696 0.7378524 0.6250122 0.6250122 ## [8] 0.6016766 0.5871069 0.5538738 ``` --- # Algorithme de Metropolis-Hastings Etat courant : `\(X_{n-1} =\)` 0.5539, loi instrumentale : `\(\mathcal{U}\)`([0.44312,0.66468]). -- <img src="slides_25_files/figure-html/unnamed-chunk-112-1.png" width="432" style="display: block; margin: auto;" /> -- Candidat à l'itération `\(n\)` : `\(Y_n\)` = 0.6417618 --- # Algorithme de Metropolis-Hastings Etat courant : `\(X_{n-1} =\)` 0.5539, loi instrumentale : `\(\mathcal{U}\)`([0.44312,0.66468]). <img src="slides_25_files/figure-html/unnamed-chunk-113-1.png" width="432" style="display: block; margin: auto;" /> Candidat à l'itération `\(n\)` : `\(Y_n\)` = 0.6417618 `\(\frac{f(Y_n)}{q(Y_n | X_{n-1})}\)` = 0.014588, `\(\frac{f(X_{n-1})}{q(X_{n-1} | Y_{n})}\)` = 0.062954, `\(\frac{f(Y_n)}{f(X_{n-1})}\)` = 0.268499 et `\(\frac{q(X_{n-1} | Y_n)}{q(Y_n | X_{n-1})}\)` = 0.863052 `\(\alpha(X_{n-1},Y_n)\)` = 0.2317289 on génère `\(U \sim \mathcal{U}([0,1])\)` : 0.2655087 --- # EP et Metropolis-Hastings Exploration de l'espace, cible : `\(f(x) = 0.3 \mathcal{N}((-3,-2)^T, I_2) + 0.7\mathcal{N}((5,3)^T, I_2)\)` **Echantillonnage préférentiel** -- <img src="slides_25_files/figure-html/unnamed-chunk-114-1.png" width="432" style="display: block; margin: auto;" /> --- # EP et Metropolis-Hastings Exploration de l'espace, cible : `\(f(x) = 0.3 \mathcal{N}((-3,-2)^T, I_2) + 0.7\mathcal{N}((5,3)^T, I_2)\)` **Echantillonnage préférentiel** <img src="slides_25_files/figure-html/unnamed-chunk-115-1.png" width="432" style="display: block; margin: auto;" /> --- # EP et Metropolis-Hastings Exploration de l'espace, cible : `\(f(x) = 0.3 \mathcal{N}((-3,-2)^T, I_2) + 0.7\mathcal{N}((5,3)^T, I_2)\)` **Metropolis-Hastings** <img src="slides_25_files/figure-html/unnamed-chunk-116-1.gif" style="display: block; margin: auto;" /> --- # EP et Metropolis-Hastings .pull-left[ <!-- --> ] .pull-right[ <!-- --> ] --- # AR et Metropolis-Hastings .pull-left[ Acceptation-rejet : exploration globale <!-- --> ] .pull-right[ Metropolis-Hastings : exploration locale <!-- --> ] --- # Algorithme de Metropolis-Hastings Preuve de convergence : -- 1. noyau de transition de la chaîne -- 2. condition suffisante d'existence d'une loi stationnaire : *detailed-balanced condition* -- 3. condition suffisante pour l'apériodicité -- 4. condition suffisante pour l'irréductibilité -- Choix de la loi `\(q\)` : deux cas particuliers - loi instrumentale indépendante de `\(X_{n-1}\)` - loi instrumentale symétrique <span style="color:#16A085">**fin du cours 12**</span> --- name: c13 # MH à loi instrumentale indépendante Cas où `\(q\)` est indépendant de `\(X_{n-1}\)` → `\(q(x | y) = q(x)\)` - la loi instrumentale ne dépend pas de l'état courant de la chaîne. - voir l'exemple de la slide 133 (loi cible `\(\text{Beta}(2.5,8)\)` et loi de proposition `\(q\)` la loi `\(\mathcal{U}([0,1])\)`) - la probabilité d'acceptation s'écrit : `\(\alpha(X_{n-1},Y_n) = \min \left(1, \frac{f(Y_n)}{f(X_{n-1})}\frac{ q(X_{n-1})}{q(Y_n)} \right)\)` - s'il existe `\(M\)` telle que `\(f\leq Mq\)` alors la chaîne est *uniformément ergodique* (propriété très forte) - comparaison avec l'algorithme d'acceptation-rejet → voir TD --- # MH à loi instrumentale symétrique Cas où `\(q\)` est symétrique en `\(x\)` et `\(y\)` → `\(q(y | x) = q(x | y)\)` - souvent on a `\(q(y | x) = q(|y - x|)\)` → exemple de la marche aléatoire symétrique - voir l'exemple de la slide 136 (loi cible `\(\text{Beta}(2.5,8)\)` et loi de proposition `\(q\)` la loi `\(\mathcal{U}([0.8X_{n-1},1.2X_{n-1}])\)`) - la probabilité d'acceptation se simplifie : `\(\alpha(X_{n-1},Y_n) = \min \left(1, \frac{f(Y_n)}{f(X_{n-1})}\right)\)` - si `\(\text{supp}(f) \subset \text{supp}(q)\)` alors la chaîne est ergodique de mesure invariante `\(f\)` - RQ : si `\(f\)` est à support dans `\(\mathbb{R}\)`, la chaîne ne pourra jamais être uniformément ergodique - .red[**attention au choix de la variance de la marche aléatoire**] --- # MH à marche aléatoire (MH-MA) **Exemple** : loi cible `\(f\)` la loi uniforme sur le disque unité, et `\(q(\cdot \mid x) = \mathcal{U}([x-\delta,x+\delta]^2)\)` pour différent choix de `\(\delta\)` et une chaîne de taille `\(N=1000\)` <img src="slides_25_files/figure-html/unnamed-chunk-121-1.gif" style="display: block; margin: auto;" /> Taux d'acceptation empirique : 98.6% - 97.6% - 85.5% --- # MH à marche aléatoire (MH-MA) **Exemple** : loi cible `\(f\)` la loi normale centrée réduite, et loi instrumentale `\(q\)` la loi normale centrée en `\(X_{n-1}\)` et de variance `\(\sigma^2\)`, pour différentes valeurs de `\(\sigma\)`. ``` r X1 <- X2 <- X3 <- rep(0,2000); for (i in 2:2000){ # chaine 1 Y <- rnorm(1,X1[i-1],0.1) alpha <- dnorm(Y)/dnorm(X1[i-1]) X1[i] <- ifelse(runif(1)<alpha,Y,X1[i-1]) # chaine 2 Y <- rnorm(1,X2[i-1],2) alpha <- dnorm(Y)/dnorm(X2[i-1]) X2[i] <- ifelse(runif(1)<alpha,Y,X2[i-1]) # chaine 3 Y <- rnorm(1,X3[i-1],15) alpha <- dnorm(Y)/dnorm(X3[i-1]) X3[i] <- ifelse(runif(1)<alpha,Y,X3[i-1]) } ``` --- # MH à marche aléatoire (MH-MA) <img src="slides_25_files/figure-html/unnamed-chunk-123-1.png" width="1080" style="display: block; margin: auto;" /> -- <img src="slides_25_files/figure-html/unnamed-chunk-124-1.png" width="1080" style="display: block; margin: auto;" /> --- # MH à marche aléatoire (MH-MA) Pour diagnostiquer une mauvaise exploration : - calculer le **taux d'acceptation** (la proportion de candidats ayant été acceptés) - ne doit être ni trop haut, ni trop bas - sur l'exemple précédent : 0.976, 0.4925 et 0.0875 pour les chaînes 1, 2 et 3 - recommandation pratique : <span style="color:#16A085"> `\(\approx\)` 0.45 en dim 1 et 2, `\(\approx\)` 0.25 si `\(p>2\)`</span> -- - étudier l'**autocorrélation** des réalisations de la chaîne (fonction `acf()` sous R, `plot_acf()` du module `statsmodels.graphics.tsaplots` sous Python) <img src="slides_25_files/figure-html/unnamed-chunk-125-1.png" width="972" style="display: block; margin: auto;" /> --- # De MH à l'échantillonneur de Gibbs - Algorithme générique qui fonctionne pour une grande classe de lois `\(f\)` - Le choix de `\(q\)` ne requiert pas une grande connaissance de `\(f\)` - On peut améliorer l'efficacité des algorithmes (e.g. la convergence) en exploitant les caractéristiques de `\(f\)` → algorithmes utilisant la forme de `\(f\)` et sa géométrie --- # Slice sampler > **Théorème** : simuler `\(X\)` selon la loi `\(f\)` est équivalent à simuler le couple `\((X,U)\)` selon la loi uniforme sur l'ensemble `\(\mathcal{D} = \{(x,u) : 0 \leq u \leq f(x)\}\)` *Preuve* -- - Pour simuler `\(X\)` selon la loi `\(f\)`, on peut alors simuler un point `\((X,U)\)` uniformément **sous le graphe** de la densité `\(f\)` et garder uniquement l'abscisse des points (la loi marginale de `\(X\)`) - C'est déjà ce que l'on fait dans la méthode d'acceptation-rejet en simulant dans une région **plus grande** et en ne gardant que les points **sous la courbe** - L'algorithme du slice sampler permet d'atteindre le même objectif mais de façon plus efficace en proposant un déplacement de type marche aléatoire sur l'ensemble `\(\mathcal{D}\)`. --- # Slice sampler Exemple de la loi `\(\text{Beta}(a,b)\)`. Avec l'acceptation-rejet (voir TD) : - on simule `\(Y \sim \mathcal{U}([0,1])\)` - puis `\(U \sim \mathcal{U}([0,M_{a,b}])\)` - et on garde les couples `\((Y,U)\)` tels que `\(U < f(Y)\)` <img src="slides_25_files/figure-html/unnamed-chunk-126-1.png" width="504" style="display: block; margin: auto;" /> → on perd `\(1-1/M\)` soit environ 59.3% des points ... --- # Slice sampler En pratique : 1. **initialisation** : `\(X_0=x_0\)` 2. **pour** `\(n=1,\dots,N\)` - `\(U_n \sim \mathcal{U}([0,f(x_{n-1})])\)` - `\(X_n \sim \mathcal{U}(\mathcal{D}_n)\)` où `\(\mathcal{D}_n = \{x : f(x) \geq u_n\}\)` <img src="slides_25_files/figure-html/unnamed-chunk-127-1.png" width="500px" style="display: block; margin: auto;" /> *(source : Andrieu et al. 2003)* <span style="color:red">**RQ. : on peut utiliser l'algorithme si `\(f\)` est connue à une constante multiplicative près**</span> --- # Slice sampler **Exemple** : on veut simuler selon la loi `\(\mathcal{N}(0,1)\)`, de densité `\(f(x) \propto e^{-x^2/2}\)` → à l'itération `\(n\)`, on a `\(X_{n-1}=x_{n-1}\)` : - on génère `\(U_n \sim \mathcal{U}([0,\exp(-x_{n-1}^2/2)])\)` - on détermine l'ensemble `\(\mathcal{D}_n = \{x : e^{-x^2/2} > u_n\}\)` -- - on simule `\(X_n \sim \mathcal{U}([-\sqrt{-2 \log(u_n)} ; \sqrt{-2 \log(u_n)}])\)` -- <img src="slides_25_files/figure-html/unnamed-chunk-128-1.png" width="504" style="display: block; margin: auto;" /> --- # Slice sampler **Exemple** : on veut simuler selon la loi `\(\mathcal{N}(0,1)\)`, de densité `\(f(x) \propto e^{-x^2/2}\)` → à l'itération `\(n\)`, on a `\(X_{n-1}=x_{n-1}\)` : - on génère `\(U_n \sim \mathcal{U}([0,\exp(-x_{n-1}^2/2)])\)` - on détermine l'ensemble `\(\mathcal{D}_n = \{x : e^{-x^2/2} > u_n\}\)` - on simule `\(X_n \sim \mathcal{U}([-\sqrt{-2 \log(u_n)} ; \sqrt{-2 \log(u_n)}])\)` <img src="slides_25_files/figure-html/unnamed-chunk-129-1.png" width="504" style="display: block; margin: auto;" /> --- # Slice sampler **Exemple** : on veut simuler selon la loi `\(\mathcal{N}(0,1)\)`, de densité `\(f(x) \propto e^{-x^2/2}\)` → à l'itération `\(n\)`, on a `\(X_{n-1}=x_{n-1}\)` : - on génère `\(U_n \sim \mathcal{U}([0,\exp(-x_{n-1}^2/2)])\)` - on détermine l'ensemble `\(\mathcal{D}_n = \{x : e^{-x^2/2} > u_n\}\)` - on simule `\(X_n \sim \mathcal{U}([-\sqrt{-2 \log(u_n)} ; \sqrt{-2 \log(u_n)}])\)` <img src="slides_25_files/figure-html/unnamed-chunk-130-1.png" width="504" style="display: block; margin: auto;" /> --- # Slice sampler **Exemple** : on veut simuler selon la loi `\(\mathcal{N}(0,1)\)`, de densité `\(f(x) \propto e^{-x^2/2}\)` → à l'itération `\(n\)`, on a `\(X_{n-1}=x_{n-1}\)` : - on génère `\(U_n \sim \mathcal{U}([0,\exp(-x_{n-1}^2/2)])\)` - on détermine l'ensemble `\(\mathcal{D}_n = \{x : e^{-x^2/2} > u_n\}\)` - on simule `\(X_n \sim \mathcal{U}([-\sqrt{-2 \log(u_n)} ; \sqrt{-2 \log(u_n)}])\)` <img src="slides_25_files/figure-html/unnamed-chunk-131-1.png" width="504" style="display: block; margin: auto;" /> --- # Slice sampler <img src="slides_25_files/figure-html/unnamed-chunk-132-1.png" width="500px" style="display: block; margin: auto;" /> .center[*(source : Bayesian Statistics Group)*] .pull-left[ .green[**Avantages**] - tirages selon des lois uniformes - pas d'hypothèses spécifiques sur `\(f\)` - `\(f\)` peut être connue à une constante près ] .pull-right[ .red[**Inconvénients**] - calcul de `\(\mathcal{D}_n\)` peut être délicat - extensions possibles dans ce cas qui peuvent être plus lourdes à implémenter ] --- # Slice sampler Exemple avec la loi `\(\text{Beta}(2,5)\)`, on a `\(\mathcal{D_n} = \{(x,u) : C x (1-x)^4 > u_n \}\)`. ``` r set.seed(0) X <- runif(1) U <- runif(1,0,dbeta(X,5,2)) *Dn <- rootSolve::uniroot.all(function(x){dbeta(x,5,2)-U},lower=0,upper=1) Dn ``` ``` ## [1] 0.4177093 0.9808084 ``` <img src="slides_25_files/figure-html/unnamed-chunk-134-1.png" width="504" style="display: block; margin: auto;" /> --- # Slice sampler ``` r n <- 1000 X <- runif(1) for (i in 2:n){ U <- runif(1,0,dbeta(X[i-1],5,2)) Dn <- uniroot.all(function(x){dbeta(x,5,2)-U[i]},lower=0,upper=1) X <- c(X,runif(1,Dn[1],Dn[2]))} ``` <img src="slides_25_files/figure-html/unnamed-chunk-136-1.gif" style="display: block; margin: auto;" /> <span style="color:#16A085">**fin du cours 13**</span> --- name: c14 # Echantillonneur de Gibbs Idée générale : en dimension `\(p\)` > 1, simuler chaque coordonnée l'une après l'autre -- **Cas `\(p=2\)`** `\((X,Y) \in \mathbb{R}^2\)` 1. **initialisation** : `\((X_0,Y_0) = (x_0,y_0)\)` 2. **pour** `\(n=1, \dots, N\)` : - `\(X_n \sim f_{X|Y}(\cdot | y_{n-1})\)` - `\(Y_n \sim f_{Y|X}(\cdot | x_{n})\)` -- **Exemple** : `\(Z=(X,Y) \sim \mathcal{N}\left(\begin{pmatrix}0 \\0 \end{pmatrix}, \begin{pmatrix}1 & \rho \\ \rho & 1 \end{pmatrix}\right)\)` --- # Echantillonneur de Gibbs ``` r n <- 10000; rho <- 0.8; X <- Y <- rep(0,n); for (i in 2:n){ # mise à jour de la 1ère coordonnée X[i] <- rnorm(1,rho*Y[i-1],sqrt(1-rho^2)) # mise à jour de la 2ème coordonnée Y[i] <- rnorm(1,rho*X[i],sqrt(1-rho^2)) } ``` .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-138-1.png" width="360" /> ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-139-1.png" width="432" /> ] --- # Echantillonneur de Gibbs ``` r n <- 10000; rho <- 0.8; X <- Y <- rep(0,n); for (i in 2:n){ # mise à jour de la 1ère coordonnée X[i] <- rnorm(1,rho*Y[i-1],sqrt(1-rho^2)) # mise à jour de la 2ème coordonnée Y[i] <- rnorm(1,rho*X[i],sqrt(1-rho^2)) } ``` .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-141-1.png" width="360" /> ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-142-1.png" width="360" /> ] --- # Echantillonneur de Gibbs En pratique, comment se passe l'itération `\(n\)` ? <ol> <li> on part du point `\((X_{n-1},Y_{n-1})\)` </ol> <img src="slides_25_files/figure-html/unnamed-chunk-143-1.png" width="576" style="display: block; margin: auto;" /> --- # Echantillonneur de Gibbs En pratique, comment se passe l'itération `\(n\)` ? <ol start="2"> <li> on génère la 1ère coordonnée conditionnellement à `\(Y_{n-1}=y_{n-1}\)` </ol> <img src="slides_25_files/figure-html/unnamed-chunk-144-1.png" width="576" style="display: block; margin: auto;" /> --- # Echantillonneur de Gibbs En pratique, comment se passe l'itération `\(n\)` ? <ol start="3"> <li> on génère la 2ème coordonnée conditionnellement à `\(X_n=x_n\)` </ol> <img src="slides_25_files/figure-html/unnamed-chunk-145-1.png" width="576" style="display: block; margin: auto;" /> --- # Echantillonneur de Gibbs → on ne fait que des mouvements *parallèles aux axes* <img src="slides_25_files/figure-html/unnamed-chunk-146-1.png" width="576" style="display: block; margin: auto;" /> --- # Echantillonneur de Gibbs - A un intérêt lorsque la loi jointe `\(f_{X,Y}\)` ou les lois marginales `\(f_X\)` ou `\(f_Y\)` ne peuvent pas être calculées explicitement -- - La chaîne `\((X_t,Y_t)_t\)` est stationnaire de loi `\(f_{X,Y}\)`, de même que `\((X_t)_t\)` et `\((Y_t)_t\)` sont stationnaires de loi `\(f_X\)` et `\(f_Y\)` resp. -- - Simuler selon les lois conditionnelles **suffit** pour obtenir un échantillon selon la loi jointe (théorème de Hammersley-Clifford) - .blue[**Question :** Est-ce le cas si on remplace 'loi conditionnelle' par 'loi marginale' ?] -- - C'est un cas particulier de l'algorithme de Metropolis-Hastings --- # Echantillonneur de Gibbs **Cas général `\(p>2\)`** `\((X_1,\dots,X_p) \in \mathbb{R}^p\)` de loi jointe `\(f\)` 1. **initialisation** : `\((X_{0,1},\dots,X_{0,p}) = (x_{0,1},\dots,x_{0,p})\)` 2. **pour** `\(n=1, \dots, N\)` : - `\(X_{n,1} \sim f_1(\cdot | x_{n-1,2},\dots,x_{n-1,p})\)` - `\(X_{n,2} \sim f_2(\cdot | x_{n,1},x_{n-1,3},\dots,x_{n-1,p})\)` - `\(\dots\)` - `\(X_{n,p} \sim f_p(\cdot | x_{n,1},\dots,x_{n,p-1})\)` - Les densités `\(f_1,\dots,f_p\)` sont appelées *lois conditionnelles complètes*. - On ne fait que des simulations en dimension 1 --- # Echantillonneur de Gibbs Pourquoi ça marche ? > *Théorème de Hammersley-Clifford.* Si le support de la loi jointe est le produit cartésien des supports des lois marginales, alors : `$$f(x_1,\dots,x_p) \propto \prod_{j=1}^p \frac{f_j(x_j | x_1,\dots,x_{j-1},z_{j+1},\dots,z_p)}{f_j(z_j | x_1,\dots,x_{j-1},z_{j+1},\dots,z_p)}$$` pour tout vecteur `\((z_1,\dots,z_p)\)` dans le support de `\(f\)`. </br> Propriétés de convergence : - Mesure invariante - Irréductibilité et récurrence - Apériodicité --- # Condition de positivité **Exemple :** loi cible `\(f(x,y) = \frac{1}{2} \left(\mathbb{1}_{[-1,0]\times[-1,0]}(x,y) + \mathbb{1}_{[0,1]\times[0,1]}(x,y)\right)\)` -- <img src="slides_25_files/figure-html/unnamed-chunk-147-1.png" width="576" style="display: block; margin: auto;" /> --- # Condition de positivité **Exemple :** loi cible `\(f(x,y) = \frac{1}{2} \left(\mathbb{1}_{[-1,0]\times[-1,0]}(x,y) + \mathbb{1}_{[0,1]\times[0,1]}(x,y)\right)\)` <img src="slides_25_files/figure-html/unnamed-chunk-148-1.png" width="576" style="display: block; margin: auto;" /> --- class: my-one-page-font # Variantes - Gibbs avec mise à jour aléatoire : à chaque itération `\(n\)`, on choisit aléatoirement la coordonnée à mettre à jour - Metropolis-within-Gibbs : si on ne sait pas simuler selon `\(f_k\)`, on peut remplacer par une étape de Metropolis-Hastings : **pour** `\(n=1, \dots, N\)` et la coordonnée `\(k\)` : - `\(Y_{n,k} \sim q_k(\cdot | x_{n,1},\dots,x_{n,k-1},x_{n-1,k},x_{n-1,k+1},\dots,x_{n-1,p})\)` - on pose `\(X_{n,k} = Y_{n,k}\)` avec probabilité `\(\alpha_k(X_{n-1,k},Y_{n,k})\)` et `\(X_{n,k}=X_{n-1,k}\)` avec probabilité `\(1- \alpha_k(X_{n-1,k},Y_{n,k})\)`, avec : `$$\alpha_k(X_{n-1,k},Y_{n,k}) = \min\left[1,\frac{\frac{f_k(Y_{n,k} | X_{n,1},\dots,X_{n,k-1},X_{n-1,k+1},\dots,X_{n-1,p})}{q_k(Y_{nk,} | X_{n,1},\dots,X_{n,k-1},X_{n-1,k},X_{n-1,k+1},\dots,X_{n-1,p})}}{\frac{q_k(X_{n-1,k} | X_{n,1},\dots,X_{n,k-1},Y_{n,k},X_{n-1,k+1},\dots,X_{n-1,p})}{f_k(X_{n-1,k} | X_{n,1},\dots,X_{n,k-1},X_{n-1,k+1},\dots,X_{n-1,p})} } \right]$$` --- # Résumé des épisodes précédents .bold[Objectif :] estimer `\(\mu = \int h(x)f(x)dx\)` avec `\(f\)` densité de probabilité -- **Méthode :** <span style="color:orange"> **obtenir des réalisations `\(X_1,\dots,X_n\)` selon la loi `\(f\)`**</span> et <span style="color:blue">**poser `\(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)\)`**</span> -- <span style="color:orange">**1. Obtenir des réalisations selon la loi `\(f\)`**</span> - *méthode 1* : simuler directement selon `\(f\)` </br> </br> </br> </br> - *méthode 2* : utiliser une loi instrumentale `\(g\)` --- # Résumé des épisodes précédents .bold[Objectif :] estimer `\(\mu = \int h(x)f(x)dx\)` avec `\(f\)` densité de probabilité **Méthode :** <span style="color:orange"> **obtenir des réalisations `\(X_1,\dots,X_n\)` selon la loi `\(f\)`**</span> et <span style="color:blue">**poser `\(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)\)`**</span> <span style="color:orange">**1. Obtenir des réalisations selon la loi `\(f\)`**</span> - *méthode 1* : simuler directement selon `\(f\)` - méthode de la fonction inverse - lien avec d'autres lois (e.g. loi conditionnelle, ...) - slice sampler - échantillonneur de Gibbs (en passant par les lois conditionnelles) - *méthode 2* : utiliser une loi instrumentale `\(g\)` - acceptation-rejet - échantillonnage préférentiel avec ré-échantillonnage - algorithme de Metropolis-Hastings - Metropolis-within-Gibbs --- name: simu_f # Résumé des épisodes précédents .bold[Objectif :] estimer `\(\mu = \int h(x)f(x)dx\)` avec `\(f\)` densité de probabilité **Méthode :** <span style="color:orange"> **obtenir des réalisations `\(X_1,\dots,X_n\)` selon la loi `\(f\)`**</span> et <span style="color:blue">**poser `\(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)\)`**</span> <span style="color:orange">**1. Obtenir des réalisations selon la loi `\(f\)`**</span> - *méthode 1* : simuler directement selon `\(f\)` <span style="color:red">→ `\(f\)` **connue exactement**</span> - méthode de la fonction inverse <span style="color:#16A085">→ **éch. i.i.d.** </span> - lien avec d'autres lois (e.g. loi conditionnelle, ...) <span style="color:#16A085">→ **éch. i.i.d.** </span> - slice sampler <span style="color:#16A085">→ **éch. MC** </span> - échantillonneur de Gibbs <span style="color:#16A085">→ **éch. MC** </span> - *méthode 2* : utiliser une loi instrumentale `\(g\)` <span style="color:red">→ `\(f\)` **connue à une constante près** </span> - acceptation-rejet <span style="color:#16A085">→ **éch. i.i.d.** </span> - échantillonnage préférentiel avec ré-échantillonnage <span style="color:#16A085">→ **éch. i.i.d.** </span> - algorithme de Metropolis-Hastings <span style="color:#16A085">→ **éch. MC** </span> - Metropolis-within-Gibbs <span style="color:#16A085">→ **éch. MC** </span> [retour](#simu_post) --- # Résumé des épisodes précédents .bold[Objectif :] estimer `\(\mu = \int h(x)f(x)dx\)` avec `\(f\)` densité de probabilité **Méthode :** <span style="color:orange"> **obtenir des réalisations `\(X_1,\dots,X_n\)` selon la loi `\(f\)`**</span> et <span style="color:blue">**poser `\(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)\)`**</span> <span style="color:blue">**2. Poser `\(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)\)` </span> ** - quelle est la qualité de cette approximation ? - convergence - biais - **variance** - coût computationnel - pour les échantillons de type chaîne de Markov : - convergence de la chaîne ? - exploration de l'espace ? - sensibilité aux conditions initiales ? --- # Diagnostics de convergence des algos MCMC - Les algorithmes MCMC produisent des chaînes de Markov **ergodiques** dont la loi stationnaire est la loi cible `\(f\)` - Comment être sûr(e)/vérifier que la chaîne est bien dans son régime stationnaire ? - **quand peut-on arrêter la chaîne ?** - **quelle est la précision obtenue?** - Trois types d'outils diagnostics pour étudier : - la convergence de la chaîne vers le régime stationnaire - la convergence des moyennes - la convergence vers un échantillonnage i.i.d. - Ce sont en général des .red[**conditions "nécessaires" mais pas suffisantes**] → permet de détecter un problème, mais ne garantit pas que tout s'est bien passé ! --- # Convergence vers la loi stationnaire (1) - En pratique (et même en théorie !) il est difficile de tester si la chaîne produit des réalisations selon `\(f\)` -- - Pour "tester" la convergence vers le régime stationnaire on peut utiliser des .red[**outils graphiques**] : - analyser le **trace plot** (évolution de la chaîne au cours des itérations) pour détecter les fortes non-stationnarités <img src="slides_25_files/figure-html/unnamed-chunk-149-1.png" width="720" style="display: block; margin: auto;" /> --- # Convergence vers la loi stationnaire (1) - En pratique (et même en théorie !) il est difficile de tester si la chaîne produit des réalisations selon `\(f\)` - Pour "tester" la convergence vers le régime stationnaire on peut utiliser des .red[**outils graphiques**] : - analyser le **trace plot** (évolution de la chaîne au cours des itérations) pour détecter les fortes non-stationnarités - lancer **plusieurs chaînes** (en parallèle) et comparer les résultats <img src="slides_25_files/figure-html/unnamed-chunk-150-1.png" width="720" style="display: block; margin: auto;" /> --- # Convergence vers la loi stationnaire (1) .red[**Attention à l'effet "we only see where we've been"**] → exemple avec loi cible `\(0.7\mathcal{N}(1,1) + 0.2\mathcal{N}(6,1)\)` et MH à loi instrumentale `\(\mathcal{N}(0,1)\)`. .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-152-1.png" width="327.6" style="display: block; margin: auto;" /> ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-153-1.png" width="324" style="display: block; margin: auto;" /> ] .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-154-1.png" width="327.6" style="display: block; margin: auto;" /> ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-155-1.png" width="324" style="display: block; margin: auto;" /> ] --- # Convergence vers la loi stationnaire (2) On peut aussi utiliser des .red[**tests de stationnarité**] : - test de Kolmogorov-Smirnov (fonction `ks.test()` sous R et `scipy.stats.kstest()` sous Python) - attention aux hypothèses du test : .red[**échantillons i.i.d.**] - test de Heidelberg et Welch (fonction `heidel.diag()` du package `coda` sous R) - test de Geweke (fonction `geweke.diag()` dans `coda`) - exploration du support de `\(f\)` → si `\(f\)` connue exactement et si `\(d\)` petit - approximation par sommes de Riemann - estimation de la densité --- class: my-one-page-font # Convergence vers la loi stationnaire (3) - En pratique, on définit aussi une période de chauffe (ou *burn-in*) correspondant aux premières itérations de l'algorithme - Ces itérations sont **supprimées** et on ne garde que les itérations suivantes <img src="slides_25_files/figure-html/unnamed-chunk-156-1.png" width="504" style="display: block; margin: auto;" /> ``` r library(coda) heidel.diag(X) ``` ``` ## ## Stationarity start p-value ## test iteration ## [,1] passed 301 0.0704 ## ## Halfwidth Mean Halfwidth ## test ## [,1] failed -0.107 0.107 ``` --- # Convergence des moyennes empiriques (1) - Si la chaîne a atteint son régime stationnaire, elle doit vérifier le *théorème ergodique* - Autrement dit, toutes les expressions de la forme `\(\frac{1}{n}\sum_{i=1}^n h(X_n)\)` doivent converger vers `\(\mathbb{E}(h(X))\)` avec `\(X \sim f\)`. - En pratique, on vérifie la stabilisation de ces quantités : - en traçant l'évolution de la moyenne cumulée - idem pour les quantiles (fonction `cumuplot()` du package `coda`) - variance inter- et intra (plusieurs chaînes) → fonction `gelman.diag()` du package `coda` .pull-left[ <img src="slides_25_files/figure-html/unnamed-chunk-158-1.png" width="576" style="display: block; margin: auto;" /> ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-159-1.png" width="576" style="display: block; margin: auto;" /> ] --- # Convergence vers un échantillonnage i.i.d. On peut "mesurer" la proximité avec un échantillon i.i.d. : - en calculant l'effective sampling size → fonction `effectiveSamplingSize()` du package `coda` - en traçant l'autocorrélation de la chaîne et en élagant (on enlève le burn-in) ``` r par(mfrow=c(1,2)) acf(X,main="Toute la chaîne") acf(X[500:2000],main="Sans le burn-in") ``` <img src="slides_25_files/figure-html/unnamed-chunk-160-1.png" width="576" style="display: block; margin: auto;" /> <span style="color:#16A085">**fin du cours 14**</span>