Statistiques computationnelles

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

# 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

---
# Sommaire

.pull-left[
**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)
  - [Cours 9](#c9) (10/02/2025)
  - [Cours 10](#c10) (11/02/2025)

**3. Méthodes MCMC**
  - [Cours 11](#c11) (24/02/2025)
  - [Cours 12](#c12) (25/02/2025)
]
.pull-right[
**3. Méthodes MCMC (suite)**
  - [Cours 13](#c13) (04/03/2025)
  - [Cours 14](#c14) (10/03/2025)
  
** 4. Statistiques bayésiennes**
  - [Cours 15](#c15) (11/03/2025)
  - Cours 16 (pas de slides) (17/03/2025)
  - [Cours 17](#c16) (18/03/2025)
  - [Cours 18](#c18) (24/03/2025)
  - [Cours 19](#c19) (25/03/2025)
]

---
name: c1
class: inverse, middle, center

# Introduction

---
class: my-one-page-font 
# C'est quoi les statistiques computationnelles ?

- 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, ...

---

# 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$`

- en général, on ne dispose que d'une seule réalisation, pour un `$\omega$` donné

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# Statistique paramétrique

---
# 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}$` ?

``` 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 ?

---
class: my-one-page-font 
# Jackknife

Comment construire de nouveaux échantillons ?

---
# 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)}$$`

- 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]$$`

**fin du cours 1 (13/01/2025)**

---
name: c2
# Résumé et remarques

- 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

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

---
# Bootstrap

---
# Bootstrap

---
# Bootstrap

---
# Bootstrap

---
# Bootstrap

---
# Bootstrap

---
# Bootstrap

---
# 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
```

---
# 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
```

---

# 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[
**on raisonne conditionnellement à `$F_n$`**
- é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[
**on raisonne conditionnellement à `$F_n$`**
- é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
]
]

.center[
`$G$` inconnue `$\longrightarrow$` `$G^*$` connue `$\longrightarrow$` `$\hat{G}^*_B$` approximation bootstrap
]

**fin du cours 2 (14/01/2025)**

---
name: c3
# Bootstrap : principe du plug-in

- on remplace `$F$` par sa version empirique `$F_n$`

---
# 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$`.

**fin du cours 3 (20/01/2025)**

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

Hypothèse : il existe une transformation `$h$` **croissante** telle que la loi de `$h(\hat{\theta})$` soit symétrique autour de `$\eta = h(\theta)$`

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[




`$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]$$`
**fin du cours 4 (21/01/2025)**

---
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]$$`

---
# 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$`

--

- 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$`

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

--

---
# Tests bootstrap

Procédure de test usuelle pour un test de niveau `$\alpha$` :

- choix des hypothèses `$H_0$` et `$H_1$` .green[**&rarr; ok**]

- choix d'une statistique de test `$T(\mathcal{X})$` .green[**&rarr; ok**]

- identification de la loi (éventuellement asymptotique) de `$T(\mathcal{X})$` .red[**&rarr; loi inconnue**]

- déterminer la région de rejet **ou** calculer la `$p$`-valeur du test .red[**&rarr; 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})$`

**Objectif** : construire une version bootstrap de `$T$` **sous `$H_0$`**

---
# 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$$`

&rarr; `$\eta$` est un *paramètre de nuisance*.

--

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)},$$

--

**&rarr; la loi de `$T$` sous `$H_0$` n'est pas connue en présence de paramètres de nuisnace**

**fin du cours 5 (27/01/2025)**

---
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 &rarr; 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 &rarr; 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**

---

# 1. Génération de variables aléatoires

---
# Génération de variables aléatoires

Lois usuelles &rarr; 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)
```

---
#  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$`.

---
#  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$`.

---
#  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$`.

---
#  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$`.

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

---
# 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$`)

**fin du cours 6 (28/01/2025)**

---
name: c7
# Méthode d'acceptation-rejet

Illustration

---
# Méthode d'acceptation-rejet

- Tirage de `$Y$` selon la loi `$g$`

---
# Méthode d'acceptation-rejet

- Tirage de `$U$` conditionnellement à `$Y=y$` selon la loi `$\mathcal{U}([0,Mg(y)])$` &rarr; on rejette

---
# Méthode d'acceptation-rejet

- Tirage de `$U$` conditionnellement à `$Y=y$` selon la loi `$\mathcal{U}([0,Mg(y)])$` &rarr; on accepte

---
# Méthode d'acceptation-rejet

A la fin :

---
# 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
```

---
# 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)$`

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

&rarr; la moitié des points simulés est 'perdue'
]

---
# Un cas particulier

Visuellement : on peut représenter les points acceptés dans le rectangle `$[-2,2] \times [0,1/2]$`

---

# 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
```
--

``` 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[
![](slides_25_files/figure-html/unnamed-chunk-61-1.png)
]

---
# Exemples d'applications

On peut quantifier l'erreur d'approximation :

``` 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, ...)

- ...

**fin du cours 7 (03/02/2025)**

---
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]$`

---
# 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]$`

---
# 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$`

---
# 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"
```

---

# 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{\mu}_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 &rarr; `$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}$` ?

---
# 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;" />
&rarr; en échantillonnant selon la loi uniforme, on obtient des points répartis de façon **uniforme** sur `$[0,1]$`.

&rarr; 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)$`.

---
# 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;" />
&rarr; 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)})
```

&rarr; 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$`) :

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)$` :

---
# Choix de la loi instrumentale

- Calcul de la variance de `$\tilde{h}_n$`

- Choix optimal pour `$g$` &rarr; celui qui minimise cette variance.

**fin du cours 8 (04/02/2025)**

---
name: c9
# Choix de la loi instrumentale

- 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)$`.

---
# 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)$`

---
# 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.

&rarr; 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.

---
# 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 :

---
# 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$` :

