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