Analyse de données 3

class: center, middle, inverse, title-slide

.title[
# Analyse de données 3
]
.author[
### <font size="5"> Charlotte Baey </font>
]
.date[
### <font size="5"> M2 MIASHS - 2025/2026 </font>
]

---

<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-col {
  width: 50%;
  float: left;
}
.right-col {
  width: 50%;
  float: right;
}

</style>

</style>

# Quelques informations pratiques

### Plan du cours
1. Simulation de variables aléatoires
2. Méthodes de Monte-Carlo
3. Modèles de mélange et algorithme EM

### Organisation

- 1 séance de 2h15 par semaine
- Répartition type : entre 1h et 1h30 de cours, et entre 45 minutes et 1h15 de TP

### Evaluation

- 1 DS intermédiaire d'une durée de 2h (sur une séance de cours)
- 1 projet à rendre (en binôme)
- 1 DS final d'une durée de 3h

.red[**Aucun document autorisé lors des examens.**]

---
# Sommaire

**1. Simulations de variables aléatoires**
  - [Cours 1](#c1) (01/09/2025)
  - [Cours 2](#c2) (08/09/2025)
  - [Cours 3](#c3) (15/09/2025)
  
**2. Méthodes de Monte-Carlo**
  - [Cours 4](#c4) (29/09/2025)
  - [Cours 5](#c5) (06/10/2025)
  - [Cours 6](#c6) (20/10/2025)
  
**3. Modèles de mélange**
  - [Cours 7](#c7) (09/11/2025)
  - [Cours 8](#c8) (01/12/2025)
  - [Cours 9](#c9) (05/12/2025)

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

# Simulation de variables aléatoires
---

# Introduction

- Pour quoi faire ? &rarr; pour simuler une expérience aléatoire et étudier le comportement d'une suite de v.a. de loi donnée

- Pour les lois usuelles, la plupart des logiciels contiennent des fonctions permettant de simuler des v.a. suivant une telle loi

- Pour les autres lois, il faut avoir recours à d'autres méthodes :

- méthode de la fonction inverse

- méthode par transformation

- méthode d'acceptation-rejet

---
# Lois usuelles
Lois usuelles &rarr; disponibles sous la plupart des langages ou logiciels

- Sous R

Fonctions commençant par 'r': `runif`, `rnorm`, `rexp`, `rchisq`, `rbinom`, `rt`, ...

``` r
rexp(n=10, rate=1/5)
```

```
##  [1]  3.1882729  6.7624105  2.7516463  2.7334782 18.7305371  0.4567041
##  [7]  1.3800451  1.5683432  4.6532437  1.3696558
```

- Sous Python

Module `scipy.stats` ou `numpy.random`

``` python
from scipy.stats import expon
expon.rvs(scale=5, size=10)
```

```
## array([6.71306787, 2.86552955, 0.40141877, 4.98535841, 3.00978684,
##        0.20228574, 0.49677082, 9.86906228, 1.41406094, 5.13284106])
```

.red[**attention aux conventions utilisées dans chaque langage !**]

---
# Reproductibilité des résultats

- Les valeurs générées sont *pseudo-aléatoires*, elles dépendent d'une séquence déterministe de nombres
- Un paramètre nommé *seed* permet d'obtenir des suites identiques

``` r
rnorm(5)
set.seed(0);rnorm(5)
set.seed(0);rnorm(5)
```

```
## [1] -1.4622811  0.2460287  2.2215117 -0.2284628  0.4604322
## [1]  1.2629543 -0.3262334  1.3297993  1.2724293  0.4146414
## [1]  1.2629543 -0.3262334  1.3297993  1.2724293  0.4146414
```
--

``` python
import numpy as np
from scipy.stats import norm
np.random.seed(0)
norm.rvs(loc=0, size=5)
np.random.seed(0)
norm.rvs(loc=0, size=5)
```

```
## array([1.76405235, 0.40015721, 0.97873798, 2.2408932 , 1.86755799])
## array([1.76405235, 0.40015721, 0.97873798, 2.2408932 , 1.86755799])
```

---
#  Méthode de la fonction inverse

Repose sur le résultat suivant :
> Soit `$X$` une variable aléatoire admettant une densité `$f$` telle que `$f>0$` p.s.. Alors sa *fonction de répartition* `$F$` admet une fonction réciproque `$F^{-1}$`. 
>De plus, si `$U \sim \mathcal{U}(]0,1[)$`, alors `$X = F^{-1}(U)$` suit la loi `$F$`.

---
#  Méthode de la fonction inverse

---
#  Méthode de la fonction inverse

---
#  Méthode de la fonction inverse

---
# Méthode de la fonction inverse

Dans le cas général, i.e. si `$X$` n'admet pas de densité, on définit l'.orange[**inverse généralisée**] `$F^-$`:

`$$\forall u \in ]0,1[, \quad F^-(u) =  \inf \ \{x \in \mathbb{R}, F(x) \geq u  \}$$`

</br>

On a alors le résultat suivant :

> Soit `$U \sim \mathcal{U}(]0,1[)$`, alors `$X = F^{-}(U)$` suit la loi `$F$`.

</br>*Preuve*

</br>
</br>
**RQ :** La fonction inverse (généralisée) de la f.d.r. s'appelle la .orange[**fonction quantile**].

---
# Méthode de la fonction inverse

Le résultat précédent conduit à l'algorithme suivant.

Soit `$F$` la fonction de répartition de la loi selon laquelle on souhaite générer des réalisations, et `$F^{-}$` son inverse (généralisée). Pour générer un échantillon de taille `$n$` selon la loi `$F$`, on réalise les étapes suivantes :

1. On génère `$n$` réalisations `$u_1,\dots,u_n$` selon la loi `$\mathcal{U}(]0,1[)$`
2. On pose `$x_i = F^{-}(u_i)$` pour `$i=1,\dots,n$`
3. On renvoie l'échantillon `$x_1,\dots,x_n$`

---
# Exemple 1 : loi à densité

Soit `$X$` une v.a. de loi de Cauchy, de densité `$f(x) = \frac{1}{\pi(1+x^2)}$`.

1. Quelle est sa fonction de répartition ? et son inverse ?
2. En déduire un algorithme pour générer un échantillon selon la loi de Cauchy.

``` r
u <- runif(10000)
Finv <- function(u){tan(pi*(u-0.5))}
x <- Finv(u)
```

---
# Exemple 2 : loi discrète

Soit `$X$` une v.a. discrète sur `$\{1,2,3,4,5\}$` dont la loi est donnée par `$p_k=\mathbb{P}(X=k)$`, avec:

`$$p_1=0.15, \quad p_2=0.25, \quad p_3=0.30, \quad p_4=0.20, \quad p_5=0.10$$`

1. Quelle est la fonction de répartition `$F$` de `$X$` ?
2. Identifier l'inverse généralisée de `$F$`.

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

1. > Utiliser la méthode de la fonction inverse pour générer un échantillon de taille `$n=1000$` de loi exponentielle de paramètre 5.
2. > Soit `$U$` une variable aléatoire uniforme sur `$[0,1]$`, et `$0 < p < 1$`. Quelle est la loi de la variable aléatoire `$X = \unicode{x1D7D9}_{U \leq p}$` ? En déduire une méthode pour simuler une variable aléatoire de loi binomiale de paramètres `$(n,p)$`, avec `$n>1$`.
3. > Soit `$X$` une variable aléatoire à valeurs dans `$\{1,\dots,K\}$`. On note `$p_k = \mathbb{P}(X=k).$` Utiliser la méthode de la fonction inverse pour proposer une méthode de simulation selon la loi de `$X$` à partir d'une loi uniforme sur `$[0,1].$`

</br>
</br>
**(Exercice 1 de la fiche 1)**

<span style="color:#16A085">**fin du cours 1**</span>

---
name: c2
# Méthode par transformation

- Repose sur les liens entre lois de probabilité.

- Si une v.a. `$X$` de loi `$f$` peut s'écrire en fonction d'une v.a. `$Y$` dont la loi est facile à simuler, on peut utiliser la relation entre `$X$` et `$Y$` pour générer un échantillon de loi `$f$` à partir d'un échantillon de loi `$g$`.

- La méthode de la fonction inverse est un cas particulier.

- Exemples :
 - transformation affines
 - sommes de v.a. i.i.d.
 - loi conditionnelle
 - loi mélange
 
---
# Transformations affines

`$X$` de loi `$f$` est une transformation affine de `$Y$` de loi `$g$`, s'il existe `$a\neq 0$` et `$b$` dans `$\mathbb{R}$` tq :
$$ X = aY + b $$
La loi de `$X$` est alors :
$$  f(x) = \frac{1}{a}g\left(\frac{x-b}{a}\right) $$
</br>

Algorithme de simulation selon `$f$`:
> 1. Générer `$y_1,\dots,y_n$` selon la loi `$g$`.

> 2. Poser `$x_i=ay_i+b$` pour `$i=1,\dots,n$`.

</br>
Quelles lois sont concernées ?

---
# Paramètres de localisation et d'échelle

On appelle .orange[**location-scale family**] une famille de distributions telles que, si `$Y$` a une loi appartenant à cette famille, alors la loi de `$X=aX+b$` appartient également à cette même famille.

Les paramètres `$a$` et `$b$` s'appellent respectivement .orange[**paramètre d'échelle**] et .orange[**paramètre de localisation**].

Quelques exemples :
 - les lois normales
 - les lois de Student
 - les lois de Fréchet
 - les lois de Cauchy
 - les lois uniformes
 - les lois de Laplace
 - ...

---
# Exemple : les lois normales

Proposer un algorithme pour générer un échantillon selon la loi normale de moyenne `$\mu$` et de variance `$\sigma^2$` à partir d'un échantillon de loi normale centrée réduite.

``` r
mu <- 5
sigma <- 2
x <- rnorm(2000)
y <- sigma*x + mu
```

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

1. > Montrer que la famille des lois exponentielles `$\{\mathcal{E}(\lambda),\lambda\in\mathbb{R}^{+*}\}$` forme une famille de lois paramétrées par un paramètre de localisation et d'échelle. Donner la valeur de ces deux paramètres pour la distribution de paramètre `$\lambda$`.

> Quelle transformation affine permet d'obtenir une variable aléatoire `$X$` de loi `$\mathcal{E}(\lambda)$` à partir d'une variable aléatoire `$Y$` de loi `$\mathcal{E}(1)$` ?

2. > Proposer un algorithme pour tirer `$X\sim\mathcal{N}(\mu,\Sigma)$`, où `$X \in \mathbb{R}^d$`, `$\mu$` est le vecteur d'espérances et `$\Sigma$` la matrice de variance-covariance.

</br>
**(Exercice 2 de la fiche 1)**

---
# Lois dérivants d'autres lois

Compléter le tableau suivant donnant le lien entre `$X_1,X_2,\dots$` i.i.d. de loi `$g$` et `$Y$` de loi `$f$`:

|f | g |
|-|-|
|$$\chi^2_n$$ | `$$\mathcal{N}(0,1)$$` |  
|$$t_n $$ | `$$\mathcal{N}(0,1)$$`  | 
|$$\mathcal{F}_{n,p} $$ | `$$\mathcal{N}(0,1)$$` | 
|$$\chi^2_{2n}$$ | `$$\mathcal{E}(1)$$` | 
|$$\Gamma(a,\beta)$$ | `$$\mathcal{E}(1)$$` |
|$$\mathcal{B}e(a,b)$$ | `$$\mathcal{E}(1)$$` |

---
# Lois dérivants d'autres lois

Compléter le tableau suivant donnant le lien entre  `$X_1,X_2,\dots$` i.i.d. de loi `$g$` et `$Y$` de loi `$f$`:

|f | g | lien|
|--|
|$$\chi^2_n$$ | `$$\mathcal{N}(0,1)$$` | `$Y = \sum_{i=1}^nX_i^2 \quad n \in \mathbb{N}^*$`|
|$$t_n $$ | `$$\mathcal{N}(0,1)$$`  | `$Y = \frac{X_{n+1}}{\sqrt{\sum_{i=1}^nX_i^2 /n}}\quad n \in \mathbb{N}^*$`|
|$$\mathcal{F}_{n,p} $$ | `$$\mathcal{N}(0,1)$$` | `$Y = \frac{\sum_{i=1}^nX_i^2/n}{\sum_{i=n+1}^{n+p}X_i^2/p}$`|
|$$\chi^2_{2n}$$ | `$$\mathcal{E}(1)$$` | `$Y = 2 \sum_{i=1}^n X_i \quad n \in \mathbb{N}^*$`|
|$$\Gamma(a,\beta)$$ | `$$\mathcal{E}(1)$$` | `$Y = \beta \sum_{i=1}^aX_i^2 \quad a \in \mathbb{N}^*$`|
|$$\mathcal{B}e(a,b)$$ | `$$\mathcal{E}(1)$$` | `$Y = \frac{\sum_{j=1}^aX_j}{\sum_{j=1}^{a+b}X_j} \quad a,b \in \mathbb{N}^*$`|

---
# Méthode par transformation

On peut également utiliser les relations entre lois jointes, lois conditionnelles et lois marginales.

- **Exemple 1 :**
La loi de Student à `$\nu$` degrés de liberté peut s'obtenir comme la loi marginale d'un couple de variables aléatoires `$(X,Y)$` telles que :
	$$X \ | \ Y=y \sim \mathcal{N}(0, \nu/y), \quad Y \sim \chi^2_\nu $$

--
.pull-left[

``` r
v <- 15; n <- 10000
Y <- rchisq(n,v)
XcondY <- rnorm(n,0,v/Y)
XY <- cbind(XcondY,Y) 
plot(XY,pch=19) 
```
]
.pull-right[
<img src="cours_files/figure-html/unnamed-chunk-17-1.png" width="432" />
]

---
# Méthode par transformation

On peut également utiliser les relations entre lois jointes, lois conditionnelles et lois marginales.

- **Exemple 2 :**
On appelle *loi mélange* (à densité) une loi qui s'écrit comme une combinaison convexe de densités : 
	`$$f(x) = \sum_{k=1}^K \alpha_k f_k(x), \quad \text{avec } \ \sum_{k=1}^K \alpha_k = 1$$`
La loi de mélange s'obtient comme la loi marginale d'un couple de variables aléatoires `$(X,Z)$` telles que :
`$$X \ | \ Z = k \sim f_k, \quad \mathbb{P}(Z=k) = \alpha_k,$$`

---
# Méthode par transformation

- **Exemple 2 (suite)** : générer un échantillon selon la loi `$f(x) =  0.4 \mathcal{N}(0,1) + 0.6 \mathcal{N}(2,0.5)$`

``` r
Z <- rbinom(n,1,0.6) # on génère une v.a. donnant le numéro de groupe
Xcond <- ifelse(Z==0,rnorm(n),rnorm(n,2,sqrt(0.5))) # loi de X sachant Z
XZ <- cbind(Xcond,Z) # on colle ensemble X et Z pour avoir la loi jointe
hist(XZ[,1],freq=F,breaks=50,main="") 
lines(seq(-4,4,0.01),0.4*dnorm(seq(-4,4,0.01)),col="red",lwd=2)
lines(seq(-4,4,0.01),0.6*dnorm(seq(-4,4,0.01),2,sqrt(0.5)),col="blue",lwd=2)
```

<span style="color:#16A085">**fin du cours 2**</span>

---
name: c3
# Méthode par transformation : Box-Muller

Méthode due à George Box et Mervin Muller (1958), permettant de générer deux échantillons indépendants selon la loi `$\mathcal{N}(0,1)$` uniquement à l'aide de tirages selon la loi `$\mathcal{U}(]0,1])$`.

Elle est utilisée dans certaines bibliothèques de calculs (notamment en C++).

Algorithme de Box-Muller :
> 1. Générer `$U_1$` et `$U_2$` indépendamment selon une loi `$\mathcal{U}(]0,1])$`
2. Poser 
`$$\begin{cases}
X_1 & = \sqrt{-2 \log U_1} \cos (2 \pi U_2) \\
X_2  & = \sqrt{-2 \log U_1} \sin (2 \pi U_2) \\
\end{cases}$$`
3. Renvoyer `$(X_1,X_2)$`, couple de v.a. indépendantes de loi `$\mathcal{N}(0,1)$`

---
# Méthode par transformation : Box-Muller

*Démonstration :*

Soient `$U_1$` et `$U_2$` deux v.a. indépendantes de loi `$\mathcal{U}(]0,1])$`.

1. Quelle est la loi de `$R= \sqrt{-2 \log U_1}$` ? et de `$\theta = 2\pi U_2$` ?

2. Exprimer `$X_1$` et `$X_2$` en fonction de `$R$` et `$\theta$`.

3. Quelle est la loi jointe de `$(X_1,X_2)$` ?

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

</br>

> Implémenter l'algorithme de Box-Muller

</br>
**(Exercice 2 (fin) de la fiche 1)**

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

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

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

---
# 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="cours_files/figure-html/unnamed-chunk-24-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="cours_files/figure-html/unnamed-chunk-25-1.png" width="576" style="display: block; margin: auto;" />

---
# Un cas particulier

`$f(x) = \frac{1}{8} |x^4 - 5x^2 + 4| \ \mathbb{1}_{[-2,2]}(x)$`

.pull-left[

``` r
Y <- runif(10000,-2,2)
U <- runif(10000,0,0.5)

accept <- (U<fdens(Y))
X <- Y[accept]

hist(X,breaks=50,freq=F)
points(x,fx,type="l",col="red")
```

``` r
mean(accept)
```

```
## [1] 0.4966
```

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

.pull-right[
![](cours_files/figure-html/unnamed-chunk-29-1.png)
]

---
# Un cas particulier

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

</br>

---
# Retour sur l'acceptation-rejet

Rappel de l'algorithme de simulation :
> Soit `$X$` une v.a. de densité `$f$` et soit `$g$` une densité de probabilité et une constante `$M \geq 1$` t.q. `$\forall x, f(x) \leq M g(x)$`. Alors pour simuler selon la loi `$f$` il suffit de :
1. simuler `$Y \sim g$`
2. simuler `$U | Y=y \sim \mathcal{U}([0,Mg(y)])$`
3. si `$0 < U < f(Y)$`, poser `$X=Y$`, sinon reprendre l'étape 1.

Autre formulation (plus générale) :
> Soit `$X$` une v.a. de densité `$f$` connue à une constante près (i.e. on connaît `$\tilde{f} = c f$`),  `$g$` une densité de probabilité et `$\tilde{M}$` t.q. `$\forall x, \tilde{f}(x) \leq \tilde{M} g(x)$`. Alors pour simuler selon la loi `$f$` il suffit de :
1. simuler `$Y \sim g$`
2. simuler `$U \sim \mathcal{U}([0,1])$`
3. si `$0 < U < \tilde{f}(Y)/\tilde{M}g(Y)$`, poser `$X=Y$`, sinon reprendre l'étape 1.

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

Soit `$f$` la densité de la loi normale centrée réduite, et `$g$` la densité de la loi de Cauchy.

1.  > Trouver une constante `$M$` telle que `$f(x) \leq M g(x)$` pour tout `$x \in \mathbb{R}$`.
	 
2. >  Proposer un algorithme d'acceptation-rejet pour simuler une loi normale à partir d'une loi de Cauchy, puis un algorithme permettant de générer une loi normale à partir d'une loi uniforme sur `$[0,1]$`.

3. > En moyenne, combien faut-il générer de variables aléatoires de loi uniforme pour obtenir une réalisation selon la loi normale ?

**(Exercice 3 de la fiche 1)**

On cherche à simuler un échantillon i.i.d. de loi `$\mathrm{Beta}(a,b)$`, avec `$a > 1$` et `$b > 1$`.
1. > Pourquoi a-t-on besoin des conditions `$a > 1$` et `$b > 1$` ? Calculer `$M_{a,b}$`, le maximum de `$f_{a,b}$`. En déduire une densité `$g$` telle que `$f_{a,b}(x) \leq M_{a,b} g(x)$`.

3. > Proposer et implémenter un algorithme de type acceptation-rejet pour simuler une variable aléatoire de loi Beta de paramètres `$a$` et `$b$`.

**(Exercice 4 de la fiche 1)**

<span style="color:#16A085">**fin du cours 3**</span>

---
name: c4
# Programme

- Finir la fiche de TD 1
- Résumer le chapitre sur la simulation de v.a.
- Commencer le chapitre sur les méthodes de Monte-Carlo

---
# Résumé du chapitre

- On a vu différentes méthodes pour simuler des v.a. selon une loi cible `$f$`

- méthode de la transformation inverse (en utilisant l'inverse de la f.d.r.)
  - méthode par transformation
  - méthode d'acceptation-rejet

- Pourquoi a-t-on besoin de savoir simuler des v.a. selon une certaine loi `$f$` ? Quelles sont les applications possibles ?

---

class: inverse, middle, center

# Les 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'aiguile 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**

---
# Monte Carlo classique

**Objectif** : calculer une intégrale de la forme `$\mu=\int h(x) f(x) dx$` où `$f$` est une densité de probabilité.

> *Définition.* Soit `$X_1,\dots,X_n$` un échantillon i.i.d. de loi `$f$`. L'estimateur de Monte Carlo de `$\mu$` est :
`$$\hat{\mu}_n = \frac{1}{n} \sum_{i=1}^n h(X_i)$$`

Propriétés :
- estimateur sans biais
- fortement consistant (par la loi forte des grands nombres)
- construction d'intervalles de confiance par 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'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, ...)

- ...

---
# Monte Carlo classique

Algorithme d'estimation par Monte-Carlo:
> 1. Générer `$X_1,...,X_n$` iid selon la loi de densité `$f$`
2. Poser `$\hat{\mu}_n = \frac{1}{n}\sum_{i=1}^n h(X_i)$`

- variance `$\text{Var} (\hat{\mu}_n) = \frac{1}{n} \left(\int h^2(x)f(x)dx - \mu^2\right) = \frac{\sigma_f^2}{n}$` (finie si `$\sigma_f^2$` est finie)
 
--

- On peut estimer `$\sigma_f^2$` par : `$$v_n = \frac{1}{n} \sum_{i=1}^n h^2(X_i) - \hat{\mu}_n^2$$`
 
--

- IC de niveau `$1-\alpha$` en utilisant le TCL appliqué à la suite de v.a. `$(Y_i)$` où `$Y_i=h(X_i)$`. Les `$(Y_i)$` sont bien i.i.d., et si `$\sigma_f^2 < +\infty$` on a :
 `$$\hat{I}_\mu =\left[\hat{\mu}_n - q_{1-\alpha/2}^{\mathcal{N}(0,1)} \sqrt{v_n/n} ; \hat{\mu}_n + q_{1-\alpha/2}^{\mathcal{N}(0,1)} \sqrt{v_n/n} \right]$$`

&rarr; pour diminuer la taille de l'IC par 10, il faut multiplier `$n$` par 100.

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

Approcher les intégrales suivantes par Monte Carlo :

1. > `$$I = \int_0^1(\cos(20x) + \sin(30x))^2 dx$$`

2. > `$$J = \int_0^1 \sin(x^2) e^{-x} dx$$`
	
3. > $$K = \int_0^{+\infty} \sin(x^3) e^{-x} e^{-\frac{x^2}{2}} dx $$

Pour 2. et 3., utiliser deux algorithmes de Monte Carlo différents.

<span style="color:#16A085">**fin du cours 4**</span>

---
name: c5
# Choix de la loi de simulation
**Exemple** : on veut calculer `$\int_0^1 e^{-x}x^{-a} dx$` pour `$0 < a < 1$`.

Monte Carlo "naïf" &rarr; `$f$` densité de la loi uniforme sur `$[0,1]$` et `$h(x) = e^{-x}x^{-a}$`.

``` r
set.seed(09102024)
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.45429531146729, variance : 0.0152093356529852"
```

Peut-on faire mieux ?

---
# Choix de la loi de simulation

A quoi ressemble la fonction `$h(x) = e^{-x}x^{-a}$` ?

---
# Choix de la loi de simulation

A quoi ressemble la fonction `$h(x) = e^{-x}x^{-a}$` ?

<img src="cours_files/figure-html/unnamed-chunk-33-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
# Choix de la loi de simulation

Autre choix : `$g(x) = (1-a) \ x^{-a} \mathbb{1}_{[0,1]}(x)$` et `$\tilde{h}(x) = e^{-x}/ (1-a).$`

<img src="cours_files/figure-html/unnamed-chunk-34-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.46719241157604, variance : 0.000158942591402943"
```

On a diminué la variance par 100 pour une même taille d'échantillon.

--
Vraie valeur : 1.49365 (intégration numérique)

---
# Choix de la loi de simulation

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

**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} = \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** : 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="cours_files/figure-html/unnamed-chunk-39-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 :

$$
`\begin{eqnarray}
\mathbb{P}(X>10) & = & \mathbb{E}(\mathbf{1}_{X>10}) \quad \text{où } X \sim \mathcal{E}(1)\\
                 & = & \int \mathbf{1}_{x>10} f(x)dx \quad \text{où } \ f(x) = e^{-x}\mathbf{1}_{x>0}\\
                 & = & \int \mathbf{1}_{x>10} \frac{f(x)}{g(x)} g(x)dx \quad \text{où} \ g(x) = \frac{1}{10}e^{-x/10}\mathbf{1}_{x>0} \\
                 & = & \int \mathbf{1}_{x>10} \frac{10 \ e^{-x}}{e^{-x/10}} g(x)dx \quad \text{où} \ g(x) = \frac{1}{10}e^{-x/10}\mathbf{1}_{x>0} \\
                 & = & \mathbb{E}\left(\mathbf{1}_{X>10} \ 10 \ e^{-9X/10}\right) \quad \text{où } X \sim \mathcal{E}(1/10)
\end{eqnarray}`
$$
En utilisant Monte-Carlo, on approche cette espérance à l'aide de la moyenne empirique d'un échantillon `$Z_1,\dots,Z_n$` de taille `$n$` de v.a. iid de loi exponentielle de paramètre 1/10.

`$$\hat{p} = \frac{1}{n} \sum_{i=1}^n 10 \ e^{-0.9 Z_i} \mathbf{1}_{Z_i > 10}$$`

---
# Echantillonnage préférentiel

Code correspondant (pour un tirage, donc une estimation) :

``` r
Z <- rexp(n,1/10)
w <- 10*exp(-9*Z/10)  # ou plus généralement : dexp(Z,1)/dexp(Z,1/10)
mean(w*(Z>10))
```

```
## [1] 4.072288e-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="cours_files/figure-html/unnamed-chunk-46-1.png" width="864" />
Avec la loi `$\mathcal{E}(1/10)$` :

<span style="color:#16A085">**fin du cours 5**</span>

---
name: c6
# Echantillonnage préférentiel

**Rappel :**

On cherche à estimer `$\mu = \int h(x)f(x)dx$` par :
`$$\tilde{\mu} = \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),$$`
avec `$Z_1,\dots,Z_n$` un échantillon i.i.d. de loi `$g$`.

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*

---
# Choix de la loi instrumentale

- Calcul de la variance de `$\tilde{\mu}$`

- Variance finie si le ratio `$f/g$` est borné (i.e. queues de distribution de `$g$` plus lourdes)

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

- 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

**Exemple** : approcher `$\int x^2 e^{-x^2/2} dx = \sqrt{2\pi} \ \mathbb{E}(X^2)$` avec `$X \sim \mathcal{N}(0,1)$` (vraie valeur : `$\sqrt{2 \pi} =$` 2.5066283)
 
- avec la loi instrumentale `$\mathcal{N}(0,0.5^2)$`. 
 
--

``` r
n <- 10000
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 : ",var(Z1^2*w1)/n))
```

```
## [1] "Estimation : 1.95613051491758, variance : 0.00728296792972418"
```
--

- avec la loi instrumentale `$\mathcal{N}(0,1.2^2)$` :

``` r
n <- 10000
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 : ",var(Z2^2*w2)/n))
```

```
## [1] "Estimation : 2.5106555058088, variance : 8.54716670546432e-05"
```

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

`$$\bar{\mu} = \frac{\sum_{i=1}^n w_i h(Z_i)}{\sum_{i=1}^n w_i}$$`
--
 
<img src="cours_files/figure-html/unnamed-chunk-54-1.png" width="648" style="display: block; margin: auto;" />

- estimateur **biaisé**
 - fortement consistant
 - variance plus faible dans certains cas

---
# Version auto-normalisée
  
  - Peut s'utiliser lorsque `$f/g$` est connue à une constante multiplicative près
 
--
 
 - Permet également de *simuler* selon la loi `$f$` : il suffit d'échantillonner selon une loi discrète à support sur l'ensemble des `$\{Z_1,\dots,Z_n\}$`, 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))
```

---
# Résumé du chapitre

Dans ce chapitre, on a vu comment approcher des intégrales par des méthodes stochastiques :

- l'approche par Monte Carlo "simple"
- l'échantillonnage préférentiel

Ces méthodes s'utilisent pour approcher une intégrale qui peut s'écrire comme une espérance :
  - soit parce qu'elle est déjà sous la forme `$\int h(x) f(x) dx$` avec `$f$` densité de probabilité
  - soit en faisant apparaître une densité `$f$`

Le choix de la loi d'échantillonnage a un impact sur la **variance** de l'estimation

<span style="color:#16A085">**fin du cours 6**</span>

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

# Modèles de mélange

<!-- ---
# Introduction

Qu'est-ce qu'un .green[**modèle de mélange**] ?

- c'est un modèle statistique, i.e. un modèle que l'on peut appliquer à un échantillon d'observations

- qui suppose que les données observées suivent une *loi mélange*, c'est-à-dire :
 - l'échantillon est constitué de `$K$` sous-populations
 - chaque sous-population suit sa propre distribution
 
 
Les modèles de mélange permettent :
 - de faire de la **classification**
 - de modéliser des jeux de données plus complexes qu'avec une approche "classique" (distribution multi-modale, asymétrique, ...)
-->

---
# Introduction -- exemple 1

- Longueur des ailes de 381 passereaux
<img src="cours_files/figure-html/unnamed-chunk-57-1.png" width="432" style="display: block; margin: auto;" />

- On distingue deux sous-populations

- En fait mâles et femelles sont mélangés, ils n'ont pas été distingués lors des mesures

- Si on savait quels étaient les individus mâles et les individus femelles, on pourrait faire l'analyse séparément par sexe.

---
# Introduction -- exemple 2

- Taux de chlorure dans le sang de 543 patients.

- L'ajustement à une loi normale semble correct sur l'histogramme, mais le QQ-plot indique la présence d'outliers mal pris en compte par un modèle gaussien

- Population constituée d'individus sains et d'individus malades, en proportions inconnues

---
# Loi mélange

On dit qu'un échantillon i.i.d. `$(X_1,\dots,X_n)$` suit un **modèle de mélange** si :
 1. la population est structurée en `$K$` sous-groupes
 2. au sein du sous-groupe `$C_k$`, `$k=1,\dots,K$`, les individus suivent une loi paramétrique `$f_k$` dont les paramètres `$\theta_k$` dépendent du groupe
 
<!--

-->

</br>
--
Pour modéliser une loi mélange :
 - on introduit `$Z_i$`, `$i=1,\dots,n$` telle que `$Z_i = k$` ssi l'individu `$i$` appartient au groupe `$k$`
 - sachant que l'individu `$i$` est dans le groupe `$k$`, on peut définir sa loi `$f_k(\cdot \mid \theta_k)$`

</br>
--
Autrement dit, on a le modèle suivant :

$$ 
`\begin{cases}
\mathbb{P}(Z_i = k) & = & p_k\\
X_i \mid Z_i = k & \sim & f_k(\cdot \mid \theta_k)
\end{cases}`
$$

---
# Loi mélange

> La loi *marginale* de `$X_i$` est 
 $$X_i \sim p_1f_1(\cdot \mid \theta_1)  + \dots + p_K f_K(\cdot \mid \theta_K) $$

On note `$p = (p_1,\dots,p_K)$`, `$\theta=(\theta_1,\dots,\theta_K)$` et `$\phi = (p,\theta)$`.

- Si les lois `$f_k$` sont discrètes, la loi de `$X_i$` l'est aussi
- Si les lois `$f_k$` sont à densité, la loi de `$X_i$` l'est aussi

- Souvent, les lois `$f_1,\dots,f_K$` appartiennent à la même famille de lois, on note alors `$f$` pour simplifier.

- Le nombre de composantes est appelé **ordre du mélange**
- La variable `$Z$` est appelé **variable latente** ou **variable cachée**

---
# Loi mélange

Dans un modèle de mélange, on ne connaît pas les groupes :
 - on ne sait pas comment sont réparties les observations dans les groupes
 - parfois, on ne sait même pas combien il y a de groupes

L'objectif d'un tel modèle peut-être :
 - de prendre en compte une structuration en sous-populations dans l'analyse des données
 - d'identifier les sous-populations (&rarr; clustering)

---
# Lois mélange : une grande diversité

(*source : T. Rebafka*)

---
# <img src="laptop-code-solid.svg" width="45" height="45" align="center"> Exercices

1. > Illustration d'un mélange de:
 - lois Gaussiennes
 - lois de Poisson
 - lois exponentielles

2. > En notant `$\mu_k$` (resp. `$\sigma_k^2$`) l'espérance (resp. la variance) de la loi de la `$k$`-ème composante (i.e. de la loi `$f_k$`), calculer l'espérance et la variance de la loi mélange.

**(Exercices 1 et 2 de la fiche 3)**

---
# Loi mélange
  
Revenons à nos observations :
  
<img src="cours_files/figure-html/unnamed-chunk-61-1.png" width="648" style="display: block; margin: auto;" />

**.center[.orange[Peut-on identifier les différentes composantes du mélange ? ]]**
  
---
# Problématiques
  
Plusieurs questions à aborder :
  
1. Identifiabilité du modèle de mélange

2. Estimation des paramètres

3. Choix du nombre de classes

---
# Identifiabilité
  
Qu'est-ce que l'identifiabilité ?
  
--
  
> **Définition :**
> Une famille de distributions, paramétrée par `$\theta \in \Theta$`, et de densité `$f(\cdot ; \theta)$` est dite .orange[**identifiable**] ssi :
  
> `$$f(x; \theta) = f(x; \theta') \ \text{presque sûrement} \Rightarrow \theta = \theta'$$`

</br>
Autrement dit, s'il existe deux jeux de paramètres `$\theta$` et `$\theta'$` tels que les densités `$f(x; \theta)$` et `$f(x; \theta')$` sont identiques presque sûrement, alors la famille *n'est pas identifiable*.

&rarr; .red[**Est-ce grave ?**]

- si deux jeux de paramètres différents donnent la même loi, alors il est **impossible** d'estimer sans ambiguité les paramètres du modèle.
- cela augmente l'incertitude sur l'estimation, voire cela rend l'estimation caduque.

---
# Identifiabilité

Quelques exemples :
- de lois identifiables :
  - lois Gaussiennes,
  - lois exponentielles,
  - lois de Poisson,
  - la plupart des lois usuelles !

- de non identifiabilité :
  - `$X = \lambda Y$`, avec `$Y \sim \mathcal{N}(0,1)$` et `$\lambda \in \mathbb{R}$`
  - les modèles de mélange !

---
# Identifiabilité

</br>
Il existe différentes situations/raisons pouvant mener à la non identifiabilité d'un modèle de mélange :

1. la loi mélange est invariante par permutation des groupes (.orange[**problème de ré-étiquetage** ou *label switching*])

2. la loi mélange est inchangée si un groupe est divisé en sous-groupes de même distribution (.orange[**problème de surdimensionnement**])

3. la loi mélange peut hériter des problèmes d'identifiabilité des lois la composant, ou de leur mélange

---
# Invariance par permutation des groupes

- La numérotation des groupes est **arbitraire**
- On peut donc permuter les étiquettes de groupes, et obtenir la même loi mélange
- Exemple : les lois mélange `$0.2 \ \mathcal{N}(0,1) + 0.8\ \mathcal{N}(1,2)$` et `$0.8\ \mathcal{N}(1,2) + 0.2 \ \mathcal{N}(0,1)$` sont identiques.

- En fait dans l'écriture :
$$ 
`\begin{cases}
\mathbb{P}(Z_i = k) & = & p_k\\
X_i \mid Z_i = k & \sim & f_k(\cdot \mid \theta_k)
\end{cases}`
$$
le choix du numéro de groupe n'a pas d'importance.

- La seule chose qui importe est d'associer le "bon" `$p_k$` avec le "bon" `$\theta_k$`.

- **La loi mélange de paramètres `$(p_1,\dots,p_K,\theta_1,\dots,\theta_K)$` est identique à la loi mélange `$(p_{\sigma(1)},\dots,p_{\sigma(K)},\theta_{\sigma(1)},\dots,\theta_{\sigma(K)})$`, où `$\sigma(\cdot)$` est une permutation de `$\{1,\dots,K\}$`.**

*(Pour vous en convaincre, vous pouvez utiliser le code de l'exercice 1 de la fiche 3 en permutant les groupes)*

---
# Invariance par permutation des groupes

Que faire ?

&rarr; on peut ajouter **des contraintes** sur les paramètres :

- imposer par exemple `$\theta_1 \leq \theta_2 \leq \dots \leq \theta_K$`

- .green[**Avantage**] : si `$\theta \in \mathbf{R}$`, cela garantit l'identifiabilité
  - .red[**Inconvénients**] : 
      - si `$\theta_k \in \mathbf{R}^q$`, comment définir un ordre ?
      - l'ajout de ces contraintes complique la tâche d'estimation (optimisation sous contraintes)
  
--

- en pratique, on se contente de la notion d'*identifiabilité à permutation des paramètres près*

<span style="color:#16A085">**fin du cours 7**</span>

---
name: c8
# Sur-dimensionnement

- Le modèle de mélange tel qu'il a été défini n'est pas identifiable :
  - il suffit d'ajouter une `$K+1$`-ème composante de poids nul 
  - ou de "couper" une des composantes en deux nouvelles composantes, avec les mêmes paramètres `$\theta$`
--

- Ex. : la loi mélange  `$0.3 \mathcal{N}(0,1) + 0.7 \mathcal{N}(1,2)$`
 est identique à la loi mélange
 `$$0.3 \mathcal{N}(0,1) + 0.7 \mathcal{N}(1,2) + 0 \ \mathcal{N}(5,10)$$` 
 et à la loi mélange
 `$$0.3 \mathcal{N}(0,1) + 0.4 \mathcal{N}(1,2) + 0.3 \mathcal{N}(1,2)$$`

- On aimerait obtenir `$K$` populations **distinctes**.

&rarr; On peut alors ajouter les contraintes :

`$$p_k > 0, \forall \ k, \quad \theta_k \neq \theta_{k'},  \forall \ k \neq k'$$`

---
# Identifiabilité : en résumé

1. Sous les contraintes :
`$$p_k > 0, \forall \ k, \quad \theta_k \neq \theta_{k'},  \forall \ k \neq k'$$`
 &rarr; on peut obtenir l'identifiabilité *à permutation des paramètres près*
  - pour la plupart des modèles de mélange 
  - **condition suffisante** : la famille de lois composant le mélange est .orange[linéairement indépendante]
  - c'est le cas de la plupart des lois usuelles (Gaussiennes, Exponentielles, Gamma, Poisson, ...)
  
2. Sous les contraintes supplémentaires (lorsqu'elles sont possibles) :
  `$$\theta_1 < \dots < \theta_K$$`
  &rarr; on peut obtenir l'identifiabilité "forte"
 
.red[**Même sous ces deux types de contraintes, il existe des lois mélange non identifiables, par ex. les mélanges de Bernoulli ou certains mélanges de lois uniformes**]

---
# <img src="pen-to-square-regular.svg" width="45" height="45" align="center"> Exercices

1. > Montrer qu'un mélange de lois de Bernoulli n'est pas identifiable.
2. > Quelle est la densité de probabilité des deux lois mélanges suivantes :
	`$$\frac{1}{2} \mathcal{U}([-2,1]) + \frac{1}{2} \mathcal{U}([-1,2]),$$`
	`$$\frac{1}{3} \mathcal{U}([-1,1]) + \frac{2}{3} \mathcal{U}([-2,2]).$$`
	Qu'en conclut-on ?
3. > Montrer qu'une loi uniforme sur un segment `$[a,b]$` peut toujours s'écrire comme un mélange de lois uniformes à supports disjoints.

**(Exercice 3 de la fiche 3)**

---
# Estimation

Soit `$X_1, \dots, X_n$`  i.i.d. de loi marginale :
$$ f(x;\theta) = \sum_{k=1}^K p_k f_k(x;\theta_k).$$
**Objectif** : estimer `$\theta = (\theta_1,\dots,\theta_K)$`, les paramètres des lois du mélange, et également `$p=(p_1,\dots,p_K)^t$`, les proportions du mélange.

On note `$\mathbf{X} = (X_1,\dots,X_n)$`. La log-vraisemblance est :

`$$\begin{align}
\ell(\phi;\mathbf{X}) & = \sum_{i=1}^n \ln f(X_i;\phi) \\
     									& = \sum_{i=1}^n \ln \left( \sum_{k=1}^K p_k f_k(X_i;\phi) \right).
\end{align}$$`

**Problème** : il n'existe pas de solution analytique au problème de maximisation de `$\ell(\phi;\mathbf{X})$` car les `$p_k$` ne sont pas estimables facilement

---
# Estimation

- Que se passerait-il si on observait aussi `$(Z_1,\dots,Z_n)$` ?
- Notre échantillon serait alors constitué des couples `$(X_i,Z_i), i=1,\dots,n$`
- On parle de **données complètes**
- La loi jointe de `$(X_i,Z_i)$` est :

`$$f(x,z;\phi) = \prod_{k=1}^K \left(p_k f_k(x;\phi)\right)^{\mathbf{1}_{z=k}}$$`

La log-vraisemblance serait alors :
`$$\begin{align}
\ \ell(\phi;\mathbf{X}, \mathbf{Z})&  =\sum_{i=1}^n \ln f(X_i,Z_i;\phi) \\
& = \sum_{i=1}^n  \sum_{k=1}^K \mathbb{1}_{Z_{i}=k} \ln \left(p_k f_k(X_i;\phi) \right)
\end{align}$$`

---
# Estimation

- Les deux quantités précédentes vont jouer un rôle particulier dans l'estimation
- Notations :
  - la log-vraisemblance des **données observées** est notée `$\ell_{\text{obs}}$` :
`$$\ell_{\text{obs}}(\phi;\mathbf{X}) = \sum_{i=1}^n \ln \left( \sum_{k=1}^K p_k f_k(X_i;\phi) \right)$$`
  - la log-vraisemblance des **données complètes** est notée `$\ell_{\text{comp}}$` :

`$$\ell_{\text{comp}}(\phi;\mathbf{X},\mathbf{Z}) = \sum_{i=1}^n  \sum_{k=1}^K \mathbf{1}_{Z_i=k} \ln \left(p_k f_k(X_i;\phi) \right)$$`
- La log-vraisemblance complète est plus facile à manipuler car on a "fait sortir" la somme du log.

- On va exploiter i) le fait que `$\ell_{\text{comp}}$` est plus facile à maximiser et ii) le lien entre ces deux fonctions, pour maximiser `$\ell_{\text{obs}}$`.

.center[.orange[&rarr; **c'est le principe de l'algorithme EM**]]

---
# Algorithme EM

- Permet d'obtenir un maximum .red[**LOCAL**] de la log-vraisemblance en présence de variables latentes
	
- Il a été introduit par Dempster, Laird et Rubin en 1977.

- Particulièrement adapté aux cas où la vraisemblance complète s'écrit plus simplement que la vraisemblance observée

- Repose sur l'idée suivante : en présence de données manquantes, une première intuition est **d'estimer ou de remplacer** ces données manquantes, puis d'estimer les paramètres du modèle à l'aide des **données "augmentées"**

---
# Algorithme EM

On a : `$\quad \quad \quad \quad \quad \quad \quad \quad \quad f(x,z;\phi) = f(x \mid z;\phi) f(z;\phi)$`

mais aussi : `$\quad \quad \ \ \ \quad \quad \quad \quad f(x,z;\phi) = f(x;\phi) f(z \mid x;\phi)$`

et donc `$\quad \quad \quad \ \quad \quad \quad \quad \ln f(x,z;\phi) = \ln f(x;\phi) + \ln f(z \mid x;\phi)$`

D'où :

`$$\newcommand{\green}[1]{\color{green}{#1}}
\newcommand{\orange}[1]{\color{orange}{#1}}
\newcommand{\red}[1]{\color{red}{#1}}
\newcommand{\blue}[1]{\color{blue}{#1}}$$`

`$$\begin{align}
\sum_{i=1}^n \ln f(X_i,Z_i;\phi) & = \sum_{i=1}^n \ln f(X_i;\phi) + \sum_{i=1}^n \ln f(Z_i \mid X_i;\phi) \\
\underbrace{\orange{\sum_{i=1}^n \ln f(X_i,Z_i;\phi)}}_{\text{vrais.} \orange{\text{ complète}}} & = \underbrace{\green{\sum_{i=1}^n \ln f(X_i;\phi)}}_{\text{vrais.} \green{\text{ observée}}} + \sum_{i=1}^n \ln f(Z_i \mid X_i;\phi) \\
\Rightarrow \ell_{\text{obs}} (\phi;\mathbf{X}) & = \ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) - \ln f(\mathbf{Z} \mid \mathbf{X}; \phi)
\end{align}$$`

---
# Algorithme EM

- Mais `$\mathbf{Z}$` est inconnu ...

- L'idée est .orange[**d'intégrer**] l'expression précédente par rapport à la loi conditionnelle de `$\mathbf{Z}$` sachant `$\mathbf{X}$`, afin d'obtenir une expression qui ne dépende plus de `$\mathbf{Z}$` (c'est-à-dire, calculer l'.red[**espérance conditionnelle**] sachant `$\mathbf{X}$`.)

- On se place à l'itération `$m+1$` de l'algorithme, avec l'estimation courante `$\phi^{(m)}$`.

- On va intégrer l'expression : `$$\ell_{\text{obs}} (\phi;\mathbf{X}) = \ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) - \ln f(\mathbf{Z} \mid \mathbf{X}; \phi)$$`
par rapport à la loi `$f(z \mid x; \phi^{(m)})$`, i.e. par rapport à la loi conditionnelle de `$\mathbf{Z}$` sachant `$\mathbf{X}$` et la valeur actuelle du paramètre `$\phi^{(m)}$`.

- On note `$\mathbb{E}_{\phi^{(m)}}\left[ \cdot \mid \mathbf{X} \right]$` l'espérance par rapport à cette loi. Autrement dit, si on considère une fonction `$h(\mathbf{X},\mathbf{Z})$`, on a :
`$$\mathbb{E}_{\phi^{(m)}}[h(\mathbf{X},\mathbf{Z})\mid \mathbf{X} = x]= \int h(x,z) f(z \mid x; \phi^{(m)}) dz$$`

---
# Algorithme EM

- Comme `$\ell_{\text{obs}} (\phi;\mathbf{X})$` est `$\sigma(\mathbf{X})$`-mesurable, on a : `$\mathbb{E}_{\phi^{(m)}}[\ell_{\text{obs}} (\phi;\mathbf{X})\mid \mathbf{X}] = \ell_{\text{obs}} (\phi;\mathbf{X}).$` En effet :

`$$\begin{align}
\mathbb{E}_{\phi^{(m)}}[\ell_{\text{obs}} (\phi;\mathbf{X})\mid \mathbf{X} = x] & = \int \ell_{\text{obs}} (\phi;x) f(z \mid x; \phi^{(m)}) dz \\
 & = \ell_{\text{obs}} (\phi;x) \int f(z \mid x; \phi^{(m)}) dz \\
 & = \ell_{\text{obs}} (\phi;x)
\end{align}$$`

- D'où :
`$$\begin{align}
\ell_{\text{obs}} (\phi;\mathbf{X}) &  = 	\orange{\mathbb{E}_{\phi^{(m)}}\left[\ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) \mid \mathbf{X} \right]}  \ \blue{-  \ \mathbb{E}_{\phi^{(m)}}\left[ \ln f(\mathbf{Z} \mid \mathbf{X}; \phi) \ \middle\vert \ \mathbf{X}\right]} \\
 & = \orange{Q(\phi ; \phi^{(m)})} \ \red{+} \ \blue{H(\phi ; \phi^{(m)})}.
\end{align}$$`

- Que peut-on dire de `$Q$` et `$H$` ?

1. on montre que `$H(\phi ; \phi^{(m)}) -  H(\phi^{(m)} ; \phi^{(m)}) \geq 0$` .red[**&rarr; au tableau**]
  2. on peut alors se contenter d'étudier `$Q$`
  
---
# Algorithme EM

Idée générale :

- A l'itération `$m+1$` de l'algorithme EM, la valeur de `$H$` augmente `$\rightarrow$` 
il suffit de s'intéresser à la fonction `$Q$`

- Chaque itération de l'algorithme EM se divise en deux étapes :
  - l'étape d'.orange[**espérance**] (étape E) qui consiste à calculer la quantité `$Q(\phi ; \phi^{m})$` en fonction de la valeur courante de l'estimation `$\phi^{m}$`
	- l'étape de .green[**maximisation**] (étape M), qui consiste à maximiser la fonction `$Q$`
	
- L'objectif est de maximiser la log-vraisemblance observée en maximisant `$Q$`.

---
# Algorithme EM

Description de l'algorithme

1. **initialisation** :  valeur initiale `$\phi^{0}$`, nombre total d'itérations `$M$`\;

2. **pour** `$m=1,\dots,M$` :
	- .orange[**Etape E (Espérance)**] :	calcul de `$Q$`
	`$$Q(\phi ; \phi^{(m-1)}) =	\mathbb{E}_{\phi^{(m-1)}}\left(\ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) \mid \mathbf{X} \right)$$`
	- .green[**Etape M (Maximisation)**] : maximisation de `$Q$`
		`$$\phi^{(m)} = \arg \max_\phi Q(\phi ; \phi^{(m-1)})$$`

<!--
---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-62-1.png" width="1000px" style="display: block; margin: auto;" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-63-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-64-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-65-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-66-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-67-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-68-1.png" width="1500px" />

---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-69-1.png" width="1500px" />
 
-->
---
# Illustration

</br>
<img src="cours_files/figure-html/unnamed-chunk-70-1.png" width="1500px" />

<span style="color:#16A085">**fin du cours 8**</span>

---
name: c9
# Algorithme EM

Description de l'algorithme

1. **initialisation** :  valeur initiale `$\phi^{0}$`, nombre total d'itérations `$M$`\;

---
# Algorithme EM pour les modèles de mélange

.large[.orange[**Etape E**]] : calcul de `$Q$`

- Dans le cas du modèle de mélange, les `$Z_i$` sont des variables aléatoires discrètes

- On montre (.red[**au tableau**]) que :
 `$$Q(\phi ; \phi^{(m-1)})	= \sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{(m-1)}) \ln p_k  + \sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{(m-1)}) \ln f_k(X_i;\theta_k)$$`
 avec
`$$t_k(X_i;\phi^{(m-1)})  = \frac{p_k^{(m-1)} f_k(X_i;\theta_k^{(m-1)})}{\sum_{k=1}^K p_k^{(m-1)} f_k(X_i;\theta_k^{(m-1)})}$$`

---
# Algorithme EM pour les modèles de mélange

.large[.green[**Etape M**]] : maximisation de `$Q$`

- La fonction `$Q$` peut se décomposer en deux parties : .blue[une qui dépend de] `$\blue{p}$` et .red[une autre qui dépend de] `$\red{\theta}$`

`$$Q(\phi ; \phi^{(m-1)}) = \blue{\sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{(m-1)}) \ln p_k } + \red{\sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{(m-1)}) \ln f_k(X_i;\theta_k)}$$`

&rarr; on maximise indépendamment pour `$\theta$` et pour `$p$`.

Si `$f_k$` est la densité Gaussienne, on obtient :
`$$\begin{align}
p_k^m & = \frac{1}{n} \sum_{i=1}^n t_k(X_i; \phi^{(m-1)})\\
\mu_k^m  & = \dfrac{\sum_{i=1}^n t_k(X_i; \phi^{m-1})\  X_i}{\sum_{i=1}^n t_k(X_i; \phi^{(m-1)})}\\
\Sigma_k^m  & = \frac{\sum_{i=1}^n t_k(X_i; \phi^{m-1})\  (X_i - \mu_k^m) (X_i - \mu_k^m)^t}{\sum_{i=1}^n t_k(X_i; \phi^{(m-1)})}
\end{align}$$`

---
# Quelques remarques

- On répète les étapes E et M jusqu'à convergence de l'algorithme, i.e. jusqu'à stabilisation de la fonction `$Q$`

- `$K$` doit être fixé en amont  `$\rightarrow$` lancer l'algorithme pour différentes valeurs de `$K$` et comparer AIC, BIC, ...

- Il existe des variantes .orange[**stochastiques**] de l'algorithme, basés sur une approximation de l'étape E par des méthodes de type Monte Carlo

- Chaque itération de l'algorithme EM permet d'augmenter la valeur de la log-vraisemblance
	- sous certaines conditions, l'algorithme .red[converge vers un point stationnaire de la vraisemblance]
	- Ce point stationnaire peut être un maximum local, un maximum global, ou un point selle
	- Dans la pratique on lance plusieurs fois l'algorithme avec des initialisations différentes

---
# Résumé

- **Objectif** : maximiser la log-vraisemblance *observée* `$\ell_{obs}(\phi;\mathbf{X}) = \sum_{i=1}^n \ln f(X_i;\phi)$`

- **Problème** : cette fonction n'est pas maximisable *explicitement*

- **Solution** : utiliser l'algorithme EM :
 - on maximise une *fonction auxiliaire* plus simple à manipuler : (à l'itération `$m$`)
 `$Q(\phi ; \phi^{(m-1)}) =	\mathbb{E}_{\phi^{(m-1)}}\left(\ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) \mid \mathbf{X} \right)\quad$` (à l'itération `$m$`)
 - cette fonction est une borne inférieure de `$\ell_{obs}(\phi;\mathbf{X})$` &rarr; en augmentant `$Q$` on augmente `$\ell_{obs}(\phi;\mathbf{X})$`

- .green[**Avantages**] :
 - sous certaines conditions, permet d'obtenir l'EMV
 - extensions stochastiques pour augmenter son implémentabilité

- .red[**Inconvénients**] :
 - dans le cas général, convergence vers un point stationnaire de la log-vraisemblance (pas forcément le max)
 - nécessite d'expliciter les étapes E et M pour chaque modèle de mélange

---
# Estimation par algorithme EM

Dans le cas Gaussien (conclusions similaires dans le cas général) :

.pull-left[

Modèle de mélange et algorithme EM :
`$$\small\begin{align*}
\hat{p}_k 		 & = \frac{\sum_{i=1}^n \red{t_k(X_i; \phi^{M})}}{n}\\
\hat{\mu}_k  & = \frac{\sum_{i=1}^n \red{t_k(X_i; \phi^{M})} X_i}{\sum_{i=1}^n \red{t_k(X_i; \phi^{M})}}\\
\hat{\Sigma}_k  & = \frac{\sum_{i=1}^n \red{t_k(X_i; \phi^{M})} (X_i - \mu_k)^t (X_i - \mu_k)^t}{\sum_{i=1}^n \red{t_k(X_i; \phi^{M})}}
\end{align*}$$`
]
.pull-right[
Si les `$Z_i$` étaient connues, EMV exact :
`$$\small\begin{align*}
\hat{p}_k 		 & = \frac{\sum_{i=1}^n \blue{\mathbf{1}_{Z_i=k}}}{n}\\
\hat{\mu}_k  & = \frac{\sum_{i=1}^n \blue{\mathbf{1}_{Z_i=k}} X_i}{\sum_{i=1}^n \blue{\mathbf{1}_{Z_i=k}}}\\
\hat{\Sigma}_k  & = \frac{\sum_{i=1}^n \blue{\mathbf{1}_{Z_i=k}} (X_i - \mu_k)^t (X_i - \mu_k)^t}{\sum_{i=1}^n \blue{\mathbf{1}_{Z_i=k}}}
\end{align*}$$`
]

où `$t_k(X_i; \phi^{M}) = \mathbb{P}(Z_i = k \mid X_i; \phi^M)$`

&rarr; si les `$Z_i$` étaient connues, on saurait à quelles classes appartiennent les observations, et on utiliserait, pour estimer les paramètres de la classe `$k$`, les observations de cette classe

&rarr; dans le cas des modèles de mélange, les `$Z_i$` ne sont pas observées : pour le calcul des paramètres associés à la classe `$k$`, on **pondère** alors chaque observation par la probabilité, calculée à l'itération `$m-1$`, pour que l'individu correspondant appartienne au groupe `$k$`, sachant ses observations.

---
# Estimation de la partition

A partir des estimations obtenues en sortie de l'algorithme EM :
- on peut déduire une règle d'affectation des individus dans les différentes classes en utilisant le principe du .blue[Maximum A Posteriori (MAP)] :

`$$i \in C_\ell \Leftrightarrow \ell = \arg \max_{k =1,\dots,K} t_k(X_i ; \hat{\phi})$$`
où `$\hat{\phi}$` est l'estimateur obtenu lorsque l'algorithme a convergé. On affecte donc l'individu `$i$` à la classe dont la probabilité a posteriori est la plus élevée.

- On utilise ensuite les estimations obtenues pour `$\theta_k$` pour interpréter les classes.

---
# Exemple

Données sur les éruptions du geyser 'Old faithful' du parc Yellowstone (disponibles sour `R` dans la libraire `MASS`):

.pull-left[
<img src="cours_files/figure-html/unnamed-chunk-71-1.png" width="432" style="display: block; margin: auto;" />
]

---
# Exemple

Données sur les éruptions du geyser 'Old faithful' du parc Yellowstone (disponibles sour `R` dans la libraire `MASS`):

```
## [1] "Poids de chaque composante : 0.3096, 0.6904"
```

```
## [1] "Moyenne dans chaque groupe : 54.26, 80.41"
```

```
## [1] "Variance dans chaque groupe : 25, 55.69"
```

---
# Classification par l'algorithme CEM

- L'algorithme EM a pour objectif principal d'estimer les paramètres du modèle, à savoir `$p_k$` et `$\theta_k$`

- On ne tient pas compte explicitement de l'objectif de classification

- On peut obtenir en pratique des classes qui ne soient pas suffisamment séparées ce qui rend plus difficile l'interprétation des résultats

- Autre approche aussi basée sur un algorithme de type EM a été proposée par Celeux et Govaert (1992), intégrant de façon explicite cet objectif de classification

`$\rightarrow$` algorithme .orange[**CEM**] (Classification - EM)

---
# Algorithme CEM

On a :
`$$\ell_{\text{obs}} (\phi;\mathbf{X}) = \ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}) - \sum_{i=1}^n \sum_{k=1}^K \mathbf{1}_{Z_i=k} \ln t_k(\mathbf{X}_i; \phi)$$`

- Si on a une partition idéale, avec des classes parfaitement séparées, alors les `$t_k(\mathbf{X}_i; \phi)$` sont binaires (il n'y a pas d'ambiguïté quant au choix de la classe)

- Le deuxième terme ci-dessus est alors .red[**nul**]

- Il s'agit en fait d'un terme d'entropie (fonction `$H$`) : mesure l'écart entre la partition définie par les `$Z_{i}$` et les probabilités a posteriori `$t_k(\mathbf{X}_i; \phi)$`.

- &rarr; la log-vraisemblance observée et la log-vraisemblance complète sont égales en cas de partition idéale

- L'algorithme CEM cherche alors à maximiser .red[**directement la vraisemblance complète**], au lieu de chercher à maximiser la log-vraisemblance observée à l'aide de la fonction `$Q$`.

---
# Algorithme CEM

Algorithme CEM

- **initialisation** : valeur initiale `$\phi^0$`, nombre total d'itérations `$M$`

- **pour `$m=1,\dots,M$`** :

- .orange[**Etape E**] : identique à l'étape E de l'algorithme EM. On obtient
		`$$t_k(X_i;\phi^{m-1}) = \frac{p_k^{m-1} f_k(X_i;\theta_k^{m-1})}{\sum_{k=1}^K p_k^{m-1} f_k(X_i;\theta_k^{m-1})}$$`

- .blue[**Etape C (Classification)**] : création de la partition.
		`$$Z^m_{i} = \arg \max_{\ell = 1, \dots, K} t_\ell(X_i;\phi^{m-1})$$`

- .green[**Etape M**] : maximisation de la vraisemblance complète
		`$$\ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}^m) = \sum_{i=1}^n  \sum_{k=1}^K \mathbf{1}_{Z^m_i=k} \ln  \left(p_k f_k(X_i;\theta_k) \right)$$`

---
# L'étape de maximisation

- L'expression de la vraisemblance complète est très proche de celle de la fonction `$Q$`, avec les `$\mathbf{1}_{Z_{i}^{m-1}=k}$` qui remplacent les `$t_k(X_i;\phi^{m-1})$`:

`$$\begin{align*}
Q(\phi ; \phi^{m-1}) & = \sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{m-1}) \ln p_k  & + & \sum_{i=1}^n \sum_{k=1}^K t_k(X_i;\phi^{m-1}) \ln f_k(X_i;\theta_k) \\
\ell_{\text{comp}} (\phi;\mathbf{X},\mathbf{Z}^m) & = \sum_{i=1}^n  \sum_{k=1}^K \mathbf{1}_{Z^m_i=k} \ln  p_k  & + & \sum_{i=1}^n  \sum_{k=1}^K \mathbf{1}_{Z^m_i=k}  \ln f_k(X_i;\theta_k)
\end{align*}$$`

- On obtient alors, dans le cas d'un mélange Gaussien:
`$$\begin{align*}
p_k^m & = \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{Z^m_i=k} = \frac{n_k^{(m)}}{n}, \quad \text{avec } n_k^{(m)} \text{ le nb d'obs. classées dans le groupe } k\\
\mu_k^m  & = \dfrac{\sum_{i=1}^n \mathbf{1}_{Z^m_i=k}  X_i}{\sum_{i=1}^n \mathbf{1}_{Z^m_i=k}} = \frac{\sum_{i=1}^n \mathbf{1}_{Z^m_i=k}  X_i}{n_k^{(m)}}\\
\Sigma_k^m  & = \frac{\sum_{i=1}^n \mathbf{1}_{Z^m_i=k} (X_i - \mu_k^m) (X_i - \mu_k^m)^t}{\sum_{i=1}^n Z^m_{ik}} = \frac{\sum_{i=1}^n \mathbf{1}_{Z^m_i=k} (X_i - \mu_k^m) (X_i - \mu_k^m)^t}{n_k^{(m)}}.
\end{align*}$$`

---
# Quelques remarques

- En pratique, l'algorithme CEM converge souvent .orange[**plus vite**] que l'algorithme EM

- Plus efficace pour réaliser un compromis biais-variance surtout pour les petits échantillons

- Les estimateurs obtenus avec l'algorithme CEM ne sont .red[**pas sans biais**], même asymptotiquement.

---
# Choix du nombre de classes

- Jusqu'à présent, on a considéré que le nombre de classes `$K$` était connu et fixé

- En pratique ce n'est pas toujours le cas &rarr; `$K$` est également un paramètre du modèle !

- Quelques pistes :
 - identifier le .orange[**nombre de modes**]

- faire un .orange[**test d'hypothèses**] portant sur le paramètre `$K$`

- comparer plusieurs choix de `$K$` à l'aide de critères de .orange[**comparaison de modèles**]

---
# Nombre de modes

- En dimension 1 :

---
# Nombre de modes

- En dimension 1 :

- Certains modes peuvent cacher plusieurs composantes
- Certaines densités de mélange sont **unimodales**
- Plus difficile en dimension > 1

&rarr; en cherchant les modes, on obtient un **minorant** du nombre de groupes

---
# Nombre de modes

Sur l'exemple des données du geyser 'Old faithful', comparaison du modèle obtenu en dimension 1 ou en dimension 2
:

---
# Test d'hypothèses

- Exemple : tester
`$$H_0 : K=K_0 \quad \text{ contre } \quad H_1 : K=K_0+1$$`
à l'aide d'un test du rapport de vraisemblances

- Pour mettre en place ce test il faut calculer la statistique de test :
`$$T_n = -2 \left( \ell_{obs}(\hat{\phi}_0;\mathbb{x}) - \ell_{obs}(\hat{\phi}_1;\mathbb{x}) \right)$$`
--

- Il faut donc calculer l'estimateur du maximum de vraisemblance sous `$H_0$` et sous `$H_1$`

- Mais, si les données proviennent du modèle sous `$H_0$` :
 - pour calculer `$\hat{\phi}_1$`, il faut ajuster un modèle à `$K_0+1$` composantes sur des données dont la loi de mélange contient `$K_0$` composantes
 - on a alors un problème d'identifiabilité
 - les conditions de régularité standards permettant d'obtenir la consistance de l'EMV .red[**ne sont pas vérifiées**]

- De même, les résultats standards sur la loi de la statistique du RV ne s'appliquent pas

---
# Test d'hypothèses

Plusieurs auteurs ont étudiés la consistance des estimateurs dans le cas où on ajuste un modèle avec 1 composante "de trop" :
- soit on converge (quand `$n$` tend vers `$+\infty$`) vers deux composantes ayant les mêmes paramètres
- soit on converge vers une composante avec un poids nul

Pour la loi (asymptotique) de la statistique de test :
- son expression exacte est difficile voire impossible à exploiter
- on peut utiliser une approche *bootstrap*

---
# Test d'hypothèses

Test `$\quad H_0 : K=K_0 \quad \text{ contre } \quad H_1 : K=K_0+1$`

**Procédure bootstrap :**

> 1. Ajuster les deux modèles à `$K_0$` et `$K_0+1$` composantes par maximum de vraisemblance pour obtenir `$\hat{\phi}_0=(\hat{p}_1,\dots,\hat{p}_{K_0},\hat{\theta}_1,\dots,\hat{\theta}_{K_0})$`, `$\hat{\phi}_1$` et `$T_n$`
> 2. Simuler `$B$` jeux de données selon la loi de mélange à `$K_0$` composantes en utilisant `$\hat{\phi}_{0}$` :
  - simuler `$Z_i^{b,*}$`, `$i=1,\dots,n$` selon la loi `$\mathbb{P}(Z_i^{b,*} = k)=\hat{p}_k$`
  - conditionnellement à la valeur de `$Z_i^{b,*}$`, simuler `$X_i^{b,*} \mid Z_i^{b,*} = k \sim f_k(\cdot \mid \hat{\theta}_k)$`
3. Sur chaque jeu de données bootstrap `$(X_1^{b,*},\dots,X_n^{b,*})$`:
    - estimer les deux modèles à `$K_0$` et `$K_0+1$` composantes (i.e. sous `$H_0$` et sous `$H_1$`) par maximum de vraisemblance.
    - calculer la statistique de test bootstrap `$T_n^{b,*}$`
4. Estimer la `$p$`-valeur bootstrap du test par `$p_{boot} = \frac{1}{B}\sum_{i=1}^B \mathbb{1}_{T_n^* \geq T_n}$`

Mise en place dans plusieurs packages `R` : `mclust`, `mixtools`, `mixR`, ...

---
# Méthodes de comparaison de modèles

On lance l'algorithme en faisant varier le nombre de classes, et on compare les résultats à l'aide de différents critères. On note `$q$` la dimension du modèle, c'est-à-dire le nombre total de paramètres.

- AIC : `$-2 \ \green{\ell_{\text{obs}} (\hat{\phi} ; \mathbf{X})} + 2 q$`

- BIC : `$-2 \ \green{\ell_{\text{obs}} (\hat{\phi} ; \mathbf{X})} + q \ln n$`

-	ICL : `$-2 \ \orange{\ell_{\text{comp}} (\hat{\phi} ; \mathbf{X}, \hat{\mathbf{Z}})} + q \ln n$`

- On a ICL = BIC `$\blue{- 2 \sum_{i=1}^n \sum_{k=1}^K  \mathbf{1}_{\hat{Z}_i=k} \ln t_k(X_i;\hat{\phi})}$` (BIC + terme d'entropie)

- L'entropie a tendance à augmenter quand le nombre de classes augmente

- ICL sélectionne des modèles avec **moins de classes** que BIC &rarr; ce critère serait meilleur pour sélectionner le nombre de classes dans une classification.

---
# Exemple

Données du geyser 'Old faithful', comparaison de deux modèles en dimension 1 :

.pull-left[
<img src="cours_files/figure-html/unnamed-chunk-76-1.png" width="360" style="display: block; margin: auto;" />

```
## [1] "BIC =  2343.61089651888"
```

```
## [1] "ICL =  2358.0937464077"
```

```
## [1] "Entropie :  7.24142494441026"
```
]
.pull-right[
<img src="cours_files/figure-html/unnamed-chunk-77-1.png" width="360" style="display: block; margin: auto;" />

```
## [1] "BIC =  2351.60546749235"
```

```
## [1] "ICL =  2458.8121392554"
```

```
## [1] "Entropie :  53.6033358815271"
```
]

---
# Exemple

Incertitudes de classification :

.pull-left[
<img src="cours_files/figure-html/unnamed-chunk-78-1.png" width="360" style="display: block; margin: auto;" />
]

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

]

---
# <img src="pen-to-square-regular.svg" width="45" height="45" align="center"> Entropie d'une classification

Pour mieux comprendre l'entropie, plaçons nous dans le cas du modèle de mélange suivant:

`$$Y \sim p_1\mathcal{P}(\lambda_1) + p_2\mathcal{P}(\lambda_2) + p_3\mathcal{P}(\lambda_3) + p_4\mathcal{P}(\lambda_4)$$`

avec `$p_1=p_2=p_3=1/5,p_4=2/5$` et `$\lambda_1=1,\lambda_2=2,\lambda_3=3,\lambda_4=4$`.

**Exercice**
> 1. Simuler un échantillon de taille `$n=20$` selon la loi de `$Y$` (voir fiche 3).
> 2. Donner l'expression de `$t_k(X_i;\phi)=\mathbf{P}(Z_i = k | X_i ; \phi)$`.
> 2. En supposant les `$(p_k,\lambda_k)$` connus, calculer pour chaque observation, la probabilité qu'elle appartienne à chacune des 4 classes.
> 3. En se basant sur ces probabilités, proposer une classification.
> 4. Calculer l'entropie de la classification, définie par :
`$$H(\phi ; \phi) = - \sum_{i=1}^n \sum_{k=1}^K  t_k(X_i;\phi) \ln t_k(X_i;\phi)$$`

> 5. Comparer à l'entropie obtenue en affectant à la classe `$k$` les observations pour lesquelles `$t_k(X_i;\phi) > 0.5$`, lorsque cela est possible