&rarr; 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))
```

---
# Version auto-normalisée

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

---
# 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$` &rarr; 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))$`

**fin du cours 9 (10/02/2025)**

---
name: c10
# 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
u <- runif(n)
h <- function(x){1/(1+x^2)}
est1 <- mean(h(u))
var1 <- mean((h(u)-est1)^2)/n
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.

---
# Variables de contrôle

L'objectif est d'utiliser une fonction auxiliaire dont on connaît l'espérance pour réduire la variance de l'estimateur.

> **Définition : ** soit `$h_0$` une fonction t.q. `$\mathbb{E}(h_0(X)) = m$` est connue et de variance finie. Soit `$c \in \mathbb{R}$`. L'estimateur par variable de contrôle est :
`$$\hat{\mu}_c = \frac{1}{n} \sum_{i=1}^n [h(X_i) - c(h_0(X_i) - m)]$$`

Propriétés :
 - Sans biais
 - fortement consistant
 - IC par TCL :
  `$$\sqrt{n} (\hat{\mu}_c - \mu) \longrightarrow \mathcal{N}(0,\sigma_c^2)$$`

---
# Variables de contrôle

Dans quel(s) cas cela fonctionne t-il ?

> **Proposition : ** l'estimateur optimal (au sens de la variance) est obtenu pour `$c^* = \frac{\text{Cov}(h(X),h_0(X))}{\text{Var}(h_0(X))}$`

``` r
n <- 1000
X <- rnorm(n)
hX <- X>2
p1 <- mean(hX)
var1 <- var(hX)/n
h0X <- X>0
copt <- 2*(1-pnorm(2))
p2 <- mean(hX - copt*(h0X-0.5))
var2 <- mean((hX - copt*(h0X-0.5) - p2)^2)/n
```

```
## [1] "Estimateur 1 (MC classique) : 0.022 variance : 2.15375e-05"
```

```
## [1] "Estimateur 2 (variable de contrôle): 0.023502 variance : 2.09642e-05"
```

Vraie valeur : 0.0227501 et `$c^* =$` 0.0455003.

---
# Variables de contrôle

- **En pratique le calcul de `$c^*$` est impossible**
 - On peut l'estimer sur un premier échantillon de taille réduite (compromis entre les simulations supplémentaires nécessaires pour estimer `$c^*$` et le gain en variance)

``` r
n <- 300
X <- rnorm(n)
hX <- X>2
h0X <- X>0
copt_hat <- cov(hX,h0X)/var(h0X)
n <- 1000
X <- rnorm(n)
hX <- X>2
h0X <- X>0
p2b <- mean(hX - copt*(h0X-0.5))
var2b <- mean((hX - copt*(h0X-0.5) - p2)^2)/n
```

Approximation de `$c^*$` : 0.0363636.

```
## [1] "Estimateur : 0.022865 variance : 2.09942e-05"
```

**fin du cours 10 (11/02/2025)**

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

&rarr; 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$`.
--

---
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-107-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)$`

&rarr; 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)$`

&rarr; 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 \}$`
 &rarr; c'est le premier temps d'atteinte de l'ensemble `$A$`
 
<img src="slides_25_files/figure-html/unnamed-chunk-108-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)}
```

---
# 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$`

- condition suffisante : `$K(x,\cdot) > 0$`

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

&rarr; 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**

&rarr; 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

&rarr; 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

- Généralisé en 1970 par W Hastings pour des distributions quelconques

<img src="slides_25_files/figure-html/unnamed-chunk-112-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.

---
# 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 <- 5000
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])}
```

**fin du cours 11 (24/02/2025)**

---
name: c12
# 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})}$`
- **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é**

---
# 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-116-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-117-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

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

---
# MH à loi instrumentale indépendante

Cas où `$q$` est indépendant de `$X_{n-1}$` &rarr; `$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 141 (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 &rarr; voir TD

---
# MH à loi instrumentale symétrique

Cas où `$q$` est symétrique en `$x$` et `$y$` &rarr; `$q(y | x) = q(x | y)$`

- souvent on a `$q(y | x) = q(|y - x|)$`
 &rarr; exemple de la marche aléatoire symétrique
 
 - voir l'exemple de la slide 143 (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 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)

---
# 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.969, 0.505 et 0.0875 pour les chaînes 1, 2 et 3
 - recommandation pratique : `$\approx$` 0.45 en dim 1 et 2, `$\approx$` 0.25 si `$p>2$`
--

- étudier l'**autocorrélation** des réalisations de la chaîne (fonction `acf()` sous R, `plot_acf()` du module `statsmodels.graphics.tsaplots` sous Python)

**fin du cours 12 (25/02/2025)**

---
name: c13

# 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$`

&rarr; 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-122-1.png" width="504" style="display: block; margin: auto;" />

&rarr; 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-123-1.png" width="500px" style="display: block; margin: auto;" />
*(source : Andrieu et al. 2003)*

**RQ. : on peut utiliser l'algorithme si `$f$` est connue à une constante multiplicative près**

---
# Slice sampler

**Exemple** : on veut simuler selon la loi `$\mathcal{N}(0,1)$`, de densité `$f(x) \propto e^{-x^2/2}$`

&rarr; à 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-124-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}$`

---
# Slice sampler

**Exemple** : on veut simuler selon la loi `$\mathcal{N}(0,1)$`, de densité `$f(x) \propto e^{-x^2/2}$`

---
# Slice sampler

**Exemple** : on veut simuler selon la loi `$\mathcal{N}(0,1)$`, de densité `$f(x) \propto e^{-x^2/2}$`

---
# Slice sampler

<img src="slides_25_files/figure-html/unnamed-chunk-128-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-130-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]))
}
```

---
# 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-134-1.png" width="360" />
]
.pull-right[
<img src="slides_25_files/figure-html/unnamed-chunk-135-1.png" width="432" />
]

---
# Echantillonneur de Gibbs

