École Normale Supérieure de Rennes Angelo Rosello

2017/2018 Sciences Fondamentales 2a

TP1 : Chaines de Markov, algorithme de Metropolis

1 Simuler une chaine de Markov
On considère l’ensemble E = {1, ..., N }, N ≥ 1.
p = (p1 , ...pN ) désigne une probabilité sur E.
1. En utilisant la commande rand (voir help rand), écrire une fonction tirage(p) qui
simule le tirage d’une variable aléatoire sur E de loi p.
2. Ecrire une fonction markov(P,mu,t) qui prend en argument une matrice stochastique
de taille N , une probabilité µ sur E et un entier t ≥ 1, et qui simule le tirage de
(X0 , X1 , ..., Xt ), où X est une CM (P, µ).
3. Exemple : tracer les 50 premiers sauts d’un tirage de la chaine de Markov suivante, issue
de 1,
1/3

1 2

1/2
2/3 1 1/2

1/2
1/2 3 4

4. Pourquoi peut-on introduire l’unique loi invariante µ de cette chaine ?

µ (vecteur ligne) est un vecteur propre à gauche pour P , associé à la valeur propre 1 :

µP = µ

autrement dit, t µ (vecteur colonne) est un vecteur propre à droite pour t P :

t
P t µ =t µ.

Déterminer µ en utilisant [V,D]=eig(P’) (voir help eig). Attention, µ doit être une
probabilité : n’oubliez pas de renormaliser pour que 4i=1 µi = 1 ! On pourra utiliser la
P

fonction sum.
5. Que sait-on du comportement de P(Xn = i) quand n → ∞ ?
On veut illustrer ce résultat de convergence.

a) En calculant les probabilités "réelles"

Que vaut P(Xn = i) en fonction de la matrice de transition P ? (on rappelle que la chaîne
démarre de l’état 1). Tracer l’évolution de P(Xn = i) pour 0 ≤ n ≤ 50 et pour chaque
i ∈ {1, 2, 3, 4}. Conclure.

b) En estimant ces probabilités (OPTIONNEL)

En utilisant markov(P,mu,t), générer K tirages (indépendants) des t premiers sauts de
la chaine :
(X0k , X1k , ..., Xtk ), 1 ≤ k ≤ K.

1
 On prendra K = 2000, t = 50. On définit ensuite
K
1 X
p̂n (i) = 1 k 0 ≤ n ≤ t, i ∈ {1, 2, 3, 4}.
K k=1 Xn =i

On pourra utiliser la commande (a==b) qui vaut 1 si a = b, 0 sinon.

En vertu de quel théorème probabiliste peut-on dire que p̂n (i) est une bonne estima-
tion de P(Xn = i) ?

Tracer l’évolution de p̂n (i) quand 0 ≤ n ≤ t, pour chacun des i ∈ {1, 2, 3, 4} et conclure.

2 Metropolis-Hastings
"Loi discrète"
On se donne µ une probabilité sur un ensemble E. On rappelle le principe de l’algorithme de
Metropolis-Hastings, partant de la matrice de transition Q d’une chaine de Markov irréductible,
apériodique sur E, dans le cas simplifié où Q est symétrique :

∀x, y ∈ E, Q(x, y) = Q(y, x).

Initialiser X0 ∈ E à une valeur quelconque.

