class: center, middle, inverse, title-slide .title[ # Statistiques computationnelles ] .author[ ###
Charlotte Baey
] .date[ ###
M1 MA - 2024/2025
] --- <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 </br> **1. Méthodes de ré-échantillonnage** - [Cours 1](#c1) (13/01/2025) - [Cours 2](#c2) (14/01/2025) - [Cours 3](#c3) (20/01/2025) - [Cours 4](#c4) (21/01/2025) - [Cours 5](#c5) (27/01/2025) **2. Méthodes de Monte-Carlo** - [Cours 6](#c6) (28/01/2025) - [Cours 7](#c7) (03/02/2025) - [Cours 8](#c8) (04/02/2025) --- 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]$$` <span style="color:#16A085">**fin du cours 1 (13/01/2025)**</span> --- 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}\)` --> <!--</br> <span style="color:#16A085">**fin du cours 1 (22/01/2024)**</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" /> --- # 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" /> --- # 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 (14/01/2025)**</span> --- name: c3 # 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> <span style="color:#16A085">**fin du cours 3 (20/01/2025)**</span> --- name: c4 # 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 ? --- # 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]$$` <span style="color:#16A085">**fin du cours 4 (21/01/2025)**</span> --- name: c5 # Construction d'intervalles de confiance #### Méthode percentile `$$\hat{I}_{perc}^*(1-\alpha) = \big[\hat{\theta}^*_{(\lceil B\alpha/2\rceil )} ; \hat{\theta}^*_{(\lceil B(1- \alpha/2)\rceil )} \big]$$` </br> </br> </br> </br> </br> --- # 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 (27/01/2025)**</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> -- 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\)`) <span style="color:#16A085">**fin du cours 6 (28/01/2025)**</span> --- name: c7 # 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;" /> --- # 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}\)` --- 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 `\(\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 `\(\mathbb{E}(h(X))\)` est : `$$\hat{h}_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\)` --- # Exemples d'applications #### Approcher la `\(p\)`-valeur d'un test ex. avec le test du rapport de vraisemblance d'un modèle gaussien : `$$H_0 : \quad (\mu,\sigma^2)=(\mu_0,\sigma_0^2)\quad\textrm{ contre }\quad H_1 : \quad (\mu,\sigma^2)\neq(\mu_0,\sigma_0^2).$$` -- On a la statistique de test suivante (voir DS de Stat Math 2024) : `$$2\ln V_n= T_n^2-n\ln T_n^2 + Z^2 + n\ln n -n,$$` avec `\(Z=\sqrt{n}\big(\frac{\hat{\mu}_n-\mu_0}{\sigma_0}\big) \quad \textrm{et}\quad T_n^2=\frac{n\hat{\sigma}^2_n}{\sigma_0^2}\)`. Problème : la loi de la statistique de test n'est pas connue (tabulée). --- # Exemples d'applications Les données : ``` ## [1] 98.23 97.91 98.24 99.00 102.64 103.44 103.81 100.64 100.86 98.79 ## [11] 103.48 98.10 101.11 98.30 96.56 98.77 100.15 101.58 97.51 100.87 ## [21] 101.99 98.59 98.72 103.97 94.75 97.92 102.30 96.88 101.44 100.12 ``` -- .pull-left[ ``` r Tn2 <- (n-1)*var(x)/4 V_obs <- Tn2 - n*log(Tn2) + n*(mean(x)-100)^2/4 + n*log(n) - n print(V_obs) ``` ``` ## [1] 1.600647 ``` ``` r N <- 100000 T <- rchisq(N,n-1) Z <- rnorm(N,0,1) V <- T - n*log(T) + Z^2 + n*log(n) - n p_val <- mean(V > V_obs) print(p_val) ``` ``` ## [1] 0.45853 ``` ] .pull-right[ <!-- --> ] --- # Exemples d'applications On peut quantifier l'erreur d'approximation : .pull-left[ ``` r N <- c(100,1000,2000,5000,10000,20000, 50000,100000,200000,300000) p_val <- rep(0,length(N)) for (i in 1:length(N)){ T <- rchisq(N[i],n-1) Z <- rnorm(N[i],0,1) V <- T-n*log(T)+Z^2+n*log(n)-n p_val[i] <- mean(V > V_obs) } ``` ] .pull-right[ <img src="slides_25_files/figure-html/unnamed-chunk-63-1.png" width="432" /> ] --- # Exemples d'application Plus généralement, les méthodes de Monte Carlo s'utilisent pour : - approcher le niveau ou la puissance d'un test - approcher une intégrale en grande dimension - faire de l'optimisation (e.g. algorithme du recuit simulé, descente de gradient stochastique, ...) - ... <span style="color:#16A085">**fin du cours 7 (03/02/2025)**</span> --- name: c8 # Monte Carlo classique Estimation de `\(\mu = \int h(x) f(x) dx\)` à l'aide d'un échantillon `\(X_1,\dots,X_n\)` i.i.d. de loi `\(f\)` : `$$\hat{h}_n = \frac{1}{n} \sum_{i=1}^n h(X_i),$$` -- - estimateur sans biais, fortement consistant - variance : `$$\text{Var} (\hat{h}_n) = \frac{1}{n} \left(\int h^2(x)f(x)dx - \mu^2\right) = \frac{\sigma_f^2}{n}$$` - `\(\sigma_f^2\)` estimée par : `$$v_n = \frac{1}{n} \sum_{i=1}^n h^2(X_i) - \hat{h}_n^2$$` - IC de niveau `\(1-\alpha\)` : `$$\hat{I}_\mu =\left[\hat{h}_n - q_{1-\alpha/2}^{\mathcal{N}(0,1)} \sqrt{v_n/n} ; \hat{h}_n + q_{1-\alpha/2}^{\mathcal{N}(0,1)} \sqrt{v_n/n} \right]$$` --- # 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-64-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-65-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-66-1.png" width="504" style="display: block; margin: auto;" /> --- # Monte Carlo classique **Exemple :** on veut calculer `\(I=\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-67-1.png" width="504" style="display: block; margin: auto;" /> --- # Monte Carlo classique <img src="slides_25_files/figure-html/unnamed-chunk-68-1.png" width="504" style="display: block; margin: auto;" /> On trouve : ``` ## [1] "Estimation :16.56502, variance :2.22047, IC = [13.644,19.486]" ``` ``` ## [1] "La vraie valeur est : 14.92257" ``` --- class: inverse, middle, center # 3. Echantillonnage préférentiel --- # Echantillonnage préférentiel **Objectif** : toujours d'évaluer `\(\mu = \int h(x) f(x) dx\)` On va **ré-écrire l'intégrale** sous la forme : `$$\mu = \int h(x) f(x) dx = \int h(x) \frac{f(x)}{g(x)} g(x) dx,$$` avec `\(g\)` densité de probabilité, telle que `\((g(x) = 0) \Rightarrow (h(x)f(x) = 0)\)`. -- Pourquoi ? - pas toujours possible de simuler selon la loi `\(f\)`, ou d'utiliser la méthode d'acceptation-rejet - ou on s'intéresse à un évènement de probabilité trop faible - dans certains cas cela permet de réduire la variance de l'estimateur --- # Echantillonnage préférentiel > *Définition.* Soit `\(Z_1,\dots,Z_n\)` un échantillon i.i.d. de loi `\(g\)`. L'estimateur par échantillonnage préférentiel de `\(\mu\)` est : `$$\tilde{h}_n = \frac{1}{n} \sum_{i=1}^n h(Z_i) \frac{f(Z_i)}{g(Z_i)} = \frac{1}{n} \sum_{i=1}^n w_i h(Z_i)$$` Propriétés : (si `\(\text{supp}(f) \subset \text{supp}(g)\)`) - estimateur sans biais - fortement consistant Remarques : - On appelle *poids d'importance* les quantités `\(w_i = \frac{f(Z_i)}{g(Z_i)}\)` - La loi `\(g\)` est appelée *loi instrumentale* --- class: my-one-page-font # Echantillonnage préférentiel **Exemple 1** : 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 ? --- # Echantillonnage préférentiel A quoi ressemble la fonction `\(h(x) = e^{-x}x^{-a}\)` ? <img src="slides_25_files/figure-html/unnamed-chunk-71-1.png" width="360" style="display: block; margin: auto;" /> --- # Echantillonnage préférentiel A quoi ressemble la fonction `\(h(x) = e^{-x}x^{-a}\)` ? <img src="slides_25_files/figure-html/unnamed-chunk-72-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 ? --- class: my-one-page-font # Echantillonnage préférentiel Echantillonnage préférentiel : `\(g(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-73-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" ``` On a diminué la variance par 100 pour une même taille d'échantillon. -- Vraie valeur : 1.49365 (intégration numérique) --- # Echantillonnage préférentiel Avant, on intégrait `\(h(x) = e^{-x}x^{-a}\)` en utilisant la loi uniforme sur `\([0,1]\)`. Maintenant, on intègre `\(\tilde{h}(x) = e^{-x}/(1-a)\)` en utilisant la loi de densité `\((1-a) \ x^{-a} \mathbb{1}_{[0,1]}(x)\)`. <img src="slides_25_files/figure-html/unnamed-chunk-75-1.png" width="864" style="display: block; margin: auto;" /> --- # 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-78-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-80-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] 5.706653e-05 ``` -- On fait plusieurs tirages (plusieurs estimations de `\(p\)`) : <img src="slides_25_files/figure-html/unnamed-chunk-83-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-85-1.png" width="864" /> Avec la loi `\(\mathcal{E}(1/10)\)` : <img src="slides_25_files/figure-html/unnamed-chunk-87-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. <span style="color:#16A085">**fin du cours 8 (04/02/2025)**</span>