.pull-left[
<img src="slides_25_files/figure-html/unnamed-chunk-137-1.png" width="360" />
]
.pull-right[
<img src="slides_25_files/figure-html/unnamed-chunk-138-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>

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

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

---
# Echantillonneur de Gibbs

&rarr; on ne fait que des mouvements *parallèles aux axes*

---
# 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$`.

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

---
# 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)$`

---
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]$$`

**fin du cours 13 (04/03/2025)**

---
name: c14
# Résumé des épisodes précédents

**Méthode :** **obtenir des réalisations `$X_1,\dots,X_n$` selon la loi `$f$`** et **poser `$\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)$`**

**1. Obtenir des réalisations selon la loi `$f$`**

- *méthode 1* : simuler directement selon `$f$`

- *méthode 2* : utiliser une loi instrumentale `$g$`

---
# Résumé des épisodes précédents

**1. Obtenir des réalisations selon la loi `$f$`**

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

**1. Obtenir des réalisations selon la loi `$f$`**

- *méthode 1* : simuler directement selon `$f$` &rarr; `$f$` **connue exactement**
 - méthode de la fonction inverse &rarr; **éch. i.i.d.** 
 - lien avec d'autres lois (e.g. loi conditionnelle, ...) &rarr; **éch. i.i.d.** 
 - slice sampler &rarr; **éch. MC** 
 - échantillonneur de Gibbs &rarr; **éch. MC**

- *méthode 2* : utiliser une loi instrumentale `$g$` &rarr; `$f$` **connue à une constante près** 
 - acceptation-rejet &rarr; **éch. i.i.d.** 
 - échantillonnage préférentiel avec ré-échantillonnage &rarr; **éch. i.i.d.** 
 - algorithme de Metropolis-Hastings &rarr; **éch. MC** 
 - Metropolis-within-Gibbs &rarr; **éch. MC**

[retour](#simu_post)
---
# Résumé des épisodes précédents

**2. Poser `$\hat{\mu} = \frac{1}{n}\sum_{i=1}^n h(X_i)$` **

- 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**]
&rarr; 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-145-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-146-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"**]
&rarr; 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-148-1.png" width="327.6" style="display: block; margin: auto;" />
]
.pull-right[
<img src="slides_25_files/figure-html/unnamed-chunk-149-1.png" width="324" style="display: block; margin: auto;" />
]
.pull-left[
<img src="slides_25_files/figure-html/unnamed-chunk-150-1.png" width="327.6" style="display: block; margin: auto;" />
]
.pull-right[
<img src="slides_25_files/figure-html/unnamed-chunk-151-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$` &rarr; 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-152-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) &rarr; fonction `gelman.diag()` du package `coda`

