Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)

Outline

Nouveaux outils informatiques

pour la Statistique exploratoire 1 Simulation de variables aléatoires
(=NOISE)
2 Méthodes de Monte Carlo et algorithme EM

Christian P. Robert 3 Méthode du bootstrap

Université Paris Dauphine 4 Statistique non–paramétrique

[Link] xian

Licence MI2E, 2006–2007

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction

Chapitre 1 : Introduction
Simulation de variables aléatoires
Besoin de “produire le hasard” par ordinateur
Evaluer le comportement d’un système complexe (programme,
Introduction réseau, file d’attente, système de particules, atmosphère,
Générateur pseudo-aléatoire épidémie, actions...)
Distributions non-uniformes (1)
Déterminer les propriétés probabilistes d’une procédure
Distributions non-uniformes (2)
statistique non-standard ou sous une loi inconnue [bootstrap]
Méthodes de Markov
Validation d’un modèle probabiliste
Approcher une espérance/intégrale sous une loi non-standard
[loi des grands nombres]
Maximiser une fonction/vraisemblance faiblement régulière
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

n= 4 n= 8 n= 16

10 15 20 25

20
30

15
20
Example (TCL pour la loi binomiale)

10
10

5
5
0

0
Si 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Xn ∼ B(n, p) , n= 32 n= 64 n= 128

14

10 15 20 25

15
0 2 4 6 8 10

10
Xn converge en loi vers la loi normale :

5
5
0

0
  0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.3 0.4 0.5 0.6 0.35 0.40 0.45 0.50 0.55 0.60 0.65

√ n→∞ p(1 − p)
n (Xn − p) N 0, n= 256 n= 512 n= 1024

30

10 20 30 40 50
5 10 15 20 25
20
5 10
0

0
0.40 0.45 0.50 0.55 0.60 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.46 0.48 0.50 0.52 0.54

Histogrammes de 1000 réalisations B(n, .5)

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

Example (Minimisation aléatoire)

On considère la fonction
65
4

h(x, y) = (x sin(20y) + y sin(20x))2 cosh(sin(10x)x)

Z
3

+ (x cos(10y) − y sin(10x))2 cosh(cos(20y)y) ,

2
1
0

à minimiser. (On sait que le minimum global vaut 0 en 0.5

1

(x, y) = (0, 0).) 0

Y
0.5

0
X
-0.
5

-0.5

-1
-1
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

Example (Minimisation aléatoire (2))

Au lieu de chercher à résoudre les équations du premier ordre

0.8
∂h(x, y) ∂h(x, y)
= 0, =0
∂x ∂y

0.6
et à vérifier les conditions du second ordre, on peut générer la suite
aléatoire dans R2
αj

0.4
θj+1 = θj + ∆h(θj , βj ζj ) ζj
2βj

où

0.2
⋄ les ζj sont uniformes sur le cercle unité x2 + y 2 = 1; -0.2 0.0 0.2 0.4 0.6

⋄ ∆h(θ, ζ) = h(θ + ζ) − h(θ − ζ);

Cas où αj = 1/10 log(1 + j) et βj = 1/j
⋄ (αj ) et (βj ) tendent vers 0

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

Problème du voyageur de commerce Problème NP-complet

Problème du voyageur de
commerce représentatif de
Problème classique d’allocation: problèmes mathématiques
durs à temps de résolution
Représentant devant visiter explosifs
un ensemble de n villes Nombre de chemins possibles
Coûts de voyages entre deux n! et solutions exactes
villes fixés [et différents] disponibles en temps O(2n )
Recherche du coût global Problème à nombreuses
minimum applications (réseaux,
conception de circuits
imprimés, séquençage de
génome, etc.) Concours Procter & Gamble
1962
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

Problème toujours ouvert Résolution par simulation

Algorithme du recuit simulé:

Répéter
Modifications aléatoires de parties du circuit de coût C0
Evaluation du coût C du nouveau circuit
Acceptation du nouveau circuit avec probabilité
 