Sachant qu’au temps n, Xn = x ∈ E,
1. Tirer un y ∈ E au hasard selon la loi : P(y) = Q(x, y)
(transition donnée par la matrice Q)
2. Si µ(y) ≥ µ(x), poser Xn+1 = y. Sinon :
(
y avec proba µ(y)/µ(x)
Xn+1 =
x avec proba 1 − µ(y)/µ(x)

1. On travaille avec E = {1, ...N }, N ≥ 1.

Ecrire une fonction metropolis(Q,mu,n) qui prend en argument une matrice stochas-
tique de taille N , une probabilité µ sur E et un entier n ≥ 1, et qui simule le tirage de
Xn , où X est la chaine correspondant à l’algorithme de métropolis.
2. Exemple : On considère la matrice Q correspondant à la marche aléatoire isotrope sur
un "tore" à 5 éléments :
1/2 1/2
1 2 3

1/2 1/2
5 4
1/2

Cette chaine est-elle irréductible ? Apériodique ?

On considère la probabilité µ = (0.1, 0.2, 0.1, 0.3, 0.3), et n = 100.

2
 On veut vérifier que la loi de Xn (où X correspond à l’algorithme de métropolis) est
proche de µ.
En utilisant metropolis(Q,mu,n), générer K tirages (indépendants)
Xnk , 1 ≤ k ≤ K
et définir
K
1 X
p̂n (i) = 1 k i ∈ {1, 2, 3, 4, 5}.
K k=1 Xn =i
qui est un estimateur de P(Xn = i). On prendra K = 1000.
Comparer p̂n (i) à µi pour 1 ≤ i ≤ 5, conclure.

On pourra (si on a du temps) réitérer cette procédure avec différents choix de proba-
bilités µ pour se convaincre de l’algorithme de Metropolis-Hastings joue bien son rôle.

"Loi continue"
On souhaite maintenant simuler une variable aléatoire réelle, admettant la densité de pro-
babilité f : Z
∀I intervalle de R, P(X ∈ I) = f (y)dy
I
R +∞
où f : R → R+ satisfait −∞ f (y)dy = 1.

Evidemment, en machine, R est remplacé par un gros ensemble fini, de pas ε > 0 très petit :
E = [ - C : epsilon : C ] = {−C, −C + ε, −C + 2ε, ..., C − ε, C}.
Le problème se ramène alors à simuler une variable aléatoire à valeurs dans E, de loi µ :
f (x) X
µ(x) = où Z = f (x).
Z x∈E

En effet, si X suit la loi µ, pour tout intervalle I,

P P
f (x) ε f (x)
µ(x) = P x∈I = P x∈I
X
P(X ∈ I) =
x∈I x∈E f (x) ε x∈E f (x)
P R
et ε x∈I f (x) est une somme de Riemann qui approche I f (y)dy, donc
R
f (y)dy
Z
P(X ∈ I) ' R I = f (y)dy.
R f (y)dy I

Pour simuler la loi µ, on propose l’algorithme de Metropolis suivant :

Initialiser X0 = 0.
Sachant qu’au temps n, Xn = x,
1. Tirer un y ∈ R au hasard selon la loi :
y =x+u
où u suit une loi U([−δ, δ]) : loi uniforme sur l’intervalle [−δ, +δ].
2. Si f (y) ≥ f (x), poser Xn+1 = y. Sinon :
(
y avec proba f (y)/f (x)
Xn+1 =
x avec proba f (y)/f (x)

3
1. Comment simuler simplement un tirage de u de loi U([−δ, δ]) ? (On rappelle qu’on dispose
de la fonction rand...)
2. Ecrire une fonction metropolis_densite(f,delta,n) qui prend en argument une den-
sité de probabilité f , un δ > 0 et un entier n ≥ 1, et qui simule le tirage de Xn , où X est
la chaine correspondant à l’algorithme que l’on vient de décrire.
3. Choisissons la densité gaussienne
1 x2
f (x) = √ e− 2 .
2π
c’est à dire f=@(x) (1/sqrt(2*pi))*exp(-x.^2/2).
On souhaite vérifier que, quand n est grand, la loi de Xn est proche de la loi admettant
la densité f . Fixons δ = 1, n = 50.

En utilisant metropolis_densite(f,delta,n), générer K tirages (indépendants)

Xnk , 1 ≤ k ≤ K.

On prendra K = 2000.
On construit ensuite un histogramme des tirages X = (Xnk )k=1...N :

— Poser I=[-5:h:5] (on prendra h = 0.1 ou 0.05...).

— H=histc(X,I) renvoie un vecteur comptant le nombre de Xnk appartenant à chacun

des intervalles délimités par I (voir help histc).
Afficher l’histogramme avec bar(I,H,’w’) : il n’est pas normalisé.

— On renormalise l’histogramme pour que l’aire totale vale 1 : H=H./(h*sum(H)).

Afficher l’histogramme renormalisé et comparer avec le graphe de f .

On pourra (si on a du temps) réitérer cette procédure avec d’autres choix de densité :

1 x2
Loi normale de variance σ 2 : f (x) = √ e− 2σ2
2πσ
1
Loi de Laplace : f (x) = e−|x|
2
1 1
Loi de Cauchy : f (x) =
π 1 + x2
On pourra également jouer sur les paramètres n, δ, K...