.pull-left[
<img src="slides_25_files/figure-html/unnamed-chunk-154-1.png" width="576" style="display: block; margin: auto;" />
]
.pull-right[
<img src="slides_25_files/figure-html/unnamed-chunk-155-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 &rarr; fonction `effectiveSamplingSize()` du package `coda`

- en traçant l'autocorrélation de la chaîne et en élagant

``` r
par(mfrow=c(1,2))
acf(X)
acf(X[500:2000])
```

---
# Estimation de la variance

Sous certaines hypothèses d'ergodicité de la chaîne de Markov, elle vérifie un TCL :

`$$\sqrt{n} (\hat{\mu}_n - \mu) \longrightarrow \mathcal{N}(0,\sigma^2_{c})$$`

- expression de `$\sigma^2_c$`
- estimation de `$\sigma_c^2$`
 - par batch means
 - par méthodes spectrales

---
name: c15
class: inverse, middle, center

# III. Introduction aux statistiques bayésiennes

---
# Le point de vue Bayésien

Combiner notre information **a priori** avec les **observations** pour mettre à jour notre connaissance **a posteriori**

---
# Le point de vue Bayésien

Mise à jour **séquentielle** de la connaissance

---
# Le point de vue Bayésien

- Mixer l'information a priori avec les observations

- Donner plus ou moins de poids à l'a priori en fonction du contexte

- Interprétation différente des résultats de l'inférence
 - en bayésien, les probabilités reflètent la *plausabilité* d'un évènement
 - en fréquentiste, les probabilités s'interprètent comme la *fréquence limite* d'un évènement répété

- Interprétation des tests d'hypothèses beaucoup plus ''naturelle''

---
# Le théorème de Bayes

Thomas Bayes (1702-1761) et indépendamment Pierre-Simon de Laplace (1749-1827)

> **Théorème de Bayes**
>
> Pour deux évènements `$A$` et `$B$` :
$$\mathbb{P}(A|B) = \frac{\mathbb{P}(B|A) \ \mathbb{P}(A)}{\mathbb{P}(B)} $$

---
# Le théorème de Bayes

- Version continue du théorème pour `$X$` et `$Y$` deux v.a. :

$$g(y|x) = \frac{f(x|y)h(y)}{\int f(x|y)h(y) dy} $$

- Bayes et Laplace appliquent ce résultat à la modélisation statistique :
 - en représentant l'*incertitude* entourant `$\theta$` par une *distribution* `$\pi$`

$$p(\theta | x) = \frac{f(x|\theta) \pi(\theta)}{\int f(x|\theta) \pi(\theta) d\theta} $$
 - le paramètre du modèle statistique est donc considéré comme **aléatoire**

---
# Modélisation Bayésienne

`$\mathcal{X} = (X_1,\dots,X_n)$` i.i.d de loi `$f$` qui dépend d'une quantité `$\theta$`

- .orange[**Modèle statistique fréquentiste**]
  - `$\theta$` est un paramètre déterministe
  - `$X_i \sim f(\cdot ; \theta)$` &rarr; `$f(\cdot;\theta)$` est la loi **marginale** des observations `$X_i$`
  - L'inférence se fait à l'aide de la **vraisemblance**

- .green[**Modèle statistique bayésien**]
  - `$\theta$` est une variable aléatoire
  - `$X_i \mid \theta \sim f(\cdot | \theta)$` &rarr; `$f(\cdot |\theta)$` est la loi **conditionnelle** des observations `$X_i$` sachant `$\theta$`
  - On ajoute une loi sur `$\theta$`, appelée **loi a priori** : `$\theta  \sim \pi$`
  - L'inférence se fait à l'aide de la **vraisemblance** .red[**et de la loi a priori**], combinées dans la .green[**distribution a posteriori**]
  - Loi a posteriori : loi conditionnelle de `$\theta$` sachant `$\{\mathcal{X}=x\}$`, notée `$p(\theta \mid x)$`.

---
# Exemple historique (Laplace, 1786)

**Question** : est-ce que la probabilité `$\theta$` d'avoir un garçon est plus grande que celle d'avoir une fille ?

- **données** : naissances à Paris entre 1745 et 1770 : 251 527 garçons et 241 945 filles

&rarr; loi des observations ? `$\dots$`

---
# Exemple historique (Laplace, 1786)

**Question** : est-ce que la probabilité `$\theta$` d'avoir un garçon est plus grande que celle d'avoir une fille ?

- **données** : naissances à Paris entre 1745 et 1770 : 251 527 garçons et 241 945 filles

&rarr; loi des observations ? `$X_i \mid \theta \sim \mathcal{B}(\theta) \Leftrightarrow X \mid \mathcal{B}(n,\theta)$` où `$X_i =1$` si garçon et `$0$` sinon et `$X = \sum_{i=1}^n X_i$`.

- **loi a priori** : `$\dots$`

---
# Exemple historique (Laplace, 1786)

**Question** : est-ce que la probabilité `$\theta$` d'avoir un garçon est plus grande que celle d'avoir une fille ?

- **données** : naissances à Paris entre 1745 et 1770 : 251 527 garçons et 241 945 filles

&rarr; loi des observations ? `$X_i \mid \theta \sim \mathcal{B}(\theta) \Leftrightarrow X \mid \mathcal{B}(n,\theta)$` où `$X_i =1$` si garçon et `$0$` sinon et `$X = \sum_{i=1}^n X_i$`.

- **loi a priori** : la probabilité d'avoir un garçon peut valoir n'importe quoi entre 0 et 1 &rarr; `$\theta \sim \mathcal{U}([0,1])$`

- **loi a posteriori** : `$\dots$`

---
# Exemple historique (Laplace, 1786)

**Question** : est-ce que la probabilité `$\theta$` d'avoir un garçon est plus grande que celle d'avoir une fille ?

- **données** : naissances à Paris entre 1745 et 1770 : 251 527 garçons et 241 945 filles

&rarr; loi des observations ? `$X_i \mid \theta \sim \mathcal{B}(\theta) \Leftrightarrow X \mid \mathcal{B}(n,\theta)$` où `$X_i =1$` si garçon et `$0$` sinon et `$X = \sum_{i=1}^n X_i$`.

- **loi a priori** : la probabilité d'avoir un garçon peut valoir n'importe quoi entre 0 et 1 &rarr; `$\theta \sim \mathcal{U}([0,1])$`

- **loi a posteriori** : `$\theta \mid x \sim \text{Beta}(x+1,n-x+1)$`

- On obtient par le calcul `$\mathbb{P}(\theta \leq 1/2 | (251 527, 241 945)) = 1.15 \cdot 10^{-42}$`

``` r
xg <- 251527
xf <- 241945
pbeta(0.5,251527+1,241945+1)
```

```
## [1] 1.146058e-42
```

---
# Le modèle Beta-binomial

- Exemple : sondage électoral .small[(inspiré de *Bayes Rule! An Introduction to Applied Bayesian Modeling*)]

- **Contexte** : on souhaite évaluer les chances d'une candidate A face à une autre à l'aide d'un sondage.

- On note `$X_1,\dots,X_n$` l'intention de vote de chaque électeur et `$\theta$` la probabilité pour qu'un électeur choisisse A.

- Modèle fréquentiste : `$X_i \sim \mathcal{B}(\theta) \ \Leftrightarrow \ \sum_{i=1}^n X_i =: X \sim \mathcal{B}(n,\theta)$`

- Modèle bayésien :
`\begin{align}
X | \theta & \sim \mathcal{B}(n,\theta) \\
\theta & \sim \pi
\end{align}`

---
# Le modèle Beta-binomial

- Modèle bayésien : `$X | \theta \sim \mathcal{B}(n,\theta), \quad \theta \sim \pi$`
 &rarr; comment choisir `$\pi$` ?

- On a déjà réalisé plusieurs sondages :

<img src="slides_25_files/figure-html/unnamed-chunk-161-1.png" width="360" style="display: block; margin: auto;" />
--

- On aimerait choisir une densité de probabilité sur `$[0,1]$` qui reflète ces informations &rarr; on peut prendre la loi `$\text{Beta}(a,b)$`

---
class: my-one-page-font
# Le modèle Beta-binomial

- La loi Beta a pour densité : `$f(x;a,b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1} (1-x)^{b-1} \mathbb{1}_{[0,1]}(x)$`

- `$X \sim \text{Beta}(a,b)$` : `$\mathbb{E}(X) = \frac{a}{a+b}$`, `$\text{Var}(X) = \frac{ab}{(a+b)^2(a+b+1)}$`

- Si `$a>1,b>1$` on peut aussi paramétrer la loi Beta par son mode `$m = \frac{a-1}{a+b-2}$` et sa concentration `$s=a+b$`

---
# Le modèle Beta-binomial

- Que choisir ici ?
- Mode observé dans les sondages : 0.375 &rarr; `$m=0.375$`

---
# Le modèle Beta-binomial

- Loi a priori : `$\pi(\theta)$` &rarr; loi `$\text{Beta}(a,b)$` avec `$a=$` 7.75 et `$b=$` 12.25 (hyperparamètres)

- Vraisemblance : on fait un nouveau sondage sur `$n=50$` personnes &rarr; `$X | \theta \sim \mathcal{B}(50,\theta)$` et on observe `$X(\omega) = 20$`.

---
# Le modèle Beta-binomial

Loi a posteriori `$\text{Beta}(20+7.75,30+12.25)$`

<img src="slides_25_files/figure-html/unnamed-chunk-165-1.png" width="504" style="display: block; margin: auto;" />
**La vraisemblance est renormalisée de telle sorte à ce que son intégrale soit égale à 1, pour pouvoir la comparer plus facilement aux densités a priori et a posteriori qui sont des densités de probabilités pour `$\theta$`, contrairement à la vraisemblance.**

---
# Le modèle Beta-binomial

- Que se passe t-il si on change de loi a priori ?

---
# Le modèle Beta-binomial

- Que se passe t-il si la taille de l'échantillon augmente ? ou diminue ?

**fin du cours 15 (11/03/2025)**

---
name: c16
# Modèle Beta-binomial

Estimateurs ponctuels : moyenne a posteriori et maximum a posterior

---
# Modèle Beta-binomial

Intervalle de crédibilité à 95%
 - en utilisant les quantiles :

``` r
qbeta(c(0.025,0.975),xobs+a,nobs-xobs+b)
```

```
## [1] 0.2126936 0.5839809
```
 - en utilisant la région de plus forte densité : (librairie `HDInterval` sous R et `PyMC` sous Python)

``` r
library(HDInterval)
hdi(qbeta,0.95,shape1=xobs+a,shape2=nobs-xobs+b)
```

```
##     lower     upper 
## 0.2071049 0.5776388 
## attr(,"credMass")
## [1] 0.95
```

---
# Modèle Beta-binomial

Représentation graphique de l'intervalle de crédibilité

---
# Modèle Beta-binomial

La différence entre l'intervalle de crédibilité basé sur les quantiles et celui basé sur l'HDI est plus visible en cas de distribution asymétrique :

- Intervalle de crédibilité donné par les quantiles : [0.043 ; 0.641]
 - Intervalle de crédibilité donné par HDI : [0.018 ; 0.591]

---
# Modèle Beta-binomial

Test d'hypothèses : `$H_0 : \theta \leq 0.5 \quad \text{vs.} \quad H_1 : \theta > 0.5$`

- risque relatif a posteriori : `$\frac{\mathbb{P}(\theta \leq 0.5 | x)}{\mathbb{P}(\theta > 0.5 | x)} = \frac{\int \mathbb{1}_{t\leq 0.5} p(\theta | x)d\theta}{\int \mathbb{1}_{t > 0.5} p(\theta | x)d\theta}$`

``` r
pbeta(0.5,20+7.75,30+12.25)/(1-pbeta(0.5,20+7.75,30+12.25))
```

```
## [1] 23.96023
```

- risque relatif a priori : `$\frac{\mathbb{P}(\theta \leq 0.5)}{\mathbb{P}(\theta > 0.5)} = \frac{\int \mathbb{1}_{t\leq 0.5} \pi(\theta)d\theta}{\int \mathbb{1}_{t > 0.5} \pi(\theta)d\theta}$`

``` r
pbeta(0.5,7.75,12.25)/(1-pbeta(0.5,7.75,12.25))
```

```
## [1] 5.624588
```

- facteur de Bayes : `$BF = \frac{\mathbb{P}(\theta \leq 0.5 | x)}{\mathbb{P}(\theta > 0.5 | x)}  / \frac{\mathbb{P}(\theta \leq 0.5)}{\mathbb{P}(\theta > 0.5)}$`

``` r
(pbeta(0.5,20+7.75,30+12.25)/(1-pbeta(0.5,20+7.75,30+12.25))) / (pbeta(0.5,7.75,12.25)/(1-pbeta(0.5,7.75,12.25)))
```

```
## [1] 4.259908
```

---
# Le facteur de Bayes

---
# Le facteur de Bayes

Echelle d'aide à la décision (Kass et Raftery, 1995)

| BF | 2 log(BF) | Force de la preuve pour H0 |
|:-------------:|:--------------:|:--------------------------------:|
| 1 < BF < 3 | entre 0 et 2 | Mérite à peine d'être mentionnée |
| 3 < BF < 20 | entre 2 et 6 | Positive |
| 20 < BF < 150 | entre 6 et 10 | Forte |
| BF > 150 | supérieur à 10 | Très forte |

---
class: inverse, middle, center

# Lien avec la théorie de la décision

---
class: my-one-page-font
# Fonction de perte

Fonction `$\ell : \Theta \times \Theta \rightarrow \mathbb{R}^+$` telle que `$\ell(\theta,\theta') = 0 \Leftrightarrow \theta = \theta'$`.

.left-code[
<img src="slides_25_files/figure-html/unnamed-chunk-178-1.png" width="360" style="display: block; margin: auto;" />
]
.right-plot[
- perte quadratique `$\ell(\theta,\theta')=(\theta - \theta')^2$` : très utilisée en pratique
 - .green[.bold[avantages :]] convexe, différentiable
 - .red[.bold[inconvénients :]] donne un poids important aux erreurs extrêmes
- perte absolue `$\ell(\theta,\theta')=|\theta - \theta'|$`
 - .green[.bold[avantages :]] convexe, moins de poids sur les extrêmes
 - .red[.bold[inconvénients :]] non différentiable
- perte de Hubert : compromis entre les deux
`$\small \ell(\theta,\theta')= \frac{1}{2} (\theta - \theta')^2 \mathbf{1}_{\theta-\theta'<k} + (k |\theta - \theta'| - \frac{1}{2}k^2)\mathbf{1}_{\theta - \theta' \geq k}$`
]

---
# Fonction de risque

Fonction `$R : \Theta \rightarrow \mathbb{R}^+$` telle que `$R(\theta,T(\mathcal{X})) = \mathbb{E}_\theta (\ell(\theta,T(\mathcal{X}))) = \int_E \ell(\theta, T(x)) dP_\theta(x).$`

**Exemple :** `$\small X \sim \mathcal{B}(n,\theta)$`, estimateurs :
- `$\hat{\theta}_1 = \frac{X}{n}$` (EMV)
- `$\hat{\theta}_2 = \frac{X + \sqrt{n}/2}{n+\sqrt{n}}$` (estimateur de Bayes minimax)
- `$\hat{\theta}_3 = \frac{1}{2}$` (estimateur constant)
- `$\hat{\theta}_4 = \frac{X+1}{n+2}$` (MAP avec loi a priori uniforme sur `$[0,1]$`)
et perte quadratique.

---
# Fonction de risque

Zoom sur les 3 risques R1, R2 et R4.

Chaque fonction de risque en fonction de `$n$`.

---
# Risque de Bayes

Pour une loi a priori `$\pi$`, `$R_{B} (\pi,T(\mathcal{X})) = \int_\Theta R(\theta,T(\mathcal{X})) \pi(\theta) d\theta$`.

Suite exemple : modèle binomial, loi a priori `$\small \mathcal{U}([0,1])$`, `$\hat{\theta}_1 = \frac{X}{n}$`, `$\hat{\theta}_2 = \frac{X + \sqrt{n}/2}{n+\sqrt{n}}$`, `$\hat{\theta}_3 = \frac{1}{2}$`, `$\hat{\theta}_4 = \frac{X+1}{n+2}$`.

---
# Vue d'ensemble

---
# Principaux résultats

> **Théorème :** Pour toute loi a priori `$\pi$`, on a  `$R_B(\pi) \leq R_M$`

&rarr; le risque de Bayes est toujours plus petit que le risque minimax

> **Théorème :** Un estimateur de Bayes unique à équivalence près, pour une loi a priori `$\pi$`, est admissible.

&rarr; un estimateur de Bayes, c'est-à-dire qui minimise l'espérance de la fonction de perte par rapport à la loi a posteriori, est admissible. .red[.bold[On dispose donc une procédure pour construire des estimateurs admissibles]]

> **Théorème :** Un estimateur admissible et de risque constant est un estimateur minimax.

&rarr; **si on peut construire un estimateur de Bayes de risque constant, alors il sera minimax**

> **Théorème :** Un estimateur de Bayes pour une loi a priori `$\pi$` dont la fonction de risque est majorée par son risque de Bayes est minimax.

**fin du cours 17 (18/03/2025)**

---
name: c18
class: inverse, middle, center

# Choix de la loi a priori

---
# Choix de la loi a priori

Plusieurs cas possibles (combinaisons possibles) :
 1. utiliser l'information dont on dispose (dires d'experts, litérature, ...)