C0 − C
exp ∧1
T

T , température, est réduite progressivement.

Solution exacte pour 15, 112 Résolution pour les 24, 978 villes
[Metropolis, 1953]
villes allemandes trouvée en 2001 suédoises en 2004 en 84.8 années
en 22.6 années CPU. CPU

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Introduction

Illustration Pricing d’options

Example (400 villes)
Calcul complexe d’espérances/valeurs moyennes d’options, E[CT ],
nécessaire pour évaluer le prix d’achat (1 + r)−T E[CT ]

Example (Options européennes)

Cas où
CT = (ST − K)+
avec

ST = S0 × Y1 × · · · × YT , Pr(Yi = u) = 1 − Pr(Yi = d) = p .

Résolution par simulation des binomiales Yi

T = 1.2
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Introduction Générateur pseudo-aléatoire

Pricing d’options (suite) Générateur pseudo-aléatoire

Elément central des méthodes de simulation : elles reposent toutes
sur la transformation de variables uniformes U (0, 1)
Example (Options asiatiques) Definition (Générateur pseudo-aléatoire)
Modèle en temps continu où Un générateur pseudo-aléatoire est une transformation
 Z + !+ déterministe Ψ de ]0, 1[ dans ]0, 1[ telle que, pour toute valeur
T
1 T
1X initiale u0 et tout n, la suite
CT = S(t)dt − K ≈ S(n) − K ,
T 0 T
n=1
{u0 , Ψ(u0 ), Ψ(Ψ(u0 )), . . . , Ψn (u0 )}
avec
a le même comportement statistique qu’une suite iid U (0, 1)
iid
S(n + 1) = S(n) × exp {∆X(n + 1)} , ∆X(n) ∼ N (0, σ 2 ) .
¡Paradoxe!
Résolution par simulation des normales ∆Xi
Sans appel au “hasard”, la suite déterministe
(u0 , u1 = Ψ(u0 ), . . . , un = Ψ(un−1 ))
doit ressembler à une suite aléatoire
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Générateur pseudo-aléatoire Générateur pseudo-aléatoire

En R, appel à la procédure

0.0 0.2 0.4 0.6 0.8 1.0

runif( )
Description:
‘runif’ generates random deviates. 500 520 540

uniform sample
560 580 600

Example:
u = runif(20) 1.5

‘[Link]’ is an integer vector, containing the random number

1.0

generator (RNG) state for random number generation in R. It can

0.5

be saved and restored, but should not be altered by the user.

0.0

0.0 0.2 0.4 0.6 0.8 1.0

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Générateur pseudo-aléatoire Générateur pseudo-aléatoire

En C, appel à la procédure

rand() / random() En Scilab, appel à la procédure

SYNOPSIS rand()
# include <stdlib.h>
rand() : with no arguments gives a scalar whose value changes
long int random(void);
each time it is referenced. By default, random numbers are
DESCRIPTION
uniformly distributed in the interval (0,1). rand(’normal’) switches
The random() function uses a non-linear additive feedback random
to a normal distribution with mean 0 and variance 1.
number generator employing a default table of size 31 long
rand(’uniform’) switches back to the uniform distribution
integers to return successive pseudo-random numbers in the range
EXAMPLE
from 0 to RAND MAX. The period of this random generator is
x=rand(10,10,’uniform’)
very large, approximately 16*((2**31)-1).
RETURN VALUE
random() returns a value between 0 and RAND MAX.

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Générateur pseudo-aléatoire Générateur pseudo-aléatoire

Example (Générateur usuel)

Le générateur congruenciel

D(x) = (ax + b) mod (M + 1).

est de période M pour les bons choix de (a, b) et se transforme en

Conclusion :
générateur sur ]0, 1[ par division par M + 2. Utiliser la fonction appropriée sur l’ordinateur ou le logiciel en
v = u*69069069 (1) service plutôt que de construire un générateur aléatoire de
1.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2

mauvaise qualité
0.8
0.6
t+1

0.4
0.2
0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

t
1.0

1.0
0.8

0.8
0.6

0.6
t+10
t+5

0.4

0.4
0.2

0.2
0.0

0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

t t
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (1) Distributions non-uniformes (1)

Autres distributions que la loi uniforme (1) Applications...

Loi binomiale, B(n, p),
Problème réglé en principe puisque
X n 
Théorème (Inversion générique) FX (x) = pi (1 − p)n−i
i
i≤x
Si U est une variable aléatoire uniforme sur [0, 1) et FX est la
fonction de répartition de la variable X, FX−1 (U ) a même loi que X et FX−1 (u) s’obtient numériquement
Preuve. On a Loi exponentielle, E xp(λ),

P (FX−1 (U ) ≤ x) = P (U ≤ FX (x)) = FX (x) FX (x) = 1 − exp(λx) et FX−1 (u) = − log(u)/λ

Note. Si FX n’est pas strictement croissante, on prend

Loi de Cauchy, C (0, 1),
FX−1 (u) = inf {x; FX (x) ≥ u}
1 1
FX (x) = arctan(x)+ et FX−1 (u) = tan(π(u−1/2))
π 2
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (1) Distributions non-uniformes (1)

Autres transformations...

[Indice]
Trouver des transformations reliant la loi d’intérêt et des lois plus Example
simples/mieux connues
Les lois de Student et de Fisher se déduisent naturellement de la
loi normale et de la loi du chi-deux.
Example (Transformation de Box-Müller)
Pour la loi normale N (0, 1), si X1 , X2 ∼ N (0, 1),
i.i.d. Example
La loi de Cauchy se déduit de la loi normale par : si
X12 + X22 ∼ χ22 , arctan(X1 /X2 ) ∼ U ([0, 2π]) i.i.d.
X1 , X2 ∼ N (0, 1), X1 /X2 ∼ C (0, 1)

[Jacobien]
Comme χ22 est identique à E xp(1/2), il vient par inversion
p p
X1 = −2 log(U1 ) sin(2πU2 ) X2 = −2 log(U1 ) cos(2πU2 )
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (1) Distributions non-uniformes (1)

Lois multidimensionnelles

Example
La loi Beta B(α, β), de densité Soit à générer dans Rp
Γ(α + β) α−1 (X1 , . . . , Xp ) ∼ f (x1 , . . . , xp )
fX (x) = x (1 − x)β−1 ,
Γ(α)Γ(β)
dont les composantes ne sont pas nécessairement indépendantes
s’obtient à partir de la loi gamma par: si X1 ∼ G a(α, 1),
X2 ∼ G a(β, 1), alors Cascade rule
X1
∼ B(α, β) f (x1 , . . . , xp ) = f1 (x1 ) × f2|1 (x2 |x1 ) . . . × fp|−p (xp |x1 , . . . , xp−1 )
X1 + X2

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (1) Distributions non-uniformes (2)

Implémentation Autres distributions que la loi uniforme (2)

Simuler pour t = 1, . . . , T
1 X1 ∼ f1 (x1 ) FX−1 rarement disponible
2 X2 ∼ f2|1 (x2 |x1 ) algorithme résident sur machine seulement pour lois usuelles
... lemme d’inversion ne s’applique qu’en dimension 1
nouvelle distribution demandant résolution rapide
p. Xp ∼ fp|−p(xp |x1 , . . . , xp−1 )
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (2) Distributions non-uniformes (2)

Méthode d’acceptation–rejet
Distribution de densité f à simuler Raison :
Loi marginale donnée par
Théorème (fondamental de la simulation)
Z ∞
I0≤u≤f (x) du = f (x)
0

0.25
La loi uniforme sur le sous-graphe et indépendance à la constante de normalisation

0.20
Sf = {(x, u); 0 ≤ u ≤ f (x)}

0.15
Example

f(x)

0.10
a comme loi marginale en x la loi Pour une loi normale, il “suffit” de simuler (u, x) au hasard dans
de densité f .
0.05
{(u, x); 0 ≤ u ≤ exp(−x2 /2)}
0.00

0 2 4 6 8 10

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (2) Distributions non-uniformes (2)

Théorème (Acceptation–rejet)
Algorithme d’acceptation-rejet La variable produite par la régle d’arrêt ci-dessous est distribuée
1 Trouver une densité g simulable telle que suivant la loi fX
f (x) Preuve (1) : On a
sup =M <∞
x g(x) ∞
X
P (X ≤ x) = P (X = Yk , Yk ≤ x)
2 Générer k=1
X∞  k−1
i.i.d. i.i.d. 1
Y1 , Y2 , . . . ∼ g , U1 , U2 , . . . ∼ U ([0, 1]) = 1− P (Uk ≤ f (Yk )/M g(Yk ) , Yk ≤ x)
M
k=1
∞ 
X k−1 Z x Z f (y)/M g(y)
1
= 1− du g(y)dy
3 Prendre X = Yk où M −∞ 0
k=1
X∞  k−1 Z x
k = inf{n ; Un ≤ f (Yn )/M g(Yn )} 1 1
= 1− f (y)dy
M M −∞
k=1
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (2) Distributions non-uniformes (2)

Propriétés
Preuve (2)

5
Fonctionne sans constante de normalisation

4
Ne nécessite pas une borne exacte M
Si (X, U ) est uniforme sur

3
Autorise le recyclage des Yk pour une autre loi f (les Yk
A ⊃ B, la distribution de (X, U )
refusés ne sont plus de loi g)

2
retreinte à B est uniforme sur B.
Demande en moyenne M va Yk pour un X (mesure

1
d’efficacité)

0
−4 −2 0 2 4

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (2) Distributions non-uniformes (2)

Théorème (Enveloppe)
Example S’il existe une densité gm , une fonction gl et une constante M
telles que
Soit f (x) = exp(−x2 /2) et g(x) = 1/(1 + x2 ) gl (x) ≤ f (x) ≤ M gm (x) ,
f (x) 2 √ alors
= (1 + x2 ) e−x /2 ≤ 2/ e
g(x) 1 Générer X ∼ gm (x), U ∼ U[0,1] ;
p
Probabilité d’acceptation e/2π = 0.66 2 Accepter X si U ≤ gl (X)/M gm (X);
3 sinon, accepter X si U ≤ f (X)/M gm (X)
donne des variables aléatoires suivant la loi f .
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Distributions non-uniformes (2) Distributions non-uniformes (2)

Algorithme du rapport d’uniformes

Slice sampler
Example
Résultat :
Simulation uniforme sur
p

0.6
{(u, v); 0 ≤ u ≤ 2f (v/u)}
Pour une loi normale, simuler

0.4
v
produit (u, v) au hasard dans

0.2
X = V /U ∼ f

0.0
Raison : 0.0 0.2 0.4 0.6

u
0.8 1.0 1.2 1.4

Changement de variable (u, v) → (x, u) de Jacobien u et loi

√ −v2 /4u2
√
marginale de x donnée par {(u, v); 0 ≤ u ≤ 2e } = {(u, v); v 2 ≤ −4 u2 log(u/ 2)}
Z √2f (x) p 2
2f (x)
x∼ u du = = f (x)
0 2
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Slice sampler
Slice sampler

0.010
Simuler pour t = 1, . . . , T
Si la simulation uniforme sur

0.008
1 ω (t+1) ∼ U[0,f (x(t))] ;

0.006
G = {(u, x); 0 ≤ u ≤ f (x)} x(t+1) ∼ UG(t+1) , où

f ( x)
2

0.004
est trop compliquée [à cause de l’inversion en x de u ≤ f (x)], on G(t+1) = {y; f (y) ≥ ω (t+1) }.
peut utiliser une marche aléatoire sur G:

0.002
0.000
0.0 0.2 0.4 0.6 0.8 1.0

x
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Justification
Preuve:

Pr((U (t+1) , X (t+1) ) ∈ A × B)

La marche aléatoire se promène uniformément sur G: Z Z Z Z
I0≤u′ ≤f (x) If (x′ )≥u′ (x′ )
si = I0≤u≤f (x) R d(x, u, x′ , u′ )
B A f (x) If (y)≥u′ dy
(U (t) , X (t) ) ∼ UG , Z Z Z
I0≤u′ ≤f (x) If (x′ )≥u′ (x′ )
alors = f (x) R d(x, x′ , u′ )
B A f (x) If (y)≥u′ dy
(U (t+1) , X (t+1) ) ∼ UG . Z Z Z
If (x′ )≥u′ (x′ )
= If (x)≥u′ dx R d(x′ , u′ )
B A If (y)≥u′ dy
Z Z
′ ′
= If (x′ )≥u′ ≥0 d(x , u )
B A

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Example (Loi normale) Note

Pour la loi normale centrée réduite, La méthode fonctionne également si on remplace f par

f (x) ∝ exp(−x2 /2), ϕ(x) ∝ f (x)

un slice sampler est Elle se généralise facilement au cas où f se décompose en

p
Y
ω|x ∼ U[0,exp(−x2 /2)] ,
f (x) = fi (x)
X|ω ∼ U[−√−2 log(ω),√−2 log(ω)] i=1
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

L’algorithme de Metropolis–Hastings
Example (Loi normale tronquée)
Si on considère à la place la loi normale N (−3, 1) tronquée à
[0, 1], de densité
Généralisation du slice sampler à des contextes où le slice sampler
exp(−(x + 3)2 /2) ne peut être utilisé facilement
f (x) = √ ∝ exp(−(x + 3)2 /2) = ϕ(x) ,
2π[Φ(4) − Φ(3)] Idée
un slice sampler est Créer une suite (Xn )n telle que, pour n “assez grand”, la densité
de la loi de Xn soit environ f
ω|x ∼ U[0,exp(−(x+3)2 /2)] ,
X|ω ∼ U[0,1∧{−3+√−2 log(ω)}]

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

L’algorithme de Metropolis–Hastings (2)

Metropolis–Hastings
Partant de X (t) = x(t) ,
1 Générer Yt ∼ q(y|x(t) ).
Soit f la densité de la distribution d’intérêt. On choisit une densité
conditionnelle
2 Prendre
q(y|x) (
Yt avec proba. ρ(x(t) , Yt ),
X (t+1) =
dite instrumentale ou propositionnelle x(t) avec proba. 1 − ρ(x(t) , Yt ),
facile à simuler
où  
partout positive là où f est positive f (y) q(x|y)
ρ(x, y) = min ,1 .
f (x) q(y|x)
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Propriétés Justification

Accepte toujours les déplacements vers des yt tels que

Loi jointe de (X (t) , X (t+1) )
f (yt ) f (xt )
≥ Si X (t) ∼ f (x(t) ),
q(yt |xt ) q(xt |yt )
n
(X (t) , X (t+1) ) ∼ f (x(t) ) ρ(x(t) , x(t+1) ) × q(x(t+1) |x(t) )
Ne dépend pas des constantes de normalisation de f et q(·|x)
[Yt accepté]
(à condition qu’elles soient indépendantes de x) Z h i 
N’accepte jamais les valeurs de yt telles que f (yt ) = 0 + 1 − ρ(x(t) , y) q(y|x(t) ) dy Ix(t) (x(t+1) )
La suite (x(t) )t peut prendre plusieurs fois la même value
[Yt rejeté]
Les X (t) sont des v.a. dépendantes (markoviennes)
liens

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Condition de balance Liens avec le slice sampling

Le slice sampler est un cas (très) particulier d’algorithme de
Metropolis-Hastings où la probabilité d’acceptation vaut toujours 1

 
1 pour la génération de U ,
f (y) q(x|y)
f (x) × ρ(x, y) × q(y|x) = f (x) min , 1 q(y|x) I0≤u′ ≤f (x) f (x)−1 I0≤u′ ≤f (x)
f (x) q(y|x) × =1
I0≤u≤f (x) f (x)−1 I0≤u≤f (x)
= min {f (y)q(x|y), f (x)q(y|x)}
[loi jointe] [loi conditionnelle]
= f (y) × ρ(y, x) × q(x|y)

Donc la loi de (X (t) , X (t+1) ) est la même que celle de 2 pour la génération de X,
(X (t+1) , X (t) ) : si X (t) a la loi f , X (t+1) aussi R
I0≤u≤f (y) I{z;u≤f (z)} (x) {z;u≤f (z)} f (z) dz
× R =1
I0≤u≤f (x) I{z;u≤f (z)} (y) {z;u≤f (z)} f (z) dz
[loi jointe] [loi conditionnelle]
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Propositions indépendantes Propriétés

Loi instrumentale q indépendante de X (t) , dénotée g comme dans

Acceptation-Rejet.
Metropolis-Hastings indépendant Alternative à Acceptation-Rejet
Partant de X (t) = x(t) , Evite le calcul de max f (x)/g(x)
1 Générer Yt ∼ g(y) Accepte plus souvent qu’Acceptation-Rejet
2 Prendre Si xt atteint max f (x)/g(x), presque identique à
 ( ) Acceptation-Rejet

Y f (Yt ) g(x(t) ) Mais la suite des xt n’est pas indépendante
t avec proba. min ,1 ,
X (t+1) = f (x(t) ) g(Yt )

 (t)
x sinon.

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Metropolis–Hastings à marche aléatoire

Example (Loi gamma)
Générer une loi Ga(α, β) à partir d’une loi instrumentale
Ga(⌊α⌋, b = ⌊α⌋/α), où ⌊α⌋ partie entière de α (qui peut être Proposition
générée comme somme d’exponentielles) Yt = X (t) + εt ,

1 Générer Yt ∼ Ga(⌊α⌋, ⌊α⌋/α) où εt ∼ g, indépendant de X (t) , et g symétrique

Loi instrumentale de densité
2 Prendre
 ( )!α−⌊α⌋
 g(y − x)

 Yt x(t) − Yt
Yt avec prob. exp
X (t+1) = x(t) α

 Motivation
x(t) sinon.
Perturbation locale de X (t) / exploration de son voisinage
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Propriétés

Metropolis–Hastings à marche aléatoire

Partant de X (t) = x(t) Monte toujours et descend parfois (cf. algorithmes de
1 Générer Yt ∼ g(y − x(t) ) gradient)
2 Prendre Dépend de la dispersion de g
   Probabilité moyenne d’acceptation

 f (Yt )
Yt
 avec proba. min 1,
f (x(t) )
, Z Z
X (t+1) = ̺= min{f (x), f (y)}g(y − x) dxdy

 [symétrie de g]

x(t) sinon
proche de 1 si g très peu dispersée [Danger!]
loin de 1 si g très dispersée [Re-Danger!]

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Example (Loi normale (2))

Statistiques après 15000 simulations

Example (Loi normale) δ 0.1 0.5 1.0

Générer N (0, 1) fondée sur une perturbation uniforme [−δ, δ] mean 0.399 −0.111 0.10
variance 0.698 1.11 1.06
Yt = X (t) + δωt
Quand δ ↑, exploration plus rapide du support de f .
Probabilité d’acceptation

400
400
250

0.5

0.5

0.5
(t)2

300
(t)
ρ(x , yt ) = exp{(x − yt2 )/2} ∧ 1.
200

300
0.0

0.0

0.0
150

200
200
-0.5

-0.5

-0.5
100

100
100
-1.0

-1.0

-1.0
50

-1.5

-1.5

-1.5
0

0
-1 0 1 2 -2 0 2 -3 -2 -1 0 1 2 3
(a) (b) (c)

3 échantillons avec δ = 0.1, 0.5 et 1.0, avec convergence des

moyennes empiriques (15000 simulations).
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Modèles à données manquantes Principe de complétion

Idée
Cas particulier de modèles où la densité à simuler s’écrit
Simuler f˜ produit des simulations suivant f
Z
f (x) = f˜(x, z)dz Si
Z
(X, Z) ∼ f˜(x, z) ,
La variable Z est alors appelée donnée manquante marginalement
X ∼ f (x)

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Example (Mélange de lois)

Data Augmentation Soit à simuler sur R2 la loi de densité proportionnelle à
Partant de x(t) , 100 n o
2 2
Y 2 /2 2 /2
1. Simuler Z (t+1) ∼ f˜Z|X (z|x(t) ) ; e−µ1 −µ2 × 0.3 e−(xi −µ1 ) + 0.7 e−(xi −µ2 )
2. Simuler X (t+1) ∼ f˜X|Z (x|z (t+1) ) . i=1

où les xi sont donnés et (µ1 , µ2 ) représente les coordonnées de la

variable aléatoire
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Echantillon de 0.3 N(2.5,1)+ 0.7 N(0,1)

0.35

3.0
Complétion (1)
0.30

Remplacer chaque somme dans la densité par une intégrale:

2.8
0.25

Z 
2 2
0.3 e−(xi −µ1 ) /2 + 0.7 e−(xi −µ2 ) /2 = I[0,0.3 e−(xi −µ1 )2 /2 ] (ui )

2.6
0.20

µ2
0.15

+I[0.3 e−(xi −µ1 )2 /2 ,0.3 e−(xi −µ1 )2 /2 +0.7 e−(xi −µ2 )2 /2 ] (ui ) dui

2.4
0.10

2.2
et simuler ((µ1 , µ2 ), (U1 , . . . , Un )) = (X, Z) par Data
0.05

Augmentation
0.00

2.0

−2 0 2 4 −0.4 −0.2 0.0 0.2 0.4

x µ1

Histogramme des xi et surface de la loi de (µ1 , µ2 ) associée

Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE) Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires Simulation de variables aléatoires
Méthodes de Markov Méthodes de Markov

Complétion (2)
Remplacer les Ui par les ξi , où Conditionnement (1)
( 2
1 si Ui ≤ 0.3 e−(xi −µ1 ) /2 , La loi conditionnelle de Z = (ξ1 , . . . , ξn ) sachant X = (µ1 , µ2 ) est
ξi = donnée par
2 sinon
2
0.3 e−(xi −µ1 ) /2
Alors Pr (ξi = 1|µ1 , µ2 ) =
0.3 e−(xi −µ1 )2 /2 + 0.7 e−(xi −µ2 )2 /2
2
0.3 e−(xi −µ1 ) /2 = 1 − Pr (ξi = 2|µ1 , µ2 )
Pr (ξi = 1|µ1 , µ2 ) =
0.3 e−(xi −µ1 )2 /2 + 0.7 e−(xi −µ2 )2 /2
= 1 − Pr (ξi = 2|µ1 , µ2 )
Nouveaux outils informatiques pour la Statistique exploratoire (=NOISE)
Simulation de variables aléatoires
Méthodes de Markov

Conditionnement (2)
La loi conditionnelle de X = (µ1 , µ2 ) sachant Z = (ξ1 , . . . , ξn ) est
donnée par
2 2
Y 2
Y 2
(µ1 , µ2 )|Z ∼ e−µ1 −µ2 × e−(xi −µ1 ) /2 × e−(xi −µ2 ) /2
{i;ξi =1} {i;ξi =2}
(   )
n1 µ̂1 2
∝ exp −(n1 + 2) µ1 − /2
n1 + 2
(   )
n2 µ̂2 2
× exp −(n2 + 2) µ2 − /2
n2 + 2

où nj est le nombre de ξi égaux à j et nj µ̂j est la somme des xi

associés à ces ξi égaux à j
[Easy!]