3. utiliser des familles de lois conjuguées

2. utiliser des a priori "impropres"

4. utiliser des a priori "non informatifs"

.red[Dans tous les cas, le choix de la loi a priori doit être .bold[justifié] et doit être fait .bold[en amont], i.e. .bold[avant] d'avoir vu et analysé les données.]

---
# 1. Prise en compte de l'information subjective

- pas toujours facile d'en déduire une **forme paramétrique** précise

- pas toujours facile de **combiner** les différentes sources d'information

&rarr; la loi a priori permet de résumer l'information disponible, mais aussi de prendre en compte l'incertitude entourant cette information

- biais cognitifs et difficulté à évaluer de petits (resp. grands) ordres de grandeur

---
# Quelques exemples de biais cognitifs

- **biais de "proximité"** : on a tendance à accorder une plus grande probabilité à un évènement qui nous concerne ou que nous avons rencontré récemment

- **biais d'ajustement et d'ancrage** : influence de la première impression ou information

- **biais de confiance** : on a tendance à sous-estimer l'intervalle de variation d'une quantité (sous-évaluation de l'incertitude)

- **biais d'équiprobabilité** : on à tendance à attribuer trop souvent la probabilité uniforme à un ensemble d'évènements. Lié aussi au biais existant lorsqu'on oublie la fréquence de base d'un évènement avant d'en évaluer des probabilités liées à cet évènement.

- **biais d'autorité** : attribuer trop de valeur au jugement d'une personne dépositaire d'une autorité

- **biais monétaire** : attribuer une plus grande efficacité aux processus/traitements plus chers

- **biais de genre**, biais socio-culturels

- ...

---
# Construction de la loi a priori

Quelques approches :
 - identifier des quantités clés (moyenne, médiane, quantiles, ...) et en déduire une loi paramétrique

- c'est ce que l'on a fait dans l'exemple Beta-binomial
  - dans ce cas il vaut mieux choisir une loi "plus large" que ce que suggèrent l'information

- évaluer la robustesse des résultats au choix de la loi a priori

- introduire un niveau hiérarchique supplémentaire dans la définition de la loi a priori (*modélisation Bayésienne hiérarchique*)

---
class: my-one-page-font
# 2. Utiliser des lois conjuguées

L'approche précédente n'est pas toujours possible :
 - Pas toujours d'information disponible ou de qualité

- Volonté de rester le plus "objectif" possible, en utilisant juste l'information du modèle

Choisir un a priori conjugué :
 - permet d'utiliser un a priori qui ne dépend que du choix du modèle (i.e. de la vraisemblance)

- limite la part de "subjectivité" au choix des paramètres de la loi a priori (et pas à sa forme paramétrique)

- facilite l'analyse car la loi a posteriori est explicite

- constitue tout de même un choix ! et à ce titre a une influence sur les résultats

- peut être incompatible avec l'information dont on dispose (dans ce cas il ne faut pas l'utiliser !)

---
# Exemples de lois conjuguées

Pour un échantillon `$X_1,\dots,X_n$` iid

- Vraisemblance : `$\mathcal{N}(\theta,\sigma^2)$`, `$\sigma$` connue &rarr; loi conjuguée Gaussienne.

- Vraisemblance `$\mathcal{N}(\mu,1/\theta)$`, `$\mu$` connue &rarr; loi conjuguée Gamma.

- Vraisemblance `$\mathcal{G}(k,\theta)$`, `$k$` connu &rarr; loi conjuguée Gamma.

- Vraisemblance `$\mathcal{B}(\theta)$` &rarr; loi conjuguée Beta.

- Vraisemblance `$\mathcal{P}(\theta)$` &rarr; loi conjuguée Gamma.

- Vraisemblance `$\mathcal{NB}(m,\theta)$` &rarr; loi conjuguée Beta.

---
# Lois conjuguées et famille exponentielle

Une densité `$f$` appartient à la **famille exponentielle** si elle s'écrit sous la forme :

`$$f_\theta(x) = h(x) \exp \left( \langle \theta , T(x) \rangle - A(\theta) \right)$$`

Si la loi des observations appartient à cette famille exponentielle, un choix  naturel de loi a priori conjuguée est :

$$\pi(\theta) = K(\mu,\lambda) \exp\left(\langle \theta , \mu \rangle - \lambda A(\theta)\right) $$

ce qui conduit à la loi a posteriori :

$$p(\theta | x) \propto \exp \left(\langle \theta, \mu + \sum_{i=1}^nT(x_i) \rangle - (\lambda + n) A(\theta) \right) $$
--

**Propriété :** si `$T(X)=X$`, la moyenne a posteriori est une combinaison linéaire de la moyenne a priori et de la moyenne empirique des observations :
`$$\mathbb{E}(\theta | x) = \frac{\mu + n\bar{X}}{\lambda + n}$$`

---
# 3. Loi a priori impropre

- Mesure sur `$\Theta$` qui n'est pas nécessairement une mesure de probabilité (ni même une mesure de masse totale finie)

- Utilisable **uniquement** si la loi a posteriori qui en découle est une mesure de probabilité.

- **Exemple** : `$X \sim \mathcal{N}(\theta,1)$` et `$\pi(\theta) = 1$`.

- .red[**La loi a posteriori peut ne pas être définie pour certaines observations**] &rarr; exemple dans le modèle Beta-binomial si prior impropre `$\text{Beta}(0,0)$`.

---
# 4. Lois a priori non-informatives

- A priori de Laplace &rarr; lois uniformes sur `$\Theta$`

- peut conduire à des a priori impropres
 - pas d'invariance par reparamétrisation

- A priori de Jeffreys &rarr; lien avec l'information de Fisher

- si `$\theta \in \mathbb{R}$` : `$\pi(\theta) \propto I^{1/2}(\theta)$`
 - si `$\theta \in \mathbb{R}^d, d> 1$` : `$\pi(\theta) \propto \det(I(\theta))^{1/2}$`

- **invariant par reparamétrisation**

> **Proposition** : Soit `$\pi$` la mesure a priori de Jeffreys pour `$\theta$`, et soit `$\eta = g(\theta)$`, avec `$g$` un `$\mathcal{C}^1$`-difféomorphisme. Alors la mesure image de `$\pi$` par `$g$` est la mesure a priori de Jeffreys pour `$\eta$`.

**fin du cours 18 (24/03/2025)**

---
name: c19
class: inverse, middle, center

# Identification de la loi a posteriori

---
# Loi a posteriori

**Cas 1 : expression exacte disponible** (famille de lois conjuguées, calcul explicite possible)

A partir de cette expression, on peut chercher à calculer :
 - moyenne, médiane, maximum a posteriori
 - quantiles pour régions de crédibilité
 - HDI
 - facteurs de Bayes

&rarr; pas toujours facile à calculer même à partir de l'expression analytique

**Cas 2 : expression connue à une constante près**

- Impossible de calculer moyenne ou médiane a posteriori ni intervalles de crédibilité, facteurs de Bayes, ...

- Toujours possible de calculer le mode a posteriori

&rarr; **on recourt à des approximations de type Monte-Carlo**

---
name: simu_post
# Simulations selon la loi a posteriori

Quelles méthodes utiliser ? [petit retour en arrière ...](#simu_f)

---
# Simulations selon la loi a posteriori

**Exemple** : exercice 3 de la fiche de TD 4, pour `$\small X | \theta \sim \mathcal{B}(n,\theta)$` et `$\small \theta \sim \mathcal{U}([1/2,1])$`. 
On a :

`$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- méthode de la fonction inverse

``` r
n <- 10; X <- 3; c <- 1-pbeta(1/2,X+1,n-X+1)
p3 <- function(x){ifelse(x<1/2|x>1, 0, dbeta(x,X+1,n-X+1)/c)}
U <- runif(10000) ; ech_finv <- qbeta(c*U+pbeta(1/2,X+1,n-X+1), X+1, n-X+1)
```

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- acceptation-rejet

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- acceptation-rejet &rarr;  `$\small p(\theta|x) \leq M (2 \cdot \mathbb{1}_{[1/2,1]}(\theta))$`

``` r
n <- 10; X <- 3; c <- 1-pbeta(1/2,X+1,n-X+1)
M <- dbeta(1/2,X+1,n-X+1)/c

Y <- runif(10000,1/2,1)
U <- runif(10000,0,M)
ech_ar <- Y[U<p3(Y)]
```

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) \propto \tilde{p}(\theta | x)$` avec `$\tilde{p}(\theta|x) =  \theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- acceptation-rejet : fonctionne aussi si on connaît la loi a posteriori à une constante près

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) \propto \tilde{p}(\theta | x)$` avec `$\tilde{p}(\theta|x) =  \theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- acceptation-rejet : fonctionne aussi si on connaît la loi a posteriori à une constante près &rarr;  `$\small \tilde{p}(\theta|x) \leq M (2 \cdot \mathbb{1}_{[1/2,1]}(\theta))$` avec `$M=1/2^{n+1}$`.

``` r
Mtilde <- 1/2^(n+1)
ptilde3 <- function(x){ifelse(x<1/2|x>1,0,x^(X)*(1-x)^(n-X))}
Y <- runif(10000,1/2,1)
U <- runif(10000,0,1)
ech_ar_tilde <- Y[U<(ptilde3(Y)/Mtilde/dunif(Y,1/2,1))]
```

---
class: my-one-page-font
# Simulations selon la loi a posteriori

- échantillonnage préférentiel

``` r
Z1 <- runif(20000,1/2,1)
w <- p3(Z1)/dunif(Z1,1/2,1)
ech_is <- sample(Z1,prob = w/sum(w),10000)
Z2 <- runif(100000,1/2,1)
w <- p3(Z2)/dunif(Z2,1/2,1)
ech_is_2 <- sample(Z2,prob = w/sum(w),10000)
```

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- Metropolis-Hastings à loi instrumentale indépendante

``` r
chaine1 <- rep(0,10000); chaine1[1] <- runif(1,1/2,1)
chaine2 <- rep(0,10000); chaine2[1] <- runif(1,1/2,1)
for (i in 2:10000){
# loi instrumentale U([1/2,1])
cand1 <- runif(1,1/2,1)
alpha1 <- min(1,p3(cand1)/p3(chaine1[i-1])*dunif(chaine1[i-1])/dunif(cand1))
chaine1[i] <- ifelse(runif(1)<alpha1,cand1,chaine1[i-1])

# loi instrumentale Beta(X+1,n-X+1)
cand2 <- rbeta(1,X+1,n-X+1)
alpha2 <- min(1,p3(cand2)/p3(chaine2[i-1])*dbeta(chaine2[i-1],X+1,n-X+1)/dbeta(cand2,X+1,n-X+1))
chaine2[i] <- ifelse(runif(1)<alpha2,cand2,chaine2[i-1])
}
```

```
## [1] "Taux d'acceptation avec loi instrumentale U([1/2,1]) : 0.2775"
## [1] "Taux d'acceptation avec loi instrumentale Beta(X+1,n-X+1) : 0.1147"
```

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- Metropolis-Hastings à loi instrumentale indépendante

---
# Simulations selon la loi a posteriori

**Exemple** : `$\small p(\theta | x) = \frac{1}{1-F_{B(x+1,n-x+1)}} \frac{\Gamma(n+2)}{\Gamma(x+1)\Gamma(n-x+1)}\theta^x(1-\theta)^{n-x} \mathbb{1}_{[1/2,1]}(\theta)$`

- Metropolis-Hastings à marche aléatoire

``` r
chaine <- rep(0,10000); chaine[1] <- runif(1,1/2,1)
for (i in 2:10000){
cand <- chaine[i-1] + rnorm(1,0,0.1)
alpha <- min(1,p3(cand)/p3(chaine[i-1]))
chaine[i] <- ifelse(runif(1)<alpha,cand,chaine[i-1])}
```

```
## [1] "Taux d'acceptation : 0.4552"
```
<img src="slides_25_files/figure-html/unnamed-chunk-197-1.png" width="720" style="display: block; margin: auto;" />

---
class: my-one-page-font
# Simulations selon la loi a posteriori

Que faire de ces échantillons (i.i.d. ou issus de chaînes de Markov) ?
- calcul de la moyenne, médiane, maximum a posteriori
- calcul d'intervalles de crédibilité, de facteurs de Bayes

**Exemple (suite) :** estimation de la moyenne a posteriori (vraie valeur : 0.5704).

|Méthode | Moyenne a posteriori | Ecart-type | Taille (eff) éch |
|:-------------------------------------|:--------------------:|:----------:|:-------------------------------------------------------:|
|Fonction inverse | 0.5707 | 0.000592 | 10000 |
|Acceptation-rejet (AR) | 0.5708 | 0.001413 | 1778 |
|AR sans la constante | 0.5722 | 0.001436 | 1735 |
|EP de taille 20000 | 0.6289 | 0.000800 | 10000 |
|EP de taille 100000 | 0.5760 | 0.000612 | 10000 |
|MH à proposition indépendante | 0.5701 | 0.001689 | 1256 |
|MH à marche aléatoire U([1/2,1]) | 0.5673 | 0.002239 | 651 |
|MH à marche aléatoire Beta(X+1,n-X+1) | 0.5727 | 0.001865 | 1059 |

---
# Simulations selon la loi a posteriori

Comparaison des résultats de la diapo précédente :
 - la **meilleure méthode** au sens de l'écart-type de l'estimateur est celle de la **fonction inverse**
 - pour cette approche, on a généré 10000 simulations iid de loi uniforme pour un obtenir échantillon iid de taille 10000 selon la loi cible
 - pour l'acceptation-rejet, on a obtenu un échantillon iid de taille env. 1750 selon la loi cible, soit env. 10 fois plus petit que pour la fonction inverse &rarr; l'écart-type est donc env. `$\sqrt{10}$` fois plus grand
 - les **deux échantillons obtenus par EP** ont en fait été générés à partir d'échantillons de taille beaucoup plus grande &rarr; ils peuvent atteindre une variance similaire à la méthode de la fonction inverse mais sont bien **plus coûteux**
 - les **méthodes MCMC**, à cause de la corrélation entre les simulations, conduisent à des **écart-types plus élevés** que dans le cas iid

---
# Simulations selon la loi a posteriori

Pour calculer l'écart-type associé aux méthodes MCMC :
 - soit on "élague" la chaine en ne gardant qu'une valeur toutes les `$m$` itérations, où `$m$` est choisi de telle sorte que l'autocorrélation entre `$X_n$` et `$X_{n-m}$` soit négligeable, et on calcule l'écart-type sur le vecteur obtenu (qui sera donc de taille plus petite)
 - soit on utilise le package `ccoda` et la fonction `summary()` qui renvoie l'écart-type estimé par des méthodes de séries temporelles

``` r
# méthode 1
chaine_elaguee <- chaine[seq(1,length(chaine),by=20)] # d'après le graphe d'autocorrélation 
sd(chaine_elaguee)/sqrt(length(chaine_elaguee))
```

```
## [1] 0.002912311
```

``` r
# méthode 2
library(coda)
summary(mcmc(chaine))
```

```
## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       0.572733       0.060705       0.000607       0.001865 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.5020 0.5247 0.5573 0.6050 0.7233
```