PRÉCIS DE SIMULATION

P. Del Moral
Centre INRIA Bordeaux Sud-Ouest
& Institut de Mathématiques de Bordeaux
Université Bordeaux I, 351, cours de la Libération
33405 Talence, France
Table des matières

1 La notion de variable aléatoire 5

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Le principe d’indépendance . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Les générateurs aléatoires . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 Quelques lois usuelles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.1 Les lois uniformes . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.2 Les lois d’événements composés . . . . . . . . . . . . . . . . . . . 13
1.4.3 Les lois géométriques . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4.4 Les lois discrètes . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.4.5 Les lois continues . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2 Quelques outils de simulation 27

2.1 La méthode d’inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.1 Le principe de base . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.2 Retour aux lois géométriques . . . . . . . . . . . . . . . . . . . . 29
2.1.3 Les lois de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.1.4 Les lois de Rayleigh-Weibull . . . . . . . . . . . . . . . . . . . . . 32
2.2 Changements de variables . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.2 Lois uniformes sur des cubes . . . . . . . . . . . . . . . . . . . . 34
2.2.3 Lois uniformes sur des surfaces . . . . . . . . . . . . . . . . . . . 36
2.2.4 Lois uniformes sur des disques . . . . . . . . . . . . . . . . . . . 41
2.2.5 L’algorithme Box-Muller . . . . . . . . . . . . . . . . . . . . . . . 44
2.3 Méthode d’acceptation-rejet . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3.1 Les lois de référence . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3.2 Les taux d’acceptation . . . . . . . . . . . . . . . . . . . . . . . . 50
2.3.3 Vérification mathématique . . . . . . . . . . . . . . . . . . . . . 51
2.4 Files d’attentes exponentielles . . . . . . . . . . . . . . . . . . . . . . . . 53
2.4.1 La loi Gamma . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.4.2 La loi de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.4.3 La statistique d’ordre uniforme . . . . . . . . . . . . . . . . . . . 57

1
3 Chaı̂nes de Markov discrètes 61
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.2 Chaı̂nes de Markov discrètes . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.2.1 Semigroupes de transitions . . . . . . . . . . . . . . . . . . . . . 64
3.2.2 Processus historique . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.2.3 Interprétation matricielle . . . . . . . . . . . . . . . . . . . . . . 67
3.3 Quelques Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.3.1 Files d’attentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.3.2 Modèle d’urne . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.3.3 Marche aléatoire sur Z . . . . . . . . . . . . . . . . . . . . . . . . 72
3.3.4 Marche aléatoire sur Zd . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.5 Marche aléatoire arrétée . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.6 Processus de branchements . . . . . . . . . . . . . . . . . . . . . 75

4 Chaı̂nes de Markov abstraites 79

4.1 Description des modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.1.1 Semigroupe des transitions . . . . . . . . . . . . . . . . . . . . . 80
4.1.2 Équations de Chapman-Kolmogorov . . . . . . . . . . . . . . . . 82
4.1.3 Processus historique . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Chaı̂nes linéaires et gaussiennes . . . . . . . . . . . . . . . . . . . . . . . 84
4.2.1 Formulation canonique . . . . . . . . . . . . . . . . . . . . . . . . 85
4.2.2 Formulation dynamique . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 Processus de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.4 Évolutions dans des milieux absorbants . . . . . . . . . . . . . . . . . . 87

5 Chaı̂nes de Markov non linéaires 91

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.2 Description des modèles . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3 Interprétations particulaires en champ moyen . . . . . . . . . . . . . . . 93
5.4 Champs moyens de type gaussien . . . . . . . . . . . . . . . . . . . . . . 94
5.5 Modèles simplifiés de gaz de McKean . . . . . . . . . . . . . . . . . . . . 96
5.6 Flots de mesures de Feynman-Kac . . . . . . . . . . . . . . . . . . . . . 96
5.6.1 Description des modèles . . . . . . . . . . . . . . . . . . . . . . . 96
5.6.2 Chaı̂nes de Markov non linéaires . . . . . . . . . . . . . . . . . . 99
5.6.3 Champs moyens de type évolutionnaire . . . . . . . . . . . . . . 100

6 L’équation de la chaleur 103

6.1 Les fluctuations browniennes . . . . . . . . . . . . . . . . . . . . . . . . 103
6.2 La loi des grands nombres . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.3 Marches aléatoires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.4 L’équation de la chaleur . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.5 Une formulation faible . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
7 Dynamiques de population avec branchements 113
7.1 Processus de branchements spatio-temporels . . . . . . . . . . . . . . . . 113
7.2 Algorithme génétique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7.2.1 Sélection/Adaptation . . . . . . . . . . . . . . . . . . . . . . . . 117
7.2.2 Mutation/Exploration . . . . . . . . . . . . . . . . . . . . . . . . 117
7.3 Modèles d’arbres généalogiques . . . . . . . . . . . . . . . . . . . . . . . 119
7.3.1 Modèles non homogènes . . . . . . . . . . . . . . . . . . . . . . . 119
7.3.2 Modèles trajectoriels . . . . . . . . . . . . . . . . . . . . . . . . . 120
7.4 Chaı̂nes renforcées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
Chapitre 1

La notion de variable aléatoire

1.1 Introduction
Le mot latin aléa, et le mot arabe hasard se réfèrent tous deux au “jeu de dès”. Cette
métaphore antique transmet ainsi l’idée que le résultat du lancer de dés est totalement
inattendu et imprévisible. La formalisation mathématique de cette notion de hasard
correspond à la notion de variable aléatoire. Ce nom probabiliste est un vrai pléonasme
mathématique combinant la variabilité d’une quantité est les fluctuations aléatoires des
résultats observés.
L’archétype des variables aléatoires est le simple jeu du pile ou face. Le résultat
X du lancer d’une pièce de monnaie est clairement l’une des deux faces, mais il est
bien difficile de prévoir laquelle. Si l’on désigne par le nombre 1 le coté pile, et par
le nombre 2 le coté face, les valeurs possibles du lancer sont données un ensemble
réduit à deux points {1, 2}. Lorsque la pièce n’est pas truquée, les deux résultats sont
équiprobables. Autrement dit, la probabilité pour que X prenne la valeur 1 coincide
avec la probabilité pour que X prenne la valeur 2. Plus formellement, si l’on note
respectivement par P(X = 1) et P(X = 2), ces deux probabilités, nous avons

P(X = 1) = P(X = 2) = 1/2

Cette formule mathématique exprime le fait que chaque face se réalise avec une
probabilité 1/2.
Le lecteur profane est en droit de se poser au moins un couple de questions
naturelles :
– Que signifie ces grandeurs numériques ?
– Existe-t-il des pièces si parfaites ?
Comme nous l’avons souligné dans l’introduction, les modèles probabilistes sont
toujours des modèles simplifiés, mais suffisament précis pour prédire les chances
d’obtenir une réalisation donnée. Pour être plus précis, supposons que l’on lance de
nombreuses fois la même pièce de monnaie. On note scrupuleusement la proportion de

5
fois où le cote pile s’est réalisé. Après un certain nombre de lancers, ces proportions se
stabilisent autour d’une valeur déterministe comprise dans l’intervalle [0, 1]. Ce nombre
correspond à la probabilité pour que le résultat du lancer soit “pile”.
Une probabilité est donc par essence une notion asymptotique, évaluant
le nombre de chances d’obtenir une réalisation d’un évènement aléatoire donné. Si les
fréquences de réalisation du coté pile convergent vers 1/100, lorsqu’on augmente le
nombre de lancers, c’est que la pièce utilisée est truquée. Le coté “face” se réalise 99
fois sur 100. Cela ne veut bien évidemment pas dire que sur 100 expériences, on aura
exactement 99 fois le coté “face” ! Néanmoins sur de longues séquences, on est à peu
sur que le coté “face” se réalisera 99 fois plus souvent.
L’existence d’une pièce parfaite est aussi utopique que l’existence d’une pièce pour
laquelle la probabilité de réalisation du coté “pile” est exactement égale à 1/π. Une
probabilité rationnelle, pouvant s’exprimer sous forme de fraction d’entiers, est elle
plus vraisemblable qu’une quelconque probabilité réelle ?
Cette simple question concernant l’existence d’une telle pièce nous conduit vers une
région mathématique encore plus mystérieuse, celle des nombres réels.
La continuité et l’extension infinie de notre univers restent deux des grandes
questions des sciences modernes. L’existence des nombres réels semble dictée par
l’analyse mathématique. Plus précisément, les nombres rationnels sont d’un usage
courant, mais leur ensemble est trop troué pour définir convenablement certains objects
limites, tels des vitesses ou des accélérations de mouvement, des volumes ou encore des
aires de surfaces planes. Sans ces nombres ils serait par exemple totalement impossible
de trouver le périmètre d’un cercle de un mètre diamètre !
Une fois admis l’existence de ces nombres réels, il devient très aisé d’imaginer des
phénomènes aléatoires à valeurs dans des espaces continus : lancers de flechettes sur
une disque, fluctuations de températures d’un liquide, perturbations de mesures dans
des capteurs électroniques, évolutions de cyclones, instants d’arrivée de clients dans
une file d’attente, etc.
L’exemple le plus simple de phénomène aléatoire réel est celui du choix d’un point
au hasard dans l’intervalle réel [0, 1]. Pour poursuivre notre discussion, on note par U
le résultat aléatoire d’une telle expérience. La probabilité de choisir un nombre entre 0
et 1/2 est égale à la probabilité de choisir ce nombre entre 1/2 et 1

P(U ∈ [0, 1/2]) = P(U ∈ [1/2, 1]) = 1/2

De même, la La probabilité de choisir un nombre entre 0 et 1/4 est égale à la probabilité

de choisir ce nombre entre 1/4 et 1/2, elle même égale à la probabilité de le choisir entre
1/2 et 3/4, ou encore entre 3/4 et 1

P(U ∈ [0, 1/4]) = P(U ∈ [1/4, 1/2]) = P(U ∈ [1/2, 3/4]) = P(U ∈ [3/4, 1]) = 1/4

Plus généralement, la probabilité de choisir un nombre dans un intervalle donné [a, b] ⊂

[0, 1], dépend uniquement de la longueur de ce dernier
Z b
P(U ∈ [a, b]) = dx = (b − a)
a

Ainsi la probabilité de choisir U dans [1/2, 1/2 + 1/1000] est de 1/1000, celle de le
choisir dans [1/2, 1/2 + 1/10000] est de 1/10000. De proche en proche, on montre qu’il
est de plus en plus improbable de le choisir de plus en plus proche de 1/2
→0
P(U ∈ [1/2, 1/2 + ]) −−−−→ P(U = 1/2) = 0

Ce résultat au coeur de la théorie de l’intégration de Lebesgue est assez déroutant. Si

l’on mesure une réalisation d’un phénomène uniforme sur [0, 1], tel le choix au hasard
d’une pièce de monnaie truquée ou non, on observera clairement la réalisation d’un
nombre quelconque entre 0 et 1. Cependant, la probabilité d’observer ce même nombre
lors d’une telle expérience est tout simplement nulle ! Ces évènements élémentaires n’ont
donc aucune chance de se réaliser. Cette propriété mathématique peut être rapprochée
de la philosophie des “matérialistes quantiques” soutenant le fait que “l’Univers n’existe
pas indépendamment de tout acte d’observation” (Fritz Rohrlich, Science 1983).
Ces pensées sont en partie soutenues par le principe d’incertitude d’Heisenberg
énoncé dès 1927. Ce principe d’indétermination est à la base des sciences physiques
modernes. Il souligne qu’il est impossible de connaitre avec précision, et à un
instant donné, le couple vitesse-position d’une particule à l’échelle atomique. Plus
précisément, les incertitudes de mesures sur l’une des deux composantes sont
inversement proportionnelles au degré de connaissance sur la seconde. En ce sens,
certains aspects déroutant de la théorie des probabilités sont très probablement le
reflet de certaines propriétés naturelles et physiques.
Les variables aléatoires uniformes sur l’intervalle [0, 1], sont le ciment des méthodes
de simulation numérique de lois de probabilités complexes. Plus précisément, on peut
démontrer que tout phénomène aléatoire à valeurs réelles s’exprime en terme
des variables uniformes sur [0, 1]. Autrement dit, on peut simuler toute loi de
probabilité sur un espace multi-dimensionnel Rd , avec d ≥ 1, par une suite finie de
variables aléatoires uniformes sur [0, 1]. Pour s’en convaincre, on peut noter que le
simple jeu de pile ou face examiné plus haut, peut être décrit par la formule suivante :

1 si U ∈ [0, 1/2[
X=
2 si U ∈ [1/2, 1]

Autrement dit, lorsque l’évènement U ∈ [0, 1/2[ se réalise, on convient que la pièce
de monnaie est tombée sur le coté pile, et sur le coté face dans le cas contraire. Par
construction, la probabilité d’obtenir le coté “pile” est bien à nouveau égale à 1/2 :

P(X = 1) = P(U ∈ [0, 1/2[) = 1/2

1.2 Le principe d’indépendance
Les variables aléatoires uniformes sur [0, 1] permettent de construire très simplement
des phénomènes aléatoires dans l’espace. Un point U = (U1 , U2 ) choisi au hasard dans
le carré unité C = ([0, 1] × [0, 1]) correspond tout simplement au choix de deux variables
U1 et U2 uniformes et indépendantes sur [0, 1]. La probabilité pour que U soit choisi
dans un petit pavé

P = ([a1 , a2 ] × [b1 , b2 ]) ⊂ C = ([0, 1] × [0, 1])

est le produit des probabilités pour que chacune des deux coordonnées U1 et U2 , tombe
respectivement sur les faces [a1 , a2 ] et [b1 , b2 ]. Par conséquent, la probabilité de choisir
un point U = (U1 , U2 ) au hasard dans le carré P est alors donnée par l’aire du pavé :

P(U ∈ P) = P(U1 ∈ [a1 , a2 ] et U2 ∈ [b1 , b2 ])

= P(U1 ∈ [a1 , a2 ]) × P(U2 ∈ [b1 , b2 ]) = (a2 − a1 ) × (b2 − b1 )

On peut étendre ces constructions au cube unité

V = ([0, 1] × [0, 1] × [0, 1]) = [0, 1]3

Choisir au hasard un point U = (U1 , U2 , U3 ) dans ce cube revient trois variables U1 , U2 ,

et U3 uniformes et indépendantes sur [0, 1]. La probabilité de choisir ce point dans
un pavé
P = ([a1 , a2 ] × [b1 , b2 ] × [c1 , c2 ])
est alors donnée par l’aire du pavé :

P(U ∈ P) = P(U1 ∈ [a1 , a2 ] et U2 ∈ [b1 , b2 ] et U3 ∈ [c1 , c2 ])

= P(U1 ∈ [a1 , a2 ]) × P(U2 ∈ [b1 , b2 ]) × P(U3 ∈ [c1 , c2 ])
= (a2 − a1 ) × (b2 − b1 ) × (c2 − c1 )

Dans notre discussion, nous avons utilisé sans y préter trop d’attention une des
notions fondamentales de la théorie des probabilités : la notion d’indépendance
entre des évènements aléatoires. Intuitivement, deux évènements aléatoires sont
physiquement indépendants lorsque la réalisation de l’un n’a aucune influence sur les
possibilités de réalisation du second.
Reprenons l’exemple du jeu de pile ou face. On note X1 et X2 les résultats de deux
lancers successifs. L’indépendance entre les lancers se traduit par le fait que le résultat
du premier lancer n’influe pas sur la réalisation du second. Par exemple, le fait d’avoir
“pile” au premier lancer, n’augmente pas les chance d’avoir à nouveau “pile” au second.
Les deux évènements

A1 = (X1 = 0) et A2 = (X2 = 0)
sont donc physiquement indépendants. Examinons de plus près la caratérisation
mathématique de cette notion d’indépendance. Nous commençons par noter que toutes
les réalisations possibles et envisageables sont équiprobables. Plus précisément, nous
avons autant de chance de voir deux fois le coté “pile”, que deux fois le coté ”face”, ou
encore deux cotés opposés dans un ordre quelconque :

P((X1 , X2 ) = (0, 0)) = P((X1 , X2 ) = (0, 1))

= P((X1 , X2 ) = (1, 0)) = P((X1 , X2 ) = (1, 1)) = 1/4

On a donc une chance sur quatre pour que les deux évènements A1 et A2 se réalisent :

P(A1 et A2 ) = P((X1 , X2 ) = (0, 0)) = 1/4

D’autre part, puisque chacun d’entre eux à une chance sur deux de se réaliser, on a
aussi
P(A1 ) × P(A2 ) = P(X1 = 0) × P(X2 = 0) = 1/2 × 1/2 = 1/4
Par conséquent, la probabilité pour que A1 et A2 se produisent simultanément est égale
au produit des probabilités pour que chacun d’eux se soit séparément réalisé

P(A1 et A2 ) = P(A1 ) × P(A2 )

En termes mathématiques, deux évènements vérifiant cette propriété sont dit

indépendants. Cette propriété élémentaire permet de calculer très simplement les
probabilités d’évènements parfois complexes. Par exemple, choisir successivement, et
au hasard une longue suite de points (U1 , U2 , . . . , Un ) dans le segment [0, 1], s’effectue
en n étapes indépendantes. On choisit tout d’abord la première coordonnée U1 ,
puis la seconde U2 , puis la troisième U3 , jusqu’à la nième Un dans l’intervalle [0, 1].
La probabilité pour que ces variables tombent successivement dans des segments
[a1 , b1 ],. . . , [an , bn ] ⊂ [0, 1] est tout simplement donnée par la formule produit

P(U1 ∈ [a1 , b1 ] , U2 ∈ [a2 , b2 ] , . . . , Un ∈ [an , bn ])

= P(U1 ∈ [a1 , b1 ]) × P(U2 ∈ [a2 , b2 ]) × . . . × P(Un ∈ [an , bn ])

= (b1 − a1 ) × (b2 − a2 ) × . . . × (bn − an )

1.3 Les générateurs aléatoires

Simuler sur ordinateur des variables aléatoires uniformes et indépendantes sur
[0, 1], revient à générer une séquence de nombres u = (un )n≥0 compris entre 0 et
1 qui ressemble à une suite de variables aléatoires U = (Un )n≥0 indépendantes et
uniformes sur [0, 1]. Cette suite u doit avoir les mêmes propriétés statistiques que U .
Par exemple, la séquence de nombres u doit donc se dérouler dans le temps de façon
totalement imprévisible. Si la suite u était périodique et de période trop faible, on
pourrait trop facilement prévoir ses valeurs après cette période, La suite u doit aussi
être uniformément répartie sur l’intervalle [0, 1]. On s’attend donc à avoir la moitié des
nombres dans [0, 1/2] et plus généralement il devrait y avoir une proportion (b − a) de
nombres dans tout intervalle [a, b] ⊂ [0, 1], etc.

Il est important de souligner qu’il n’existe à ce jour aucune méthode exacte pour
générer une telle suite de nombres u. Pire encore, la plupart des générateurs de
nombres uniformes connus sont basés sur des systèmes déterministes chaotiques. Ces
algorithmes sont loin d’être aléatoire ! Ils possèdent toujours une période mais cette
dernière est suffisamment grande pour ne pas être vue par l’utilisateur. Ces systèmes
s’expriment le plus souvent par une suite arithmétique linéaire congruentielle de la
forme suivante :

un = xn /m avec xn = (a xn−1 + b) mod(m)

La valeur initiale u0 est appelée la graine, et les paramètres (a, b, m) sont des entiers
bien choisis. Lorsque

a = 5, b = 1, , m = 24 et x0 = 1

la suite xn = (5 xn−1 + 1) mod(16) obtenue par cet algorithme est donnée par

1, 6, 15, 12, 13, 2, 11, 8, 9, 14, 7, 4, 5, 10, 3, 0, 1, 6, 15, 12, 13, 2, 11, 8, 9, 14, . . .

La période de cette suite est égale à 16. Cette période est bien entendu maximale car
l’opération “modulo 16” ne permet de générer que des nombres compris entre 0 et 16.
Lorsque b = 0, 0 < x0 < m, et m est un nombre premier, le petit théorème de Fermat
indique que la période de la suite est égale à (m − 1).
Il existe une foule de générateurs pseudo-aléatoires fondés sur ces opérations
arithmétiques, quasiment autant que de choix des valeurs (a, b, m) ! Pour donner
quelques exemples, certains logiciels utilisent

(a, b, m) = (1313 , 0, 259 ) ou (a, b, m) = (75 , 0, 231 − 1)

Le générateur Unix rand48 correspond aux valeurs

a = 25214903917, b = 11 et m = 248

et celui du CRAY

a = 44485709377909, b=0 et m = 248

On peut aussi trouver dans la littérature d’autres générateurs à programmer soi-même

([?, ?]).
 Il existe de nouveaux types de générateurs basés sur les lois d’incertitudes inhérentes
à la mécanique quantique. L’idée prometteuse de ces générateurs est à peu près la
suivante. Le noyau de certaines particules telles le Krypton-85 contient une horloge
aléatoire. Lorsque cette dernière sonne, un éléctron est expulsé. Les règles de la
physique quantique montrent qu’il est impossible de prédire ces instants. Ces horloges
physiques sont par essence aléatoires ! Les futurs ordinateurs quantiques seront capable
de “simuler” de vrais aléas physiques. Mieux encore, la physique quantique prédit que
les différents états quantiques cohérents se réalisent dans des dimensions multiples. Ces
ordinateurs ultrapides seraient donc en mesure de traiter en parallèle des informations
variées dans chaque univers !

1.4 Quelques lois usuelles

Comme leur nom l’indique, les variables discrètes représentent des phénomènes
aléatoires ne prenant qu’un nombre fini, ou au plus dénombrable de valeurs. Les
archétypes de ces variables sont bien entendu les jeux de “pile ou face”, le “lancé
de dés”, ou encore les célèbres “jeux de loto” nationaux, ou européens.
Ces variables discrètes permettent aussi de formaliser tout type de succession
d’évènements aléatoires : tirage de boules dans une urne, séquences de pannes de
machines, suites d’échecs ou de succés lors d’une expérience répétée d’évènements,
sélection aléatoire d’objets, nombre de naissances ou de décés d’individus sur des sites
donnés, nombres d’accidents d’avions par année, etc.
Ces phénomènes aléatoires discrets obéissent à des règles probabilistes simples et
précises. Ces dernières sont bien souvent de nature combinatoire, ce qui n’a rien de très
séduisant ni le profane, ni pour la plupart des probabilistes avisés.

1.4.1 Les lois uniformes

Les phénomènes aléatoires les plus simples à décrire sont ceux ne prenant
qu’un nombre fini de valeurs, toutes équiprobables. Ces phénomènes formalisés
correspondent à la notion de variables discrètes et uniformes. Aucune de ces valeurs
possibles n’est privilégiée ! La description mathématique de leur lois est “simplement”
de type combinatoire.
Examinons pour commencer le cas des variables aléatoires uniformes X sur un
espace réduit à deux points E = {1, 2}. Le simple jeu du “pile ou face” (avec un pièce
non truquée) peut bien entendu être décrit par de telles variables. Le résultat du lancer
X peut alors prendre l’une des deux valeurs 1 ou 2, avec la même probabilité
P(X = 1) = P(X = 2) = 1/2
Pour générer sur ordinateur un tel choix, on utilise une variable uniforme U sur [0, 1],
et l’on pose simplement
X = 1 × 1[0,1/2[ (U ) + 2 × 1[1/2,1] (U ) (1.1)
où 1A désigne la fonction indicatrice d’ un ensemble A

1 si u ∈ A
1A (u) =
0 si u 6∈ A

L’expression (1.1) permet de détecter lequel des deux évènements

A1 = (U ∈ [0, 1/2[) ou A2 = (U ∈ [1/2, 1])

s’est produit lors de la simulation de la variable U . Dans le cas du jeu de “pile ou face”,
si le nombre U tombe dans l’intervalle [0, 1/2), on convient que le résultat du lancer
est “pile” (et “face” dans le cas contraire).
On peut éviter de tester si la simulation du nombre U est tombée soit à gauche, soit
à droite de 1/2, en posant tout simplement

X = 1 + [2U ]

où [a] désigne la partie entère d’un nombre réel a. Il est aisé de vérifier les équivalences
entre les évènements suivantes

(0 ≤ U < 1/2) ⇐⇒ (X = 1) et (1/2 ≤ U < 1) ⇐⇒ (X = 2)

Dans le le jeu du “lancé de dés” avec un dés non truqué, le résultat du lancer X
peut prendre l’une des six valeurs 1, 2, 3, 4, 5, ou 6 avec la même probabilité

P(X = 1) = P(X = 2) = P(X = 3) = P(X = 4) = P(X = 5) = P(X = 6) = 1/6

Pour générer sur ordinateur un lancer de dés aléatoire, on utilise une seule variable
uniforme U sur [0, 1], et l’on pose simplement

X = 1 + [6U ]

Comme précédemment, nous avons les équivalences évènementielles

        
1 1 2 5
U ∈ 0, ⇔ (X = 1) , U ∈ , ⇔ (X = 2) , . . . , U ∈ ,1 ⇔ (X = 6)
6 6 6 6
Plus généralement, la loi d’une variable aléatoire X prenant uniformément ses
valeurs dans un ensemble fini E est déterminée pour chaque valeur possible x de E
par la formule
P(X = x) = 1/Card(E)
où Card(E) désigne le cardinal de l’ensemble E, c’est à dire le nombre d’éléments dans
l’ensemble E. Dans le jeu de pile ou face, et celui du lancer de dès, l’ensemble des
valeurs possibles a pour cardinal respectivement 2, et 6

E1 = {1, 2} et E2 = {1, 2, . . . , 6} =⇒ Card(E1 ) = 2 et Card(E2 ) = 6

Pour générer sur ordinateur un tel phénomène aléatoire, on utilise à nouveau une seule
variable uniforme U sur [0, 1], et l’on pose
X = 1 + [Card(E) × U ]
Il est donc d’autant plus improbable de choisir un élément donné que le cardinal de
l’ensemble est grand. On a par exemple une chance sur 10.000 de choisir exactement un
élément donné dans un ensemble a 10.000 élements. Dans un jeux de loto où l’on tire
sans remise 6 boules parmi 49, le résultat du tirage X est à valeurs dans un ensemble
ayant approximativement 14 millions d’éléments ! Bien qu’il soit bien plus artistique
de simuler physiquement, ou mécaniquement les tirages successifs des boules dans une
urne, l’on peut numéroter chaque possibilité par une nombre compris entre 1 et 14
millions. Ce jeu de loto simplifié se résume alors à choisir un nombre ente 1 et 14
millions. Pour le simuler sur ordinateur, il suffit de poser
X = 1 + [14 × 106 × U ]
Le joueur gagne s’il devine le nombre X choisi aléatoirement ente 1 et 14 millions par
l’ordinateur.

1.4.2 Les lois d’événements composés

Introduction
Lors d’une expérience aléatoire, un évènement donné, disons A, se réalise avec une
probabilité p connue
P(A) = p ∈ [0, 1]
Par exemple, lors d’un lancer de dés non truqué le chiffre 6 a exactement une chance sur
6 de se réaliser. Si l’on désigne par X le résultat du lancer, la probabilité de réalisation
de l’évènement A = {X = 6} est clairement égale à un sixième
P(A) = P(X = 6) = 1/6
L’expérience se répète dans le temps. Au fil du temps, l’évènement en question a de
fortes chance de se réaliser au moins une fois. Il est par exemple bien improbable que
le chiffre 6 ne soit jamais sorti lors d’un million de lancers de dès. Sur une infinité de
lancers, il est même certain qu’il apparaı̂sse au moins un fois !
Plusieurs questions naturelles se posent alors :
Combien de temps doit on attendre pour que l’évènement qui nous intéresse se
réalise avec une probabilité supérieure à 1/2, voire supérieure à 99/100 ? En reprenant
l’exemple précédent, combien de fois doit on lancer les dés pour que le chiffre 6
apparaisse avec une probabilité supérieure à 9/10 ? Combien d’années doit on jouer
au jeu de loto décrit précédemment, pour avoir une chance de gagner supérieure à 1/2 ?
Dans un autre ordre d’idées, combien de joueurs doivent participer à ce loto, pour avoir
au moins un gagnant avec une probabilité supérieure à 9/10 ?
Associations évènementielles
Avant de répondre précisément aux questions précédentes, il convient de faire un
certain nombre de remarques. Tout d’abord, il est important de souligner que tout
évènement A ayant une probabilité donnée de se réaliser peut s’exprimer
en terme d’une variable uniforme U sur [0, 1]. Plus précisément, on convient que
l’évènement A se réalise lorsque la simulation de la variable uniforme U tombe dans
le segment [0, p]. On met ainsi en correspondance deux évènements ayant la même
probabilité de se réaliser
{“A se réalise”} = {U ∈ [0, p]}
Pour s’en convaincre, il suffit de noter que l’on a
P(A) = p = P(U ∈ [0, p[)
Dans l’exemple du lancer de dés, on peut associer les deux évènement équiprobables
{X = 6} et {U ∈ [0, 1/6]}. Le chiffre 6 a autant de chances de se réaliser que la variable
simulée U a de chances de tomber dans l’intervalle [0, 1/6] de longeur 1/6
{X = 6} = {U ∈ [0, 1/6]} =⇒ P(X = 6) = P(U ∈ [0, 1/6]) = 1/6

Description des modèles

En utilisant ces associations évènementielles, réaliser l’expérience plusieur fois,
revient à simuler une suite de variables aléatoires (Un )n≥1 uniformes et indépendantes
sur [0, 1]. A chaque étape n, on vérifie si l’évènement A s’est réalisé en détectant si la
variable Un est tombée dans l’intervalle [0, p]. Dans une séquence de 4 lancers de dés,
la suite d’évènements
{U1 6∈ [0, 1/6]}, {U2 ∈ [0, 1/6]}, {U3 ∈ [0, 1/6]}, {U4 ∈ [0, 1/6]}
correspond à la situation où le chiffre 6 est sorti à chacun des instants, sauf au premier.
De même, la suite d’évènements
{U1 6∈ [0, 1/6]}, {U2 6∈ [0, 1/6]}, {U3 6∈ [0, 1/6]}, {U4 ∈ [0, 1/6]}
correspond à la situation où le chiffre 6 n’est sorti qu’au dernier instant. En utilisant
les propriétés d’indépendance entre évènements, nous avons
P(“le chiffre 6 ne sort jamais pendant les n premiers lancers”)

= P(U1 6∈ [0, 1/6], U2 6∈ [0, 1/6] , . . . , Un 6∈ [0/1/6])

= P(U1 6∈ [0, 1/6]) × P(U2 6∈ [0, 1/6]) × . . . × P(Un 6∈ [0/1/6])

= (5/6)n
En passant aux évènements complémentaires, la probabilité pour que le chiffre 6
apparaı̂sse au moins une fois est donné par

P(“le chiffre 6 sort au moins une fois lors de n lancers”)

= 1 − P(“le chiffre 6 ne sort jamais pendant n lancers”)

= 1 − (5/6)n

Ces formules élémentaires nous permettent de répondre à toutes nos questions. Avec
plus d’une chance sur deux, le chiffre 6 sort au moins une fois après 4 lancers. En effet,
nous avons

1 − (5/6)n ≥ 1/2 ⇔ (5/6)n ≤ 1/2

⇔ n log (5/6) ≤ log (1/2)
⇔ n log (6/5) ≥ log (2)
⇔ n ≥ log (2)/ log (6/5) ' 3, 80

De même, avec plus de 99 chances sur 100, le chiffre 6 sort au moins une fois après 26
lancers

1 − (5/6)n ≥ 99/100 ⇔ (5/6)n ≤ 1/100

⇔ n log (5/6) ≤ log (1/100)
⇔ n log (6/5) ≥ log (100)
⇔ n ≥ log (100)/ log (6/5) ' 25, 26

Plus généralement, la probabilité pour qu’un événement quelconque A se réalise au

moins une fois, ou jamais lors de n expériences aléatoires indépendantes est donnée par
la formule
P(“A se réalise au moins une fois lors de n expériences”)

= 1 − P(“A ne s’est jamais réalisé après n expériences”)

= 1 − (1 − p)n

L’exemple du loto
Les probabilités de gain au loto sont bien sur plus faibles que celle de voir apparaı̂tre
un 6 lors d’un lancer de dès. Néanmoins les raisonnements sont les mêmes, il suffit de
remplacer la probabilité 1/6 par la probabilité de gagner au loto égale à 1/(14 × 106 ).
On obtient ainsi les formules suivantes
P(“La grille choisie sort au moins une fois lors de n tirages de loto”)

= 1 − P(“La grille choisie ne sort jamais lors de n tirages de loto”)

 n
1
=1− 1− 14×106

On en conclut qu’il nous faut jouer plus de 10 millions de fois pour que notre grille ait
une chance sur deux de sortir au moins une fois
 n  n
1 1
1− 1− ≥ 1/2 ⇔ 1− ≤ 1/2
14 × 106 14 × 106
 
1
⇔ n log 1 − ≤ log (1/2)
14 × 106
14 × 106
 
⇔ n log ≥ log (2)
14 × 106 − 1
14 × 106
 
⇔ n ≥ log 2/ log ' 9, 8 × 106
14 × 106 − 1

En supposant que le joueur avisé joue une seule fois par semaine, il lui faudrait jouer
10 millions de semaines, soit plus de 192 mille ans pour espérer gagner au moins une
fois avec une chance sur deux !
Par des raisonnements analogues, on vérifie qu’il est nécessaire de jouer plus de 65
millions de fois pour que notre grille ait plus de 99 chances sur 100 de sortir au moins
une fois. De façon équivalente, il faudrait à peu près 65 millions de joueurs de loto
acharnés pour être sur que l’un d’entre eux puisse gagner avec une probabilité 99/100
 n  n
1 1
1− 1− ≥ 99/100 ⇔ 1− ≤ 1/100
14 × 106 14 × 106
 
1
⇔ n log 1 − ≤ log (1/100)
14 × 106
14 × 106
 
⇔ n log ≥ log (100)
14 × 106 − 1
14 × 106
 
⇔ n ≥ log (100)/ log ' 65, 7 × 106
14 × 106 − 1

Ces chiffres pour le moins astronomiques n’ont jamais détourné les joueurs de ces jeux
de loto. Ces derniers semblent toujours apprécier l’espoir de gagner, avec de faibles
mises des sommes d’argent parfois inversement proportionnelles aux chances de gain.
1.4.3 Les lois géométriques
Introduction
Les phénomènes aléatoires les plus intéressants se déroulent dans le temps.
Ils se composent d’une successions d’évènements élémentaires aussi complexes
qu’imprévisibles. Pour donner quelques exemples d’actualié, les formations de cyclones,
ou de gaz au niveau atomique, ou encore l’évolution d’actifs financiers sont le fruit de
fluctuations aléatoires météorologiques, quantiques, ou financières se déroulant dans le
temps.
La formalisation mathématique de ces phénomènes aléatoires complexes est fondée
sur des études physiques ou statistiques précises. Les modèles probabilistes formalisés
s’expriment alors sous la forme de processus stochastiques diffusifs ou à sauts. Nous
reviendrons sur ces modèles lorsque nous aborderons la simulation de chaı̂nes de
Markov.
Avant d’étudier en profondeur ces phénomènes aléatoires, il est souvent très utile
de restreindre notre analyse à certaines questions plus élémentaires, tout en simplifiant
à l’extrême ces modèles :
Quelle est la probabilité pour qu’un cyclone se forme à une date donnée ? Quelle
sont nos “chances” pour qu’un appareil électronique tombe en panne après une centaine
d’utilisations ? Quelle est la probabilité pour qu’une action financière s’éffondre à une
date donnée ?
La réponse précise à ces questions nécessite bien évidemment une analyse fine des
processus de turbulence météorologique, des fluctuations thermiques et quantiques des
composants électroniques, ou encore une étude profonde des fluctuations du marché
financier. Néanmoins, une première étude statistique naive permettent d’estimer très
groissièrement les probabilités pour que l’évènement qui nous intéresse se réalise.
Supposons donc qu’un instrument électronique ait une chance sur mille de tomber
en panne à chaque utilisation. Quelle est la probabilité pour qu’il tombe en panne
après une centaine d’utilisations ? Quelle la loi de probabilité complète de ses instants
de pannes ? De même, si un cyclone peut se former chaque jour sur le golfe du Mexique
avec une chance sur mille, quelle est la probabilité d’avoir une catastrophe cyclonique
exactement après six mois ? Enfin, si un actif financier peut s’éffondrer tous les jours
avec un chance sur mille, quelle est la lois exacte des instants de faillites ? Dans tous
les cas, on supposera que les réalisations ou non des évènements en questions à chaque
étape sont mutuellement indépendants.
Nous insiterons sur le fait que les lois de probabilités réelles de ces phénomènes
aléatoires dépendent du paramètre temporel. De plus, dans la plupart des phénomènes
physiques, les états du système sont loins d’être indépendants dans le temps. Les
mécanismes de dépendance sont souvent dictés par des lois fondamentales de la
physique. Par conséquent, les modèles statistiques que nous venons de décrire sont
très peu réalistes, mais ils suffisent pour se donner une première idée du comportement
des systèmes.

Description des modèles

Toutes les questions précédentes correspondent à la donnée d’un événement A
pouvant se réaliser à chaque unité de temps, avec une probabilité fixée

P(A) = p ∈ [0, 1]

On simule la réalisation ou non de cet évènement en convenant que l’évènement A se

réalise, si variable uniforme U générée sur [0, 1] par l’ordinateur tombe dans l’intervalle
[0, p[, de longeur p. Cette stratégie est résumée par la formule

A = {U ∈ [0, p[}

L’évènement contraire Ac , correspondant à la situation où A ne se réalise pas, est

donnée par
Ac = {U 6∈ [0, p[} = {U ∈ [p, 1]}
Dans l’exemple de l’intrument électronique, l’évènement A correspond au
disfonctionnement du système. On convient que la panne est effective lorsque la variable
générée par l’ordinateur tombe dans le segment [0, p[, avec p = 1/1000. Autrement dit,
nous avons

A = {“L’appareil tombe en panne”} = {U ∈ [0, 1/1000[}

et
Ac = {“L’appareil fonctionne”} = {U ∈ [1/1000, 1]}
On vérifie que les probabilités d’occurence de ces évènements correspondent bien aux
données du problème

P(“L’appareil tombe en panne”) = 1/1000

et
P(“L’appareil fonctionne”) = 1 − 1/1000 = 999/1000

La loi géométrique
Pour simuler l’évolution temporelle du phénomène aléatoire, on génère sur
l’ordinateur une suite (Un )n≥1 de variables uniformes sur [0, 1], tout en vérifiant au
fil du temps si l’une des variables tombe dans [0, p[. On note T la première fois où
cela se produit. Cet instant aléatoire T correspond au plus petit entier n, pour lequel
Un tombe dans [0, p[. On note ce temps aléatoire

T = inf {n ≥ 1 : Un ∈ [0, p[}

Plus précisément, si la première simulation U1 tombe dans [0, p[, on convient que
l’évènement A s’est réalisé dès le premier instant. On pose dans ce cas

T =1

Lorsque U1 est dans [p, 1], on simule une nouvelle variable uniforme U2 sur [0, 1]. Si
cette variable tombe dans [0, p[, alors on convient que A s’est réalisé lors de la seconde
expérience, et on pose
T =2
Dans le cas contraire, les variables U1 et U2 sont tombées dans l’intervalle [p, 1]. On
génère alors une nouvelle variable U3 sur [0, 1]. Si cette variable tombe dans [0, p[, alors
on convient que A s’est réalisé après trois expériences, et on pose

T =3

On continue ce procédé jusqu’à ce que l’évènement souhaité se réalise enfin. Ces

premiers instants T de réalisation de l’évènement A sont des variables aléatoires à
valeurs entières. D’après la discussion précédente, nous avons les correspondances
évènementielles suivantes :

{T = 1} = {U1 ∈ [0, p[}

{T = 2} = {U1 6∈ [0, p[, U2 ∈ [0, p[}
{T = 3} = {U1 6∈ [0, p[, U2 6∈ [0, p[, U3 ∈ [0, p[}
..
.
{T = n} = {U1 6∈ [0, p[, U2 6∈ [0, p[, . . . , Un−1 6∈ [0, p[, Un ∈ [0, p[}
..
.

À titre d’exemple, la situation où notre appareil électronique tombe en panne apr̀es
quatre utilisations est décrite par l’évènement

{T = 4} = {U1 6∈ [0, 1/1000[, U2 6∈ [0, 1/1000[, U3 6∈ [0, 1/1000[, U4 ∈ [0, 1/1000[}

D’après les propriétés d’indépendance entre les évènements, la probabilité pour que cet
évènement se réalise est donnée par

P(T = 4) = P(U1 6∈ [0, 1/1000[, U2 6∈ [0, 1/1000[, U3 6∈ [0, 1/1000[, U4 ∈ [0, 1/1000[)
= P(U1 6∈ [0, 1/1000[) × P(U2 6∈ [0, 1/1000[) × P(U3 6∈ [0, 1/1000[)
×P(U4 ∈ [0, 1/1000[)
3
= P(U 6∈ [0, 1/1000[) × P(U ∈ [0, 1/1000[)
 3
1 1
= × 1−
1000 1000
Plus généralement, la probabilité pour qu’il tombe en panne après n utilisations est
donnée par la formule
 n−1
1 1
P(T = n) = × 1−
1000 1000

On a donc approximativement neuf chances sur 10.000 pour que notre appareil tombe
en panne exactement après une centaine d’utilisations
 99
1 1
P(T = 100) = × 1− ' 9 × 10−4
1000 1000

Ces lois des premiers instants où un évènement se produit sont appelées des lois
géométriques. Pour des évènements quelconques A ayant une probabilité P(A) = p
de se réaliser, le premier instant de réalisation

T = inf {n ≥ 1 : Un ∈ [0, p[}

est distribué selon la loi géométrique de paramètre p définie par la formule

∀n ≥ 1, P(T = n) = p (1 − p)n−1 .

1.4.4 Les lois discrètes

Introduction
Il arrive bien souvent que certaines réalisations d’expériences aléatoires soient
bien plus probables que d’autres. Pour illustrer cette remarque, reprenons l’exemple
réchauffé du choix d’une boule dans une urne contenant des boules blanches et noires.
Supposons que l’urne contienne 100 boules noires, et une seule boule blanche. Dans ce
cas, nous avons 99 chances sur 100 de choisir une boule noire, et seulement une chance
sur 100 de trouver la boule blanche.
Il va va de même, lors d’un lancer de dés truqué où les six faces n’ont plus la même
fréquence d’apparition. Nous profiterons de cet exemple du lancer de dès pour souligner
que certains évènements aléatoires sont par essence plus probables que d’autres. Lors
d’un lancer de deux dès non truqués, il est plus probable que la somme des chiffres
obtenus soit égale à 9, plutôt qu’à 10. Pour vérifier cette assertion, il suffit de noter que
l’on peut obtenir 10 de trois façons

5 + 5 = 6 + 4 = 4 + 6 = 10

et le chiffre 9 de 4 façons distinctes

4+5=5+4=6+3=3+6=9
On a donc
P(la somme des deux lancers vaut 10) = 3/36 = 1/12
et
P(la somme des deux lancers vaut 9) = 4/36 = 1/9
Les processus de naissances et morts d’individus dans des environnements plus ou
moins accueillants sont aussi sujet à de telles variations. Par exemple, les naissances
dans des régions désertiques sont moins fréquentes que celles des capitales européennes.
Ces taux de naissances sont souvent reliés à la qualité d’adaptation d’un individu dans
des sites donnés.
Dans un tout autre registre, l’élaboration de portefeuilles financiers est souvent
fondée sur des critères de qualité rendant compte des tendances à la hausse ou à la
baisse d’une collection d’actions. Ainsi, un bon agent financier cherchera à miser le
plus souvent sur des actifs ayant de fort potentiel de gain. Le comportement d’un
agent financier peut ainsi être décrit par des sélections aléatoires d’actifs, certains
plus probables que d’autres en fonction d’une conjoncture financière, ou d’un climat
économique donné.

Description des modèles

D’un point de vue pratique, comment simuler sur un ordinateur de tels phénomènes
aléatoires ? L’idée est relativement simple. Nous sommes en présence d’une variable
aléatoire X pouvant prendre un certain nombre de valeurs, que nous noterons
{x1 , . . . , xd }, avec une probabilité donnée

∀i ∈ {1, . . . , d} P(X = xi ) = pi ∈ [0, 1]

L’interprétation de ces états dépend du problème étudié. Dans l’exemple de l’urne, la

variable X représente la couleur de la boule. Deux seules couleurs sont présentent dans
l’urne, et l’espace d’état est donc réduit à deux points {x1 , x2 }. Dans cette situation,
on peut convenir que x1 représente le choix d’une boule blanche, et x2 celui d’une boule
noire. En reprenant la composition del’urne décrite ci-dessus, nous avons

P(X = x1 ) = 1 − P(X = x2 ) = 1/100

Pour réaliser cette expérience X en terme d’une variable uniforme U sur [0, 1], on
décompose cet intervalle en deux segments de longeurs respectives 1/100 et 99/100

[0, 1] = [0, 1/100]∪]1/100, 1]

On met ensuite en correspondance les évènements équiprobables

{X = x1 } = {U ∈ [0, 1/100]} et {X = x2 } = {U ∈]1/100, 1]}

Autrement dit, l’on convient que l’on choisit une boule blance lorsque le nombre simulé
U tombe dans le sous-intervalle [0, 1/100]. Cette expérience aléatoire X peut s’écrire
sous la forme synthétique suivante

X = x1 1[0,1/100] (U ) + x2 1]1/100,1] (U )

On rapelle que 1I (U ) représente l’indicatrice de l’ensemble I


1 si U ∈ I
1I (U ) =
0 si U 6∈ I

Dans l’exemple du lancer de dès, la variable X représente le résultat de l’expérience.

Dans cette situation, la variable X peut prendre l’une des valeurs

x1 = 1 , x2 = 2 , x3 = 3 , x4 = 4 , x5 = 5 ou x6 = 6

Un exemple de dès truqué est donné par

P(X = 1) = 1/2 , P(X = 2) = 1/4

et

P(X = 3) = 1/8 , P(X = 4) = 1/16 , P(X = 5) = 1/32 et P(X = 6) = 1/32

Dans cet exemple, la chiffre 1 apparait très souvent avec une chance sur 2, alors que
la face 6 est plus improbable, elle se réalise en général une seule fois sur 32. Comme
précédemment, on réalise ce lancer de dès truqué en terme d’une variable uniforme U
sur [0, 1], en décomposant l’intervalle [0, 1] en six segments de longueurs respectives
1/2, 1/4, 1/8, 1/16, 1/32

[0, 1] = I1 ∪ I2 ∪ I3 ∪ I4 ∪ I5 ∪ I6
           
1 1 3 3 7 7 15 15 31 31
= 0, ∪ , ∪ , ∪ , ∪ , ∪ ,1
2 2 4 4 8 8 16 16 32 32
On effectue ensuite les associations évènementielles suivantes

{X = i} = {U ∈ Ii }

pour chaque i ∈ {1, 2, . . . , 6}. Le résultat du lancer de dés X peut s’écrire sous la forme
synthétique suivante

X = 1 × 1I1 (U ) + 2 × 1I2 (U ) + . . . + 6 × 1I6 (U )

Plus généralement, pour simuler une variable X distribuée dans un ensemble fini
{x1 , . . . , xd } avec la loi de probabilité donnée

∀i ∈ {1, . . . , d} P(X = xi ) = pi ∈ [0, 1]

on décompose tout d’abord l’intervalle [0, 1] en d petits segments de longueurs
respectives p1 , p2 , . . . , pd

[0, 1] = I1 ∪ I2 ∪ I3 ∪ . . . ∪ Id
= [0, p1 ] ∪ ]p1 , p1 + p2 ] ∪ ]p1 + p2 , p1 + p2 + p3 ] ∪ . . . ∪ ]p1 + . . . + pd−1 , 1]

Pour chaque indice i ∈ {1, . . . , d}, on fait les associations évènementielles suivantes :

{X = xi } = {U ∈ Ii }

Autrement dit, on convient que la variable aléatoire X prend la valeur xi , lorsque le

nombre simulé U dans l’intervalle [0, 1] prend une valeur dans l’intervalle Ii . En termes
plus synthétiques, la variable X peut s’écrire sous la forme suivante
d
X
X= xi 1Ii (U )
i=1

1.4.5 Les lois continues

Introduction
Dans le chapitre précédent nous avons présenté des techniques de simulation
d’expériences aléatoires prenant un nombre de valeurs au plus dénombrable. Ces
expériences sont parfois l’expression simplifiée de phénomènes aléatoires plus complexes,
à valeurs dans des espaces continus, ou de dimensions multiples. C’est le cas des
arrivées de clients dans des files d’attente, les instants de désintégration de particules
quantiques, la répartition de la chaleur, les impacts de météorites sur le globe terrestre,
les fluctuations de prix d’actifs financiers dans des marchés boursiers, les explorations
aléatoires de régions par des populations d’individus, etc.
Tous ces aléas sont de nature très variée : répartitions gaussiennes, inter temps
exponentiels, statistiques d’ordre uniformes, distributions log-normales, densités de
type Gamma, lois de Cauchy, de Rayleigh, de Pareto, de Weibull, distributions
uniformes sur des surfaces, ou des volumes quelconques, etc.
Chacune de ses lois traduit une répartition statistique particulière des aléas d’une
expérience. La loi gaussienne peut s’interpréter comme la répartition aléatoire de la
chaleur dans une droite. Cette distribution en forme de cloche souligne les différents
degrés de chaleur en chaque site. Les régions éloignées du point de chauffage initial sont
plus froides. Les températures autour de la source de chaleur se répartissent en forme
de cloche. Les distributions uniformes sur des sphères peuvent s’interpréter comme des
lieux d’impact de météorites sur une planète. Ces lois peuvent être étonnament reliées à
des projections de variables aléatoires gaussiennes sur chaque axe du repère considéré !
Les lois exponentielles peuvent s’interpréter comme des horloges aléatoires représentant
des instants d’arrivée d’individus dans des files d’attente, des durées de traitement
d’information dans des réseaux de télécommunication, des durées de fonctionnement
de machines, des désintégrations de particules, etc.
Comment simuler de tels phénomènes aléatoires sur un ordinateur ? Ces aléas
continus peuvent-ils toujours s’exprimer en terme de variables uniformes sur l’intervalle
[0, 1] ? Existe-t-il des procédés de simulation de lois universels ? Avant de répondre
à ces questions, examinons tout d’abord l’expression mathématique des phénomènes
continus. Pour des expériences aléatoires réelles X, telles les répartition de température
sur une barre rectiligne, les instants d’arrivées, de pannes, ou de ruptures, la probabilité
pour que le résultat aléatoire soit compris entre deux valeurs a, et b, correspond à une
aire entre le graphe d’une fonction positive

p : x ∈ R −→ p(x) ∈ [0, ∞[

et le segment [a, b].

p(x)

a b x

Fig. 1.1 – Intégrale entre a et b

Cette aire est appelé l’intégrale de la densité de probabilité p entre le point a et

le point b, et on la note
Z b Z
P(a ≤ X ≤ b) = p(x) dx = p(x) dx
a [a,b]

Cette interprétation probabiliste offre une interprétation physique à la théorie de

l’intégration. Par exemple, la somme des aires sur deux intervalles disjoints [a1 , b1 ], et
[a2 .b2 ] correspond à la probabilité pour que le résultat de l’expérience X prenne ses
valeurs dans le premier intervalle [a1 , b1 ], ou dans le second [a2 , b2 ]
Z
P(X ∈ [a1 , b1 ] ou X ∈ [a2 , b2 ]) = p(x) dx
[a1 ,b1 ]∪[a2 ,b2 ]
Z b1 Z b2
= p(x) dx + p(x) dx
a1 a2
= P(X ∈ [a1 , b1 ]) + P(X ∈ [a2 , b2 ])

Horloges exponentielles
Comme nous l’avons souligné dans l’introduction, l’évolution aléatoire du nombre de
clients dans une file d’attente est clairement le fruit d’un processus d’arrivées aléatoires
d’individus au fil du temps. Ces dates d’arrivées peuvent être représentées par des
variables aléatoires T à valeurs dans la demi-droite [0, ∞[. Leur nature est souvent
exponentielle. Plus précisement, la probabilité pour qu’une personne arrive dans une
file d’attente à un instant T compris entre t1 et t2 est donnée par l’intégrale d’une
fonction exponentielle entre les dates en question
Z t2
P(t1 ≤ T ≤ t2 ) = λ e−λs ds = e−λt1 − e−λt2
t1

En particulier, nous avons

P(T ≥ t) = lim P(t ≤ T ≤ s) = e−λt

s↑∞

On utilise souvent la notation synthétique suivante

P(T ∈ dt) =déf. P(T ∈ [t, t + dt[) = λ e−λt dt

D’aprés la remarque précédente, la probabilité pour qu’une personne arrive entre deux
instants très rapprochés t et t + dt, sachant qu’elle n’est pas encore arrivée avant t, est
égale au produit λdt

P(T ∈ [t, t + dt[ et T ≥ t)

P(T ∈ dt | T ≥ t) =déf. P(T ∈ [t, t + dt[ | T ≥ t) =
P(T ≥ t)
P(T ∈ [t, t + dt[)
= = λ dt
P(T ≥ t)

Le paramètre λ > 0 est appelé l’intensité de la variable aléatoire. Notons que pour
t1 = 0, et lorsque t2 ↑ ∞, nous avons

P(T < ∞) = 1 − lim e−λt2 = 1

t2 →∞
 Les amateurs de l’intégration par parties, pourrons noter que l’instant moyen
d’arrivée d’un client dans une file d’attente est inversement proportionnel à l’intensité
Z ∞ Z ∞ Z ∞
−λs ∂ −λs
E(T ) = s × P(T ∈ ds) = sλe ds = − s (e ) ds
0 0 0 ∂s
1 ∞
Z
1
= −[s e−λs ]∞0 + λ e−λs ds =
λ 0 λ

Le paramètre λ correspond donc à la fréquence d’arrivée des clients par unité de temps,
disons par minute. Lorsque λ = 10 clients par minute, l’instant d’arrivée moyen d’un
client est de 1/10 de minute. Autrement dit, les clients arrivent en moyenne toutes les
six minutes dans la file d’attente.
La simulation de ces instants d’arrivée aléatoire est à nouveau l’expression d’une
variable uniforme U sur [0, 1] ! Pour être plus précis, commençons par poser
1
T =− log U
λ
On vérifie sans trop de peine que l’on a

{t1 ≤ T ≤ t2 } = {e−λt2 ≤ U ≤ e−λt1 }

Cette association événementielle nous conduit au résultat recherché

P(t1 ≤ T ≤ t2 ) = P(e−λt2 ≤ U ≤ e−λt1 ) = e−λt1 − e−λt2

Pour simuler l’arrivée d’un client dans une file d’attente, on prend donc tout simplement
le logarithme d’un nombre uniforme sur [0, 1] fourni par l’ordinateur, et l’on multiplie
le résultat par −1/λ.
Chapitre 2

Quelques outils de simulation

2.1 La méthode d’inversion

2.1.1 Le principe de base
Derrière la méthode de simulation d’horloges exponentielles décrite dans la
section précédente, se cache une technique universelle. Pour être plus précis, nous
commencerons par noter que la fonction de répartition de la loi exponentielle de
paramètre λ est donnée par la formule explicite
Z x
F (x) = P(T ≤ x) = λ e−λt dt
0
Z x
∂  −λt 
= − e dt = −[e−λt ]x0 = 1 − e−λx
0 ∂t

L’inverse de cette fonction est facilement calculable en chaque point u ∈]0, 1[

F (x) = u ⇐⇒ e−λx = (1 − u)
1
⇐⇒ x = F −1 (u) = − log (1 − u)
λ
Cet inverse ressemble étrangement à l’expression de la variable exponentielle en terme
d’une variable uniforme présentée dans la section précédente.
1
T = − log U
λ
En effet, les variables U et (1 − U ) sont toutes deux uniformes sur l’intervalle [0, 1] ! Il
semble donc qu’en posant
T = F −1 (U )
l’on obtienne un procédé de simulation universel d’une variable aléatoire ayant pour
fonction de répartition F . Pour vérifier cette intuition, notons F la fonction de

27
répartition d’une variable réelle quelconque X définie par

F (x) = P (X ≤ x) .

Lorsque cette fonction est inversible, on a clairement pour chaque x ∈ R

Z F (x)
−1


P F (U ) ≤ x = P(U ≤ F (x)) = du = F (x).
0

Il en découle que la variable F −1 (U ) est bien distribuée selon la même loi que la variable
aléatoire X.
Inversement, si X est une variable aléatoire réelle ayant une fonction de répartition
strictement croissante F : R → [0, 1] alors F (X) est une variable uniforme sur [0, 1].
On vérifie aisément ce résultat en notant que pour chaque u ∈ [0, 1], nous avons

P(F (X) ≤ u) = P(X ≤ F −1 (u)) = F ◦ F −1 (u) = u

Cette propriété est souvent utilisée en pratique non pas pour simuler des variables
uniformes mais plutôt pour vérifier si un échantillon de données statistiques
indépendantes suit bien une loi donnée !
Illustrons graphiquement cette technique d’inversion. La nature de la fonction de
répartition décrit les accroissement de masse de probabilité entre les points. Pour des
variables aléatoires réelles admettant une densité de probabilité p, la accroissement de
la fonction de répartition
Z
F (x) = P(X ≤ x) = p(y) dy
]−∞,x]

entre les points x1 et x2 correspondent à l’intégrale de la densité de probabilité sur le

segment [x1 , x2 ]
Z x2
F (x2 ) − F (x1 ) = P(X ∈ [x1 , x2 ]) = p(x) dx
x1

La fonction de répartition d’une densité de probabilité p ayant deux modes séparés

présentera deux paliers de croissance sur chacun d’entre eux. De plus les variations
d’aires sont toujours plus faibles sur les bords de l’axe des réels. Ceci se traduit par des
fonctions de répartition F (x) convergeant vers 0, lorsque x ↓ −∞, et vers 1, lorsque
x ↑ +∞.

La figure précédente illustre bien la répartition sur les masses de probabilités d’une
série de points choisis au hasard par le méthode d’inversion

X1 = F −1 (U1 ) , X2 = F −1 (U2 ) , . . . , Xn = F −1 (Un )

 p(x)

x
X1,...,Xn,...

F9x)

U1,...,Un,...

Fig. 2.1 – Méthode d’inversion

où U1 , U2 , . . . , Un désignent une suite de variables uniformes sur [0, 1].

Dans les sections suivantes, nous illustrerons la performance de méthode d’inversion
de la fonction de répartition avec les algorithmes de simulation des lois géométriques,
de Cauchy, de Rayleigh, et de Weibull. Nous essayerons dans chaque situation d’offrir
une preuve mathématique rigoureuse des formules d’inversion associée à ces diverses
lois. Le lecteur ayant aucune connaissance sur les techniques de changement de variable
dans les intégrales à de fortes chances de s’ennuyer, voire de s’ennerver ! Afin d’éviter
de telles crises, notre recommendation dans ce cas est d’éviter tout simplement ces
démonstrations, et de passer directement à la description des algorithmes de simulation
de ces lois.

2.1.2 Retour aux lois géométriques

Les lois géométriques correspondent au nombre d’expériences nécessaires pour faire
apparaı̂tre un événement aléatoire donné. Nous avons définie et étudiée la simulation
de ces lois dans la section 1.4.3.
Ces distributions temporelles ou de longeurs discrètes sont des modèles simplifiés de
phénomènes aléatoires naturels souvent plus complexes. En biologie, ces lois sont parfois
utilisés pour simuler des phénomènes de morphogenèse dans la croissance des conifères,
des longueurs d’exons de gènes. En physique, la loi géométrique est plutôt utilisée
pour décrire grossièrement des instants de désintégration de particules radioactives,
ou d’autre type d’absorption de particules dans des puits de potentiels. En économie,
ces lois peuvent représenter des dates de hausse ou de chute des valeurs de produits
financiers.
Cette distribution géométrique est l’analogue de la loi exponentielle de paramètre
λ = − log (1 − p).
En effet, si X est une variable exponentielle de paramètre λ alors, pour tout n ≥ 1,
nous avons
Z n  
P(n − 1 ≤ X < n) = λ e−λ x dx = e−(n−1)λ 1 − e−λ = p (1 − p)n−1 .
n−1

Par conséquent la variable Y = 1 + [X] est distribuée selon une loi géométrique de
paramètre p. En utilisant les techniques de simulation de lois exponentielles décrites
précédemment, on vérifie sans trop de peine que la variable
Y = 1 + [log (U )/log (1 − p)]
avec U uniforme sur [0, 1], suit encore une loi géométrique de paramètre p. On pourra
noter que ce procédé de simulation est bien plus rapide que celui présenté dans la
section 1.4.3. Une seule variable uniforme permet de décrire l’apparition d’un événement
aléatoire de probabilité donnés, sans simuler les échecs !

2.1.3 Les lois de Cauchy

La loi de Cauchy est la hantise des probabilistes. Cette distribution de probabilités
particulière n’admet aucune moyenne, et a fortiori aucun moments ! Il est donc
impossible d’effectuer des calculs élémentaires, et y développer une quelconque
intuition.
Cette distribution d’aléas est souvent associée à des phénomènes aléatoires très
oscillants. En optique, cette loi de probabilité est parfois utilisée pour modéliser les
variations de l’indice de réfraction d’un milieu liées à la fréquence du rayonnement
monochromatique passant de l’air à ce milieu. En économie, ces lois représentent plutôt
les fluctuations d’amplitudes des logarithmes des rapports journaliers de prix d’actifs
financiers.
Pour être plus précis, une variable de Cauchy X de paramètre σ > 0, est distribuée
sur tout segment [a, b] de la droite réelle selon la loi
Z b
σ
P(X ∈ [a, b]) = 2 2
dx.
a π (x + σ )

La fonction de répartition F de X est définie pour chaque x ∈ R par la formule

Z x
1 x
Z
σ dy
F (x) = 2 2
dy =
−∞ π(y + σ ) σ −∞ π((y/σ)2 + 1)
En effectuant le changement de variable

z = y/σ ∈] − ∞, x/σ] =⇒ dz = dy/σ

la fonction F est définie plus simplement par la formule

Z x/σ
dz
F (x) = 2 + 1)
−∞ π(z
Il reste à poser
sin θ
z = tan θ =déf. avec θ ∈] − π/2, arctan (x/σ)]
cos θ
Dan ce cas, nous avons
∂z cos2 θ + sin2 θ dz
dθ = 1 + tan2 θ dθ ⇒


dz = dθ = 2
= dθ
∂θ cos θ 1 + z2
et par conséquent
Z arctan (x/σ)
1 1 1 x
F (x) = dθ = + arctan
π − π2 2 π σ

Cette fonction de répartition est inversible, et son inverse F −1 est donné pour chaque
u ∈ [0, 1] par la formule
 
−1 1
F (u) = σ tan π (u − ) .
2
La simulation d’une loi de Cauchy est donc claire, on simule tout d’abord une
variable uniforme U sur l’intervalle ]0, 1[, et on pose
 
−1 1
X = F (U ) = σ tan π (U − )
2
D’après ce qui précède, cette variable est bien une variable aléatoire de Cauchy de
paramètre σ > 0. On pourra noter que la variable aléatoire angulaire
1
V = π (U − )
2
est uniforme sur l’intevalle ] − π/2, π/2[. L’algorithme de simulation peut donc
s’interpréter comme la tangente d’un angle choisi au hasard dans l’intervalle ]−π/2, π/2[

X = σ tan V

Le paramètre de dilatation σ vient augmenter ou diminuer l’amplitude de cette

tangente.
2.1.4 Les lois de Rayleigh-Weibull
Les lois de Rayleigh-Weibull sont des distributions de probabilités sur la demi-droite
des réels positifs. Tout comme les lois exponentielles, ces distributions permettent de
décrire des instants, des durées, des longueurs, ou des tailles de files d’attentes aléatoires.
Cependant, les lois de Rayleigh sont plutôt associées à des phénomènes aléatoires
temporels bien plus rapides, ou à des longueurs aléatoires de plus faibles amplitudes.
Elle sont ainsi utilisées pour décrire des bruits de mesure dans des recepteurs de
transmission, des instants de dégradation de matériel, des durées d’appels dans certains
réseaux de télécommunication, des temps d’accés sur des sites internet, etc.
Une variable aléatoire X de Rayleigh de paramètre σ > 0 est distribuée sur la demi
droite [0, ∞[ selon la loi

x − x22
P(X ∈ dx) = e 2σ 1[0,∞[ (x) dx.
σ2
La fonction de répartition F de X est définie pour chaque x ∈ [0, ∞[ par la formule
explicite suivante
Z x
x2
   
∂ y2
F (x) = P(X ≤ x) = − e− 2σ2 dy = 1 − exp − 2
0 ∂y 2σ

Cette fonction est clairement inversible, et son inverse est donné par la fonction

u ∈ [0, 1[−→ F −1 (u) = σ − log (1 − u) ∈ [0, ∞[

p

La simulation d’une variable de Rayleigh de paramètre σ > 0 consiste donc à prendre

la racine du logarithme d’une variable U uniforme sur [0, 1]
p
X = σ − log U

Les lois de Weibull sont de simples extensions des lois de Raighley. Plus précisément,
une variable aléatoire X de Weibull de paramètres a, b > 0, est distribuée selon la loi
a a−1 −(x/b)a
P(X ∈ dx) = x e 1[0,∞[ (x) dx.
ba
Avec ces notations, une loi de Rayleigh de paramètre σ est donc une loi de Weibull de
paramètres
a = 2 et b = σ
Par des raisonnements analogues aux précédents, la fonction de répartition F de X est
définie explictement pour chaque x ∈ [0, ∞[ par la formule suivante
Z x
∂ a a
F (x) = (−e−(z/b) ) dz = 1 − e−(x/b)
0 ∂z
Cette fonction est clairement inversible et son inverse F −1 est donné pour chaque
u ∈ [0, 1[ par
F −1 (u) = b (− log (1 − u))1/a .
Si U est une variable aléatoire uniforme sur ]0, 1] alors

X = b × (− log U )1/a

est une variable aléatoire de Weibull de paramètres a, b > 0. Pour conclure, on pourra
noter que les variables exponentielles de paramètre λ = 1 données par

Y = − log U

et les variables de Weibull de paramètres (a, b) sont reliées par la formule suivante

X = b × Y 1/a

En terme de fiabilité, la variable exponentielle de référence correspond à des instants

de pannes purement accidentelles. Les paramètres (a, b) des variables de Weibull
représentent les degrés d’usure d’un appareil. Par exemple, les durées de fonctionnement
d’un appareil plus usé sont d’autant plus courtes que le paramètre b est faible, ou a est
grand.
Tout une panoplie de variables réelles peuvent être simulées en utilisant ce principe
d’inversion de la fonction de répartition. Son application nécessite néanmoins d’avoir
une description explicite de l’inverse de cette fonction. Pour plus d’information, nous
renvoyons le lecteur à l’ouvrage [?].

2.2 Changements de variables

2.2.1 Introduction
Il existe une variété d’expériences aléatoires X prenant leurs valeurs dans des espaces
multi-dimensionnels Rd , avec d ≥ 1. C’est le cas des répartitions de température sur
une surface ou dans un volume, les évolutions de photons dans des tissus lors d’examens
d’imagerie médicale, les perturbations de mesures sur chaque composante d’un capteur
éléctronique, les impacts de météorites sur une surface, etc. Dans la plupart de ces
situations, la probabilité pour que le résultat aléatoire soit compris dans un pavé A ⊂ Rd
correspond à un volume entre le pavé en question, et le graphe d’une fonction positive

p : x ∈ Rd −→ p(x) ∈ [0, ∞[

Cette aire est appelé l’intégrale de la densité de probabilité p sur la pavé A, et on la

note Z
P(X ∈ A) = p(x) dx
A
2.2.2 Lois uniformes sur des cubes
Comment choisir un mot au hasard dans un texte, ou dans un manuscript de 200
pages d’épaisseur ? Existe-t-il un procédé de simulation permettant de colorier un point
choisi au hasard dans un tableau ? Les pixels d’une image numérique peuvent-ils être
allumés au hasard, et indépendemment les uns des autres ? Peut-on faire apparaitre un
personnage de jeu vidéo dans une pièce virtuelle, de façon totalement imprévisible ? Un
adepte du pointillisme peut-il peindre par petites touches totallement hasardeuses un
tableau ? Peut-on choisir au hasard un point dans un glaçon cubique ?
D’un point de vue probabiliste, toutes ces petites questions reviennent à trouver
un algorithme de simulation d’une variable aléatoire uniforme sur des rectangles, des
cubes, ou tout autre pavé d’un espace euclidien Rd , de dimension d ≥ 1.
Nous commencerons par examiner le cas des pavés dans le plan. La densité de
probabilité d’une variable aléatoire uniforme X = (X1 , X2 ) sur un rectangle quelconque

A = ([a1 , b1 ] × [a2 , b2 ])

est donnée par la fonction

1
p : x ∈ C 7→ p(x) = 1A (x)
Aire(A)
La probabilité pour que le point aléatoirement choisi X tombe dans un sous-ensemble
régulier B ⊂ A est simplement donnée par le rapport des aires
Aire(B)
P(X ∈ B) =
Aire(A)
La fonction p étant une fonction indicatrice sur le carrée, de hauteur 1/Aire(A), on
retrouve le fait que la probabilité de l’évènement {X ∈ B} est le volume du pavé de
base B et de hauteur 1/Aire(A).
Il est très aisé de simuler le point X = (X1 , X2 ) à l’aide d’un couple de variables
uniformes et indépendants (U1 , U2 ) sur le segment [0, 1]. Il suffit de poser

X1 = a1 + (b1 − a1 ) U1 et X2 = a2 + (b2 − a2 ) U2

Pour vérifier cette assertion, on considère la fonction θ qui à un point (u1 , u2 ) ∈ [0, 1]2
associe le point

(x1 , x2 ) = θ(u1 , u2 ) = (a1 + (b1 − a1 )u1 , a2 + (b2 − a2 )u2 ) ∈ A = ([a1 , b1 ] × [a2 , b2 ])

Cette tranformation envoie le carré unité [0, 1]2 sur le carré

A = θ([0, 1]2 )

Les points (U1 , U2 ) choisis sur [0, 1]2 sont ainsi envoyés sur des points (X1 , X2 ) =
θ(U1 , U2 ) du carré A. Mieux encore, la fonction θ est une bijection entre ces deux
ensembles. Autrement dit, θ est une une correspondance bi-univoque entre [0, 1]2 et A.
L’inverse de θ est donné par la formule suivante
 
−1 x1 − a1 x2 − a2
(u1 , u2 ) = θ (x1 , x2 ) = ,
b1 − a1 b2 − a2

Ceci nous permet simplement de vérifier que l’on a

P(X ∈ B) = P(θ(U1 , U2 ) ∈ B) = P((U1 , U2 ) ∈ θ−1 (B)) = Aire(θ−1 (B))

Il nous reste cependant à vérifier que l’on a

Aire(B)
Aire(θ−1 (B)) =
Aire(A)

Avant de vérifier cette assertion, nous commencerons par noter que la transformation
θ est loin de préserver les aires des carrés. Par exemple, le carré unité [0, 1]2 d’aire 1
est envoyé sur un carré

A = θ([0, 1]2 ) = ([a1 , b1 ] × [a2 , b2 ])

dont l’aire est donnée par la formule

Aire(A) = (b1 − a1 )(b2 − a2 )

Pour être plus précis, nous avons les associations d’aires infinitésimales suivantes

∂u1 ∂u2 ∂u1 ∂u2

du1 du2 − × dx1 dx2 (2.1)
∂x1 ∂x2 ∂x2 ∂x1

Pour plus détails concernant ces associations d’aires, nous renvoyons le lecteur à la
section 2.2.3 consacrée au variations d’aires engendrées par des transformations planes.
Dans notre cas, nous avons
∂u1 1 ∂u2 1 ∂u1 ∂u2
= , = , et =0=
∂x1 (b1 − a1 ) ∂x2 (b2 − a2 ) ∂x2 ∂x1

Ceci nous conduit aux associations d’aires infinitésimales

1
du1 du2 × dx1 dx2
Aire(A)

On obtient finallement le résultat recherché

Aire(θ−1 (B)) = P((U1 , U2 ) ∈ θ−1 (B))

Z Z
1 Aire(B)
= du1 du2 = × dx1 dx2 =
θ−1 (B) Aire(A) B Aire(A)
 La simulation d’une variable uniforme sur le cube unité [0, 1]d consiste à prendre
un vecteur
(U1 , . . . , Ud )
formé de d variables aléatoires indépendantes (Uk )k=1,...,d de loi uniforme sur [0, 1]. Par
exemple, choisir un point au hasard (U1 , U2 , U3 ) dans un cube [0, 1]3 revient à choisir
chacune des coordonnées Uk uniformément sur chacune des trois arêtes de base [0, 1].
De même, pour choisir un point uniformément distribué sur un pavé de la forme

[a1 , b1 ] × · · · × [an , bn ]

il suffit de poser

(X1 , . . . , Xn ) = (a1 + (b1 − a1 )U1 , . . . , an + (bn − an )Un )

Pour vérifier que chacune des composantes Xk est uniformément distribuée sur chaque
arête [ak , bk ], on note simplement que pour tout segment [αk , βk ] ⊂ [ak , bk ], nous avons
  
αk − ak βk − ak βk − αk
P (Xk ∈ [αk , βk ]) = P Uk ∈ , =
bk − ak bk − ak bk − ak

2.2.3 Lois uniformes sur des surfaces

Les régles de l’intégration traduisent les aires globales de surfaces A ⊂ R2 en termes
d’aires infinitésimales Z
Aire(A) = du1 du2
A
La mesure du1 du2 désigne la mesure de Lebesgue sur le plan, cette mesure s’exprime
simplement par le fait que l’aire d’un rectangle est donnée par la formule classique
apprise au lycée
Z
Aire([a, b] × [c, d]) = du1 du2 = (b − a) × (d − c)
[a,b]×[c,d]

On doit interpréter du1 du2 comme une aire virtuelle d’un petit carré infinitésimal dont
les cotés sont de longueurs respectives du1 et du2 .

L’intégrale consiste simplement à sommer les aires tous les petits carrés
infinitésimaux contenus dans la surface en question. Considérons désormais une
transformation θ suffisament régulière du plan dans lui même. Cette transformation
peut s’interpréter comme une déformation physique et élastique d’une surface plane.
Chaque point (u1 , u2 ) est envoyé sur un nouveau point

(x1 , x2 ) = θ(u1 , u2 )
 du1 du2

du2

du1

Fig. 2.2 – Aire de carrés infinitésimaux

Les petits carrés infinitésimaux

[u1 , u1 + du1 ] × [u2 , u2 + du2 ]

de surface du1 du2 autour de points (u1 , u2 ) sont alors transformés en petits carrés

[x1 , x1 + dx1 ] × [x2 , x2 + dx2 ]

de surface dx1 dx2 , autours des points transformés (x1 , x2 ) = θ(u1 , u2 ). Les variations
d’étirement du plan subies par l’action de la transformation élastique θ se traduisent
par une variation des longueurs des cotés des carrés. Plus précisément, dans la base des
longueurs de référence (du1 , du2 ) sur chaque coordonnées, nous avons

∂x1 ∂x1 ∂x2 ∂x2

dx1 = du1 + du2 et dx2 = du1 + du2
∂u1 ∂u2 ∂u1 ∂u2
L’objectif est de calculer l’aire des rectangles

[x1 , x1 + dx1 ] × [x2 , x2 + dx2 ]

en fonction des variations de longueurs infinitésimales

∂xi
, avec i, j ∈ {1, 2}
∂uj

Autrement dit, il convient de calculer l’aire du parallélogramme engendré par les

vecteurs ! !
∂x1 ∂x2
dx1 = ∂u1 et dx2 = ∂u1
∂x1 ∂x2
∂u2 ∂u2

dans la base de référence (du1 , du2 ).

 dx2/du2
dx1 dx2
dx2
dx1/dx2

du2 dx1

du1
dx2/du1 dx1/du1

Fig. 2.3 – Aire infinitésimale dx1 dx2

\
Si (dx 1 , dx2 ) désigne l’angle formé par ces deux vecteurs, nous avons

\ h
sin (dx 1 , dx2 ) = r 2  2
∂x2 ∂x2
∂u1 + ∂u2

où h désigne la hauteur du parallélogramme. Il est aisé de vérifier que l’aire de ce

parallélogramme est donnée par la formule produit classique
s
∂x1 2 ∂x1 2
  
h × (base dx1 ) = h × +
∂u1 ∂u2
s 2  2 s
∂x2 2 ∂x2 2
  
∂x1 ∂x1 \
= + × + × sin (dx 1 , dx2 )
∂u1 ∂u2 ∂u1 ∂u2
 \
Il nous reste donc à évaluer le sinus de l’angle (dx 1 , dx2 ). Pour simplifier nos calcul, on
considérer les vecteurs normalisés sur le cercle unité
dx1 dx2
dy1 = r 2  2 et dy2 = r 2  2
∂x1 ∂x1 ∂x2 ∂x2
∂u1 + ∂u2 ∂u1 + ∂u2

On remarque que l’on a

\
(dx \
1 , dx2 ) = (dy1 , dy2 )

En identifiant le cercle unité à l’ensemble des nombres complexes de module 1, nous

avons
∂x1 ∂x1
∂u1 ∂u2
dy1 = eia1 = cos a1 + i sin a1 = r 2  2 + i r  2  2
∂x1 ∂x1 ∂x1 ∂x1
∂u1 + ∂u2 ∂u1 + ∂u2
∂x2 ∂x2
∂u1 ∂u2
dy2 = eia2 = cos a2 + i sin a2 = r 2  2 + i r  2  2
∂x2 ∂x2 ∂x2 ∂x2
∂u1 + ∂u2 ∂u1 + ∂u2

où a1 , et a2 désignent respectivement les angles entre l’axe des abscisses et le vecteurs
dy1 , et dy2 .

Avec ces notations, nous avons

\
(dx \
1 , dx2 ) = (dy1 , dy2 ) = a2 − a1

Il reste alors à remarquer que l’on a

cos (a2 − a1 ) + i sin (a2 − a1 ) = ei(a2 −a1 ) = eia2 e−ia1 = dy2 × dy1 (2.2)

où dy1 désigne le conjugué de dy1

dy1 = e−ia1 = cos a1 − i sin a1

En identifiant les parties complexes dans le produit (2.2), on en conclut que

sin (a2 − a1 ) = sin a1 cos a2 − cos a1 sin a2

et finallement
 
\ 1 1 ∂x1 ∂x2 ∂x1 ∂x2
sin (dx 1 , dx2 ) = r 2 r −
∂x1
2 
∂x1 ∂x2
2 
∂x2
2 ∂u2 ∂u1 ∂u1 ∂u2
∂u1 + ∂u2 ∂u1 + ∂u2
 1 dx2

dy2

dx1
dy1
a2
a1

−1 1

−1

Fig. 2.4 – Angles a1 , a2 sur le cercle trigonométrique

L’aire recherchée est donc donnée par la formule

 
∂x1 ∂x2 ∂x1 ∂x2
h × (base dx1 ) = −
∂u2 ∂u1 ∂u1 ∂u2

En conclusion, nous avons les associations d’aires infinitésimales

 
∂x1 ∂x2 ∂x1 ∂x2
dx1 dx2 − du1 du2
∂u2 ∂u1 ∂u1 ∂u2

Le terme du1 du2 représente l’unité d’aire infinitésimale dans laquelle nous avons effectué
nos calculs. Ces formules témoignent du fait que l’aire d’une surface D suffisament
régulière Z
Aire(D) = du1 du2
D
est transformée par θ en une surface θ(D) d’aire
Z Z
∂x1 ∂x2 ∂x1 ∂x2
Aire(θ(D)) = dx1 dx2 = − du1 du2
θ(D) D ∂u2 ∂u1 ∂u1 ∂u2
Plus généralement, pour des fonctions intégrables sur la surface θ(D), nous avons la
formule de changement de variables
∂θ1 ∂θ2 ∂θ1 ∂θ2
Z Z
f (x1 , x2 ) dx1 dx2 = f (θ(u1 , u2 )) − du1 du2
θ(D) D ∂u2 ∂u1 ∂u1 ∂u2

avec la transformation

(u1 , u2 ) ∈ D 7→ θ(u1 , u2 ) = (θ1 (u1 , u2 ), θ2 (u1 , u2 )) = (x1 , x2 ) ∈ θ(D)

2.2.4 Lois uniformes sur des disques

Comment placer au hasard le diament d’une platine sur un disque vinyl ? Peut-on
élaborer une pizza en plaçant tous les ingrédients de façon totalement hasardeuse ?
Comment simuler un mauvais joueur de flechette sur ordinateur ? Un géomètre a-t-il la
possibilité de choisir un point X = (X1 , X2 ) au hasard dans le disque unité

D = (x1 , x2 ) ∈ R2 ; x21 + x22 < 1 ?



Comme dans le cas des choix de points au hasard dans des pavés examiné dans la
section précédente, nous commencerons par noter que la probabilité de choisir (X1 , X2 )
dans un sous ensemble régulier B ⊂ D correspond simplement au rapport des surfaces
de B et de D
Z
Aire(B)
P ((X1 , X2 ) ∈ B) = avec Aire(D) = dx1 dx2 = π
Aire(D) D

Tout point (x1 , x2 ) du disque peut être paramétrisé par la donnée d’un angle a et d’un
rayon r. Plus précisément, nous avons la formule de changement de coordonnée polaires

x1 = r cos a
x2 = r sin a

L’application

θ : D0 = D − ([0, 1[×{0}) −→ ]0, 1[×]0, 2π[

(x1 , x2 ) 7→ θ(x1 , x2 ) = (r, a)

transforme donc de façon bi-univoque le disque unité D, auquel on a oté soigneusement

le petit segment ([0, 1[×{0}), en un joli petit rectangle

θ(D0 ) =]0, 1[×]0, 2π[

Si (X1 , X2 ) est uniforme sur D, quelle est la loi du point

θ(X1 , X2 ) = (R, A)
 x2

r sin(b)
r
a

−1 0 r cos(a) 1 x1

Fig. 2.5 – Coordonnées polaires

sur le pavé ]0, 1[×]0, 2π[ ? Pour répondre à cette question, il convient de noter les
associations d’éléments de surface infinitésimale
∂x1 ∂x2 ∂x1 ∂x2
dx1 dx2 − × drda
∂r ∂a ∂a ∂r
Dans notre situation, nous avons
∂x1 ∂x2
= cos a , = sin a
∂r ∂r
et
∂x1 ∂x2
= −r sin a , = r cos a
∂a ∂a
Par conséquent, nous avons
∂x1 ∂x2 ∂x1 ∂x2
= r cos2 a + sin2 a = r


−
∂r ∂a ∂a ∂r
et donc
dx1 dx2 r drda
La densité de probabilité du point (X1 , X2 )
1
pX1 ,X2 (x1 , x2 ) = 1D (x1 , x2 ) dx1 , dx2
π
est donc transformée par θ en la densité de probabilité du point θ(X1 , X2 ) = (R, A)
1
pR,A (r, a) = 1 0 (r, a) r drda
π θ(D )
Pour approfondir notre étude du point (R, A), il convient de noter que la fonction
pR,A peut s’exprimer sous la forme d’un produit de densités
 
1
pR,A (r, a) = 2r 1[0,1] (r) dr ×
 
1[0,2π] (a) da
2π
Ceci nous informe que les deux composantes de rayon R et d’angle A sont
indépendantes, de densités respectives
1
pR (r) = 2r 1[0,1] (r) dr et pA (a) = 1 (a) da
2π [0,2π]
L’angle A est donc uniforme sur [0, 2π], le rayon aléatoire R est distribuée selon une loi
rendant plus probables les rayons proches de 1.
Si U désigne une variable aléatoire uniforme sur [0, 1] alors la variable (2πU ) à la
même loi que A, en effet pour tout a ∈ [0, 2π]
Z a
 a a 1
P (2πU ≤ a) = P U ≤ = = P (A ≤ a) = da0 .
2π 2π 2π 0
√
D’autre part, si U 0 désigne une variable uniforme sur [0, 1] alors U 0 a la même loi que
R. En effet pour chaque r ∈ [0, 1[ on a bien
√  Z r
∂ 2
U 0 ≤ r = P U 0 ≤ r2 = r2 = P (R ≤ r) =


P (s ) ds.
0 ∂s
En choisissant U et U 0 indépendantes, on a donc montré les équivalences en loi
loi
√
θ(X1 , X2 ) = (2πU, U 0 )

Par conséquent, si U et U 0 sont deux variables uniformes et indépendantes sur [0, 1]

alors le point (X1 , X2 ) de coordonnées
√
X1 = √U 0 cos(2πU )


X2 = U 0 sin(2πU )

est uniformément réparti sur le disque unité D. Bien entendu, nous pouvons vérifier
directement ce résultat en utilisant la transformation entre les les variables uniformes
(u, u0 ) sur [0, 1]2 , avec u ∈
6 { 14 , 34 }, et les points (x1 , x2 ) sur le disque donnée par la
formule
1 x2
u= arctan et u0 = x21 + x22
2π x1
Un simple calcul de dérivées permet d’obtenir les équations suivantes
∂u x2 1 x2 1
= − 2 2
=− 2
∂x1 2πx1 1 + (x2 /x1 ) 2πx1 1 + (x2 /x1 )2
∂u 1 1 1 1
= 2
=
∂x2 2πx1 1 + (x2 /x1 ) 2πx1 1 + (x2 /x1 )2
et
∂u0 ∂u0
= 2x1 , = 2x2
∂x1 ∂x2
On obtient alors aisément les les associations d’éléments de surface infinitésimale
recherchés :
∂u ∂u0 ∂u ∂u0 1
dudu0 − × dx1 dx2 = dx1 dx2
∂x1 ∂x2 ∂x2 ∂x1 π

Notons pour conclure que l’on peut facilement décrire les loi des coordonnées
cartésiennes X1 et X2 . Par exemple la loi du demi-cercle est la loi de X1
Z 1
1
P (X1 ∈ dx1 ) = P (X1 ∈ dx1 , X2 ∈] − 1, 1[) = ( 1D (x1 , x2 ) dx2 ) 1]−1,1[ (x1 ) dx1
π −1
Z √1−x2
1 2
1
q
= ( dx ) 1 (x ) dx = 1 − x21 1]−1,1[ (x1 ) dx1
π −√1−x21
2 ]−1,1[ 1 1
π

D’après ce qui précède on en conclut que la variable aléatoire

√
X1 = U 0 cos(2πU )

est distribuée selon la loi du demi-cercle.

2.2.5 L’algorithme Box-Muller

Les fluctuations aléatoires de nature gaussienne font vraisemblablement partie des
phénomènes aléatoires les plus fréquement observés dans la nature. Cette fréquence
s’explique en grande partie par le théorème central de la limite. Ce fameux résultat
probabiliste affirme que tout phénomène aléatoire résultant d’une accumulation de
petites fluctuations indépendantes de même nature est nécessairement de nature
gaussienne ! On utilise donc de tels aléas gaussiens en engéniérie pour décrire les
fluctuations d’erreurs de capteurs électroniques, ou tout autre type d’erreurs de
modèles. En physique, les répartitions de chaleur, ou de fluides peuvent s’interpréter
comme des excitations et des collisions de particules. Les modèles cinétiques utilisés
dans ce contexte sont souvent formés de particules browniennes évoluant dans l’espace
selon des processus de diffusion de type gaussien.
Il existe une bonne demi-douzaine d’algorithmes permettant de simuler des
phénomènes gaussiens. Le plus connu est sans nul doute l’algorithme de Box-Muller. La
célébrité de cette technique provient de sa rapidité d’exécution. En effet, l’algorithme
de Box-Muller permet de construire deux variables indépendantes et gaussiennes
simplement à l’aide de deux variables uniformes et indépendantes sur [0, 1].
Nous avons choisi de présenter cet algorithme en utilisant l’étude de la uniforme sur
le disque unité D décrite dans la section 2.2.4. On note (R, A) les coordonnées polaires
d’un point choisi au hasard dans D. Nous avons montré précédemment que ces variables
aléatoires sont indépendantes, et distribuées selon la densité
 
R,A
  1
p (r, a) = 2r 1]0,1[ (r) 1 (a)
2π ]0,2π[
On opère maintenant une transformation non linéaire sur le rayon
p
S = −4 log R

On rappelle que l’on peut simuler le rayon R à l’aide d’une variable uniforme V sur
]0, 1[ en posant √
R= V
Par conséquent, le point S peut être simulé en posant
p p
S = −2 log R2 = −2 log V

La loi de S est donnée par un calcul élémentaire sur la fonction de répartition

s2
e− 4
∂r2
Z
−s2 /4 2
P(S ≥ s) = P(R ≤ e )= dr = e−s /2
0 ∂r
Par conséquent S est distribuée sur la demi-droite réelle R+ = [0, ∞[ avec la densité
∂ s2
pS (s) = − P(S ≥ s) = s e− 2 1R+ (s)
∂s
On en déduit que les variables (S, A) sont indépendantes distribuées sur (R+ ×]0, 2π[)
avec la densité    
S,A
2
− s2 1
p (s, a) = s e 1R+ (s) 1 (a)
2π ]0,2π[
Examinons la distribution dans l’espace du point

X1 = S cos A
X2 = S sin A
Pour déterminer cette loi de probabilité, il convient tout d’abord de noter que le point
(X1 , X2 ) est le transformé du point (S, A)
(X1 , X2 ) = θ(S, A)
avec la fonction
   
π 3π
θ(s, a) = (s cos a, s sin a) ∈ R2 − { {0} × R}


θ : (s, a) ∈ R+ × ]0, 2π[− ,
2 2
On obtient aisément les formules d’inversion
 q 
θ(s, a) = (x1 , x2 ) ⇐⇒ a = arctan (x2 /x1 ) 2 2
et s = x1 + x2

Les variations des coordonnées (a, s) en fonction du couple (x1 , x2 ) sont données ci-après
 
∂a 1 x2 ∂a 1 1
= 2
× − 2 , = 2
×
∂x1 1 + (x2 /x1 ) x1 ∂x2 1 + (x2 /x1 ) x1
et
∂s x1 ∂s x2
=p 2 , =p 2
∂x1 x1 + x22 ∂x2 x1 + x22
On en déduit les associations de densité de probabilité
1 1 ∂s ∂a ∂s ∂a
q
1
−s2 /2 2 2
se dsda x21 + x22 e− 2 (x1 +x2 ) − dx1 dx2
2π 2π ∂x1 ∂x2 ∂x2 ∂x1
1 1 2 2
= e− 2 (x1 +x2 )
2π
Autrement dit, X1 et X2 sont deux variables gaussiennes centrées et normées
indépendantes. L’algorithme de Box-Muller correspond à l’expression de ce résultat
en terme de variables uniformes sur [0, 1]. Si (U, V ) sont deux variables indépendantes
et uniformes sur [0, 1] alors
p
X1 = −2 log U cos (2πV )
et p
Y2 = −2 log U sin (2πV )
sont deux variables gaussiennes centrées et normées indépendantes. Pour générer sur
ordinateur un couple de variables indépendantes gaussiennes (X1 , X2 ) de moyennes
(m1 , m2 ) et d’écarts type (σ2 , σ2 ) on posera tout simplement
p
X1 = m1 + σ1 −2 log U cos (2πV )
et p
X2 = m2 + σ2 −2 log U sin (2πV ).
Exercice 2.2.1 Nous proposons de montrer que la projection de la loi gaussienne dans
R2 sur le cercle unité devient uniforme.
Soit (X1 , X2 ) un couple de v.a. p gaussiennes, indépendantes, centrées et normées
sur R. Montrer que la v.a. R = X12 + X22 est indépendante du point aléatoire P =
( √ X2 1 2 , √ X2 2 2 ), qui est uniformément distribué sur le cercle C de rayon 1 (centré
X1 +X2 X1 +X2
en l’origine) de R2 . Calculer la loi exacte de R.

(X1,X2)
x1

!1 0 1 x2

Fig. 2.6 – Point uniforme sur le cercle

Exercice 2.2.2 Soit (X1 , X2 , X3 ) un triplé de p v.a. gaussiennes, indépendantes,

centrées et normées sur R. Montrer que la v.a. R = X12 + X22 + X32 est indépendante
du point aléatoire P = ( √ 2 X1 2 2 , √ 2 X2 2 2 , √ 2 X3 2 2 ), qui est uniformément
X1 +X2 +X3 X1 +X2 +X3 X1 +X2 +X3
distribué sur la sphère S2 de rayon 1, et centré en l’origine, de R3 . Calculer la loi exacte
de R.
 x3
(X1,X2,X3)

0
!1 1 x1

x2

Fig. 2.7 – Point uniforme sur la sphère

Exercice 2.2.3 Soit (X1 , . . . , Xn ) une suite de n v.a. q

gaussiennes, indépendantes,
Pn 2
centrées et normées sur R. Montrer que la v.a. Rn = i=1 Xi est indépendante
du point aléatoire Pn = ( qPX
n
i
)1≤i≤n , qui est uniformément distribué sur la sphère
j=1 Xj2
S n−1 de rayon 1, et centré en l’origine, de Rn . Calculer la loi exacte de Rn , ainsi que
le volume de la sphère S n−1 .
Pn
Exercice 2.2.4 1. Vérifier que la somme χ2n = 2
i=1 Xi , de n carrés de v.a.
gaussiennes centrées et normées (Xp )1≤p≤n , suit une loi gamma de paramètres
(1/2, n/2)
 n/2
χ2n 1 1 n
P (dx) = 1R+ (x) e−x/2 x 2 −1 dx
2 Γ(n/2)
Montrer que E(χ2n ) = n, et Var(χ2n ) = 2n. En statistique mathématique, les lois
gamma de paramètres (1/2, n/2), sont appelés des Chi-deux à n degrés de liberté.
2. Montrer que la v.a. χ2n est indépendante du vecteur normalisé
p
(Xi / χ2n )1≤i≤n
qP
n 2
et vérifier que χn = i=1 Xi est distribué selon la loi de probabilité

1 x2
Pχn (dx) = 1R+ (x) n
−1
xn−1 e− 2 dx
2 2 Γ(n/2)
2.3 Méthode d’acceptation-rejet
Nous allons décrire dans cette section est une méthode de simulation universelle.
Cette technique est fondée sur un mécanisme d’apprentissage naturel et très simple.
On propose dans un premier temps des variables aléatoires distribuées selon une loi de
référence dominante et facile à simuler. Dans un second temps on accepte ou on refuse
ces variables. Le critère d’acceptation est simplement basé sur une comparaison de cette
mesure de référence avec la distribution de la variable que l’on souhaite simuler. Cette
technique universelle permet notamment de simuler tout phénomène aléatoire à valeurs
dans une régions compacte d’un espace topologique, et distribué sur cette région selon
une densité de probabilité bornée. On illustrera ce procédé pour la simulation de lois
telles que les lois conditionnelles ou la loi du demi-cercle et on présentera l’algorithme
polaire de simulation d’une loi gaussienne.

2.3.1 Les lois de référence

Commençons par nous donner un couple de densités de probabilités (p, q) sur un
espace Rd , avec d ≥ 1. On dit que la densité p domine q lorsque l’on a
q(x) ≤ m × p(x)
pour une certaine constante positive m > 0. Notons que dans cette situation, le
support de q doit être contenu dans celui de p. Autrement dit, lorsque p domine p
on a nécessairement
p(x) = 0 =⇒ q(x) = 0
Pour fixer les idées, on notera que les couples de densités suivantes vérifient cette
propriété
√
1) p(x) = 21 1[−1,1] q(x) = π2 1 − x2 1[−1,1] (x)(x)

1 1
2) p(x1 , x2 ) = 4 1[−1,1]2 (x1 , x2 ) q(x1 , x2 ) = π 1{(y1 ,y2 )∈R2 : y12 +y22 <1} (x1 , x2 )

x2
3) p(x) = 1
π (x2 +1)
q(x) = √1
2π
e− 2

Le premier exemple correspond à la loi du demi-cercle q, et la densité uniforme p sur

l’intervalle [−1, 1]. Dans cette situation, nous avons
q(x) 4
m= sup =
x∈[−1,1] p(x) π

Dans le second cas, p est le densité uniforme sur le carré [−1, 1]2 , et q la mesure uniforme
sur le disque unité. Dans cette situation, on peut noter que l’on a encore
q(x1 , x2 ) 4
m= sup =
(x1 ,x2 )∈[−1,1]2 p(x1 , x2 ) π
Le dernier exemple avec les densités gaussienne et de Cauchy, souligne le fait que
la propriété précédente est automatiquement satisfaite pour des densités strictement
positives. Dans cette situation, il convient de noter que
r
∂ π 2
(q/p) (x) = x(1 − x2 ) e−x /2
∂x 2
La fonction q/p est donc croissante sur ] − ∞, −1], décroissante sur [−1, 0], puis à
nouveau croissante sur [0, 1], et enfin décroissante sur [1, ∞[. On a de plus
r
2π
lim (q/p)(x) = 0 , (q/p)(−1) = = (q/p)(1) , et enfin (q/p)(0) = 1
x→−/+∞ e
Le maximum de (q/p) est donc atteint aux points x ∈ {−1, 1}. Autrement dit, nous
avons r
2π
m = sup (q/p)(x) =
x∈R e
On peut aussi noter que
R les restrictions q(x) d’une densité donnée p(x) sur des
d
régions A ⊂ R (telles que A q(y)dy > 0) vérifient clairement cette propriété
1
p(x) = 0 =⇒ q(x) = R 1A (x) p(x) = 0
A p(y)dy

Un contre exemple est facile à construire, le couple de densités uniformes suivant

1 1
p(x) = 1 (x) et q(x) = 1 (x)
2 [−1,1] 4 [−2,2]
ne vérifie pas cette propriété simplement à cause du fait que le support [−1, 1] de la
densité p ne contient pas le support [−2, 2] de p.

2.3.2 Les taux d’acceptation

Simulons une variable aléatoire X1 distribué selon la densité dominate p. Par
construction, nous avons
1 q(X1 )
0≤ × ≤1
m p(X1 )
Autrement dit, le segment [0, 1] est partagé en deux intervalles
   
1 q(X1 ) 1 q(X1 )
[0, 1] = 0, × ∪ × ,1
m p(X1 ) m p(X1 )
On simuleh alors un nombre
i uniforme U1 sur [0, 1]. Si ce nombre U1 appartient au premier
1 q(X1 )
segment 0, m × p(X1 ) , on garde X1 , et l’on pose

Y = X1
Dans le cas contraire, on recommence toute l’opération. On simule une nouvelle variable
aléatoire X2 distribué selon la densité dominate p, et
h une variable i uniforme U2 sur [0, 1].
1 q(X2 )
Si ce nombre U2 appartient au premier segment 0, m × p(X2 ) , on garde X2 , et l’on
pose
Y = X2
Dans le cas contraire, on reproduit entièrement le procédé. On simule une nouvelle
variable aléatoire X3 distribué selon la densité dominate p,het une variable
i uniforme U3
1 q(X3 )
sur [0, 1]. Si ce nombre U3 appartient au premier segment 0, m × p(X3 ) , on garde X3 ,
et l’on pose
Y = X3
Dans le cas contraire, on réitère le procédé, jusqu’au premier instant Tm où l’on finit
par accepter la variable simulée selon la loi dominate

Y = XTm

En terme mathématiques, ce procédé de simulation revient à produire une suite de

variables Xn indépendantes et distribuées selon la densité dominante, ainsi qu’une suite
Un de variables indépendantes et uniformes sur [0, 1]. On note ensuite
 
1 q(Xn )
Tm = inf n ≥ 1 : Un ≤ ×
m p(Xn )

1 q(Xn )
le premier instant où l’on a Un ≤ m × p(X n)
. La magie de ce procédé de simulation est
résumée dans les deux lignes suivantes. Nous venons de construire un couple (XTm , Tm )
de variables indépendantes telles que :
– La variable aléatoire XTm est distribuée selon la densité q.
1
– Le temps aléatoire d’acceptation Tm suit une loi géométrique de paramètre m .
Autrement dit, nous avons pour tout n ≥ 1

1 n−1
 
1
P(Tm = n) = 1− .
m m

En particulier le nombre moyen de boucles nécessaires pour obtenir une réalisation

de la variable XTm est E(Tm ) = m.

2.3.3 Vérification mathématique

La vérification mathématique du tour de magie présenté dans la section précédente
ne nécessite aucun raisonnement analytique. La dynamique de la preuve est purement
ptobabiliste !
 On commence par remarquer que pour chaque variable proposée selon la loi
dominante, disons X
h n , la probabilité
i pour que la variable uniforme Un tombe dans
1 q(Xn )
le premier segment 0, m × p(Xn ) est donnée par la probabilité conditionnelle
 
1 q(Xn ) 1 q(Xn )
P Un ≤ × | Xn = ×
m p(Xn ) m p(Xn )
En intégrant sur toutes le valeurs que peut prendre Xn , on obtient
  Z Z
1 q(Xn ) 1 q(xn ) 1 1
P Un ≤ × = × p(xn ) dxn = q(xn ) dxn =
m p(Xn ) m p(xn ) m m
Autrement dit, lors du déroulement dynamique de l’algorithme stochastique, la
probabilité d’un succés est égale à 1/m. Inversement, la probabilité d’échec et de refus
d’une variable est donnée par la formule suivante
 
1 q(Xn ) 1
P Un ≥ × =1−
m p(Xn ) m
Dans la pratique, il est donc très important de bien choisir la loi dominante p telle que
le nombre m soit le plus proche possible de 1. L’idéal étant bien entendu de choisir
p = q, de sorte à avoir m = 1 ; mais la distribution cible q est souvent trop complexe à
simuler. D’où l’intérêt de ce tour de magie.

Pour chaque région régulière A ⊂ Rd , et pour chaque entier n ≥ 1, nous avons

aussi
1 q(Xn ) 1 q(Xn )
P(Xn ∈ A , Un ≤ × ) = = E( 1A (Xn ) P(Un ≤ × | Xn ) )
m p(Xn ) m p(Xn )
 
1 q(Xn )
= E 1A (Xn ) ×
m p(Xn )
Z Z
1 q(xn ) 1
= × p(xn ) dxn = q(xn ) dxn .
A m p(xn ) m A
On utilise maintenant les décompositions évènementielles suivantes
{XTm ∈ A} = {Il existe un instant n ≥ 1 tel que Tm = n et Xn ∈ A}
= ∪n≥1 {Xn ∈ A , Tm = n}
et
{Xn ∈ A , Tm = n}
n o
1 q(Xn ) 1 q(Xk )
= Xn ∈ A , Un ≤ m p(Xn ) et pour les instants précédents 1 ≤ k < n on a Uk > m p(Xk )
n o  
1 q(Xk )

1 q(Xn ) n−1
= Xn ∈ A , Un ≤ m p(Xn ) ∩ ∩k=1 Uk > m p(X )
k
(2.3)
D’après les deux premières remarques, et en utilisant l’indépendance des couples de
variables aléatoires (Uk , Xk ), on obtient

1 n−1
Z  X   Z
1
P(XTm ∈ A) = q(x)dx 1− = q(x)dx
m A m A
n≥1

En prenant A = Rd dans (2.3) on trouve finallement

 n−1 
1 n−1
   
1 q(Xn ) Y 1 q(Xk ) 1
P(Tm = n) = P Un ≤ P Uk > = 1− .
m p(Xn ) m p(Xk ) m m
k=1

L’indépendance entre les variables XTm et Tm provient simplement du fait que

P(XTm ∈ A , Tm = n) = P(Xn ∈ A , Tm = n)
1 n−1
Z   
1
= q(x)dx 1−
A m m
= P(XTm ∈ A) P(Tm = n)
1
En posant z = m, on a enfin
 
X ∂  X
E(Tm ) = n(1 − z)n−1 z = −z (1 − z)n 
∂z
n≥1 n≥1
 
∂ 1 1
= −z −1 = =m
∂z z z

2.4 Files d’attentes exponentielles

Cette section concerne la simulation d’arrivées aléatoires de clients virtuels dans une
file d’attente d’un magasin tout aussi virtuel. Ce type de phénomène aléatoire se prête à
diverses interprétations. En effet, les instants d’arrivées des clients correspondent tout
simplement à des dates aléatoires. Ces dernières peuvent alors s’interpréter comme
des durées de fonctionnement de machines, ou des instants de panne, des durées de
connexions, des temps de traitement d’information, des arrivées d’automobiles dans les
bouchons du périphérique parisien, etc.
La simulation sur ordinateur de tels processus est assez claire. Il suffit de simuler
les arrivées successives de chacun des clients. Les écarts temporels entre les arrivées
de deux clients successifs peuvent être de nature très variée. Ces instants peuvent être
distribués sur l’axe des temps selon des lois exponentielles, des lois de Rayleigh, de
Weibull, etc.
La section 2.4.1 concerne la simulation de formations aléatoires des files d’attente
exponentielles. Ces processus d’arrivées de clients à des dates aléatoires exponentielles
sont d’un usage très fréquent en pratique. Ces horloges aléatoires apparaissent
notamment de façon naturelle dans la construction de processus markoviens à temps
continu sautant de temps à autre, et à des dates exponentielles sur d’autres site de
l’espace d’état. Dans la section 2.4.2, nous aborderons la distribution des tailles de files
d’attente exponentielles.

2.4.1 La loi Gamma

La simulation d’une file d’attente où les clients arrivent aléatoirement à une
fréquence d’exponentielle est très simple à réaliser. En effet, il suffit de se donner une
suite de variables aléatoires indépendantes Sn de même loi exponentielle

P(Sn ∈ ds) = λ e−λ s

1[0,∞[ (s) ds avec λ>0

Dans la section 1.4.5 consacrée aux horloges exponentielles, nous avons vu que le
paramètre λ représente l’intensité d’arrivée des clients. Plus précisément, nous avons
E(Sn ) = 1/λ. La fréquence d’arrivée des clients est donc d’autant plus grande que le
paramètre λ est petit.
Chaque variable Sn représente le temps d’attente aléatoire entre l’arrivée du (n −
1)ième client, et le suivant. Autrement dit, si l’on note Tn la date d’arrivée du nième client,
nous avons
Tn = Tn−1 + Sn = S1 + . . . + Sn−1 + Sn avec T0 = 0
Deux questions naturelles viennent à l’esprit : Quelles est la loi des instants d’arrivée
du groupe de clients, et quelle est la loi d’arrivée du nième individu ?
Pour calculer la loi du vecteur aléatoire (T1 , . . . , Tn ), on se donne une fonction f
suffisament régulière de [0, ∞[n dans R, et l’on commence par remarquer que l’on a
Z
E(f (T1 , . . . , Tn )) = f (t1 , . . . , tn ) p(T1 ,...,Tn ) (t1 , . . . , tn ) dt1 . . . dtn
= E(f (S1 , S1 + S2 , . . . , S1 + . . . + Sn ))
Z ∞ Z ∞
= ... f (s1 , s1 + s2 , . . . , s1 + . . . + sn ) λn e−λ (s1 +...+sn )
ds1 . . . dsn
0 0

Pour identifier la densité de probabilité p(T1 ,...,Tn ) du vecteur aléatoire (T1 , . . . , Tn ), on

effectue le changement de variable


 t1 = s1
 t2 = s1 + s2



t3 = s1 + s2 + s3
 ..
.




tn = s1 + s2 + . . . + sn

Par construction le vecteur des temps (t1 , . . . , tn ) est formée de nombres croissants.
Plus précisément, nous avons

(t1 , . . . , tn ) ∈ Cn = {(t1 , . . . , tn ) ∈ Rn+ ; t1 < . . . < tn }.

Inversement, nous avons 


 s1 = t1
 s2 = t2 − t1



s3 = t3 − t2
 ..



 .
sn = tn − tn−1


En effectuant un simple changement de variable dans chaque intégrale, on obtient les

associations infinitésimales

ds1 = dt1 , ds2 = dt2 , ... dsn = dtn

et finallement
Z
E(f (T1 , . . . , Tn )) = f (t1 , . . . , tn ) λn e−λ tn
dt1 . . . dtn
Cn

En résumé, la loi du vecteur (T1 , . . . , Tn ) est donnée par

P ((T1 , . . . , Tn ) ∈ d(t1 , . . . , tn )) = λn e−λ tn

1Cn (t1 , . . . , tn ) dt1 . . . dtn

Pour calculer la loi des arrivées alátoires Tn du nième client, on commence par noter que

P(Tn ∈ dtn ) = P (Tn ∈ dtn , (T1 , . . . , Tn−1 ) ∈ Cn−1 (tn )) (2.4)

"Z #
= λn e−λ tn dtn × dt1 . . . dtn−1
Cn−1 (tn )

avec les ensembles

n−1

Cn−1 (tn ) = (t1 , . . . , tn−1 ) ∈ R+ ; 0 < t1 < . . . < tn−1 < tn
= {(t1 , . . . , tn−1 ) ∈ Cn−1 tels que tn−1 < tn }

Il nous reste donc à calculer les intégrales de la forme suivante :

Z
In (t) = 1Cn (t) (t1 , . . . , tn ) dt1 . . . dtn
Z
= dt1 . . . dtn
0<t1 <...<tn−1 <tn <t
Z "Z t
# Z t
= dt1 . . . dtn−1 dtn = In−1 (tn ) dtn
0 0<t1 <...<tn−1 <tn 0
Une simple récurrence nous permet alors de vérifier que l’on a
Z t Z t
t2
I1 (t) = ds = t ⇒ I2 (t) = I1 (s) ds =
0 0 2!
Z t
t3
⇒ I3 (t) = I2 (s) ds =
0 3!
..
.
tn
⇒ In (t) =
n!
En utilisant (2.4), on obtient finallement
P(Tn ∈ dt) = In−1 (t) × λn e−λ t 1[0,∞[ (t) dt
(t)n−1
= × λn e−λ t 1[0,∞[ (t) dt
(n − 1)!
(λt)n−1
= λ e−λ t × 1 (t) dt
(n − 1)! [0,∞[
Lorsque λ = 1, les arrivées aléatoires du (n + 1)ième client sont distribuées selon la loi
Gamma de paramètre (n + 1) définie par la formule suivante
tn
P(Tn+1 ∈ dt) = e−t 1 (t) dt (2.5)
n! [0,∞[

2.4.2 La loi de Poisson

Les distributions de Poisson sont des lois de comptage du nombre de fois où un
évènement aléatoire se réalise sur un intervalle de temps donné. Ces lois sont souvent
utilisées pour estimer le nombre de clients en attente dans un réseau de communication,
le nombre de pannes d’une machine dans une journée, le nombre de mutations
génétiques au cours de l’évolution d’une espèce, etc. En termes mathématiques, on
dit qu’une variable aléatoire N à valeurs entières suit une loi de Poisson de paramètre
γ lorsque l’on a
γn
P(N = n) = e−γ
n!
pour chaque entier n ∈ N. Le paramètre γ correspond à la valeur moyenne de la variable
de Poisson N
 
X γ n X γ n−1
E(N ) = e−γ n =γ  e−γ =γ
n! (n − 1)!
n≥1 n≥1

La forme de la loi Gamma introduite en (2.5) est très proche de la loi de Poisson.
L’inversion des rôles des paramètres t et n est basée sur les formules suivantes
Z +∞
γ n −γ n
 
∂ −t t
e = −e dt
n! γ ∂t n!
et
n tn tn−1
 
∂ −t t
−e = e−t − e−t .
∂t n! n! (n − 1)!
Lorsque les variables aléatoires Sn sont distribuées selon une loi exponentielle de
paramètre λ = 1 alors d’après (2.5) on en déduit que
Z +∞ Z +∞
γ n −γ −t t
n tn−1
e = e dt − e−t dt
n! γ n! γ (n − 1)!
= P(Tn+1 > γ) − P(Tn > γ)
= P(Tn+1 > γ et Tn ≤ γ).

La dernière égalité provient du fait que l’on a

A = {Tn > γ} ⊂ B = {Tn+1 > γ} et P(B) = P(B − A) + P(A).

Ainsi, en posant
N = inf{n ≥ 1 ; Tn+1 > γ}
on a clairement
{N = n} = {Tn ≤ γ et Tn+1 > γ}
Cette simple équivalence évènementielle nous conduit au procédé de simulation d’une
variable de Poisson :
Si Sn sont des variables aléatoires indépendantes de loi exponentielle de paramètre
λ = 1 alors pour tout γ > 0 la variable aléatoire entière suivante

N = inf{n ≥ 1 ; S1 + . . . + Sn+1 > γ}

est distribuée sur N selon une loi de Poisson de paramètre γ.

2.4.3 La statistique d’ordre uniforme

Examinons de plus près la répartition des dates d’arrivées relatives de (n + 1)
visiteurs Gamma décrites dans la section 2.4.1
 
T1 T2 Tn
(V1 , . . . , Vn ) =déf. , ,..., (2.6)
Tn+1 Tn+1 Tn+1

Ces dates sont rangées dans un ordre croissant dans l’intervalle [0, 1]. Le dernier visiteur
sert de référence temporelle. Comment se distribuent ces dates aléatoires relatives ?
Les premiers visiteurs arrivent-ils généralement plus rapidement que les derniers ?
Intuitivement non, toutes ces dates semblent uniformément réparties sur l’intervalle
[0, 1]. Si cela était le cas, il serait équivalent de simuler une suite indépendantes de
dates uniformément réparties sur [0, 1], puis de les ranger ! Cependant, les coûts de
calculs pour ordonner cette longue séquence d’instants, uniformes sur [0, 1], seraient
bien plus élevés que ceux nécessaires à la simulation des arrivées exponentielles relatives
décrites plus haut. De plus, si ces deux phénomènes aléatoires étaient équivalents,
cela impliquerait que les dates d’arrivées relatives (2.6) des visiteurs Gamma seraient
indépendantes de l’instant d’arrivée Tn+1 du dernier visiteur ! Bien que les conséquences
de notre hypothèse d’équivalence statistique semblent déroutantes, nous allons vérifier
par de simples calculs que notre intuition est bien cohérente.
Pour poursuivre notre étude, il convient de se rappeller que la distribution (2.4) du
vecteur (T1 , . . . , Tn+1 ) peut s’écrire sous la forme suivante :
P ((T1 , . . . , Tn+1 ) ∈ d(t1 , . . . , tn+1 ))

= 1Cn+1 (t1 , . . . , tn+1 ) λn+1 e−λ tn+1 dt1 . . . dtn+1

 tn+1 (λ tn+1 )
n

= n! 1Cn (tn+1 ) (t1 , . . . , tn ) t−n
n+1 dt1 . . . dtn λ e−λ n! 1[0,∞[ (tn+1 ) dtn+1

= P ((T1 , . . . , Tn ) ∈ d(t1 , . . . , tn ) | Tn+1 = tn+1 ) × P(Tn+1 ∈ dtn+1 )

avec les ensembles désormais classiques
Cn (t) = {(t1 , . . . , tn ) ∈ Rn+ ; 0 < t1 < . . . < tn < t}.
Afin de trouver la loi du vecteur (V1 , . . . , Vn ), on se donne une fonction suffisament
régulière et bornée sur [0, 1]n , et l’on note que
   
T1 T2
E f Tn+1 , Tn+1 , . . . , TTn+1
n
| Tn+1 = tn+1
Z  
t1 t2 tn
= f , ,..., n! t−n
n+1 dt1 . . . dtn
Cn (tn+1 ) tn+1 tn+1 tn+1
Un simple changement de variable
 
t1 t2 tn
=⇒ dv1 . . . dvn = t−n


v1 = , v2 = , . . . , vn = n+1 dt1 . . . dtn
tn+1 tn+1 tn+1
permet de vérifier que l’espérance conditionnelle précédente est indépendante de la date
d’arrivée du (n + 1)ième visiteur ! Plus précisement, on obtient
   
T1 T2
E f Tn+1 , Tn+1 , . . . , TTn+1
n
| Tn+1 = tn+1
Z
= f (v1 , . . . , vn ) n! dv1 . . . dvn
Cn (1)

D’aprés les calculs d’intégrales effectués à la page 56, on peut observer que l’on a
Z
1
dv1 . . . dvn = = P((U1 , . . . , Un ) ∈ Cn (1))
Cn (1) n!
où (U1 , . . . , Un ) désignent une suite de variables aléatoires indépendantes et
indentiquement distribuée selon une loi uniforme sur l’intevalle [0, 1]. Le vecteur des
dates relatives  
T1 T2 Tn
(V1 , . . . , Vn ) =déf. , ,...,
Tn+1 Tn+1 Tn+1
est donc indépendant de la date Tn+1 , et il se distribue sur l’intervalle [0, 1] comme
une suite ordonnée de nombres aléatoires choisis indépendamment et uniformément sur
[0, 1]. En résumé, nous avons montré qu’il est statistiquement équivalent de simuler et
ordonner une suite de variables uniformes et indépendantes sur [0, 1], ou de simuler
une suite d’arrivées exponentielles et relatives de visiteurs.
Chapitre 3

Chaı̂nes de Markov discrètes

3.1 Introduction
Un processus aléatoire est un phénomène dont une partie de l’évolution temporelle
est aléatoire. On rencontre ces processus dans divers domaines de la physique, ou
des sciences de l’ingénieur. Par exemple la répartition et l’évolution de la chaleur
dans un corps, la turbulence atmosphérique ; c’est aussi le cas des temps d’arrivée
d’appels téléphoniques, les erreurs de mesures dues aux perturbations thermiques dans
des capteurs électroniques de type radar ou sonar, ou encore l’évolution des cours de
marchés boursiers.
La théorie de processus aléatoires est une théorie mathématique très riche, offrant de
nombreuses interactions avec diverses branches des mathématiques, telles la théorie des
graphes, l’analyse fonctionnelle, la théorie des opérateurs, ainsi que la théorie ergodique
des systèmes dynamiques.

Formellement, un processus aléatoire est une succession de variables aléatoires

(Xn )n≥0
X0 −→ X1 −→ . . . −→ Xn −→ Xn+1 −→ . . .
L’indice n ∈ N représente la paramètre temporel. Les variables Xn peuvent
être à valeurs dans un espace discret E, dans des espaces euclidiens E = Rd , ou
dans des espaces plus complexes tels des espaces de chemins ou d’excursions.

Ces modèles trajectoriels sont utiles en biologie dans l’analyse d’évolutions ancestrales
de population. On rencontre aussi des processus à valeurs dans des espaces de matrices
en physique statistique, ou en analyse d’images. L’étude de tels modèles sort bien
entendu du cadre de cet ouvrage.
Les états Xn peuvent évoluer aléatoirement au cours du temps selon des mécanismes
plus ou moins complexes. Dans certains cas, les épreuves Xn sont indépendantes ; c’est
le cas des séquences de lancers de dés, ou les successions de jets de pièces de monnaies.

61
 Dans d’autres situations, le processus aléatoire est donné par une équation
physique récursive de la forme suivante

Xn = Fn (Xn−1 , Un )

La v.a. X0 désigne la condition initiale du système, Fn des fonctions de

dérive déterministes, et enfin Un des variables aléatoires “indépendantes” des
précédents états X0 , . . . , Xn−1 . De telles séquences sont appelées des chaı̂nes
de Markov.

En traitement du signal, de tels systèmes peuvent représenter l’évolution temporelle

d’un cible dans l’espace. Dans ce contexte, les variables aléatoires Un ont une double
dimension. Elles correspondent à la fois, aux erreurs de modélisation, ainsi qu’aux
stratégies inconnues de guidage.
En mathématiques financières, l’évolution des prix d’actifs dans des marchés
boursiers sont aussi modélisés par des chaı̂nes de Markov. Dans ce contexte, les variables
Un représentent les fluctuations, et la volatilité stochastique du marché financier.
Nous reviendrons sur des exemples plus précis dans la suite.

3.2 Chaı̂nes de Markov discrètes

Les chaı̂nes de Markov les plus élémentaires sont bien entendu celles dont les états
aléatoires ne prennent qu’un nombre fini, ou au plus dénombrable de valeurs. Cette
section est consacrée à l’étude de ces chaı̂nes élémentaires. Nous insisterons sur les
réalisations “canoniques”, et dynamiques, de ces processus aléatoires. La seconde partie
de cette section concerne l’étude des semigroupes d’évolution des ces chaı̂nes de Markov.
Afin de satisfaire la curiosité du lecteur, et souligner les points d’interaction avec
d’autres domaines scientifiques, sous soulignerons les interactions entre la théorie
des chaı̂nes de Markov, l’analyse fonctionnelle, l’algèbre matricielle, et la théorie des
graphes. Enfin, nous illustrerons ces nouvelles notions sur une variété d’exemples précis
liés à des phénomènes aléatoires issus de la physique, de la biologie, et des sciences de
l’ingénieur.
Dans la suite de cette section, Xn désigne une chaı̂ne de Markov discrète à valeurs
dans un espace E, au plus dénombrable. Ses transitions de probabilités seront données
par une collection (Mn (x, .))x∈E de mesures de probabilités sur E.
 Autrement dit, les probabilités de passage d’un état Xn−1 = x à un nouvel
état aléatoire Xn sont données par l’application suivante

y ∈ E 7→ Mn (x, y) = P(Xn = y | Xn−1 = x) ∈ [0, 1]

On utilise parfois la notation

PXn |Xn−1 (y|x) = P(Xn = y | Xn−1 = x)

On désigne par la suite ηn = PXn , la loi de l’état Xn de la chaı̂ne, à chacun des

instants n ∈ N.

On remarquera que la loi de la trajectoire (X0 , . . . , Xn ) de l’origine jusqu’à

l’instant n est alors décrite par la formule multiplicative

P(X0 ,X1 ,...,Xn ) (x0 , x1 , . . . , xn ) = P(X0 = x0 , X1 = x1 , . . . , Xn = xn )

= η0 (x0 )M1 (x0 , x1 ) . . . Mn (xn−1 , xn )

pour toute trajectoire (xp )0≤p≤n ∈ E n+1 .

On peut clairement étendre la notion de chaı̂ne de Markov précédente à des

modèles markoviens Xn prenant leurs valeurs dans des espaces En liés au paramètre
temporel ! Dans ce contexte, Mn+1 (xn , xn+1 ) désigne la probabilité de passer d’un
état xn ∈ En vers un état xn+1 ∈ En+1 . Plus formellement, nous avons à nouveau

Mn+1 (xn , xn+1 ) = P(Xn+1 = xn+1 | Xn = xn )

et

P(X0 ,X1 ,...,Xn ) (x0 , x1 , . . . , xn ) = P(X0 = x0 , X1 = x1 , . . . , Xn = xn )

= η0 (x0 )M1 (x0 , x1 ) . . . Mn (xn−1 , xn )

Bien que cette extension ne puisse paraı̂tre purement formelle, sans d’autres intérêts que
mathématiques, cette notion apparaı̂t naturellement dans la représentation de chaı̂nes
trajectorielles.
3.2.1 Semigroupes de transitions

Les probabilités de transitions Mn d’une chaı̂ne de Markov sur un espace E

sont associées à deux opérateurs naturels :
1. Le premier agit à droite sur les fonctions bornées sur E. A chacune de ces
fonctions f , on associe la fonction bornée Mn [f ] définie par la formule
suivante
X
Mn [f ](x) =déf. E(f (Xn )|Xn−1 = x) = Mn (x, y) f (y)
y∈E

2. Le second agit à gauche sur les mesures de probabilités sur E. A chacune

de ces mesures η, on associe la mesure de probabilité (ηMn ) définie par
X
(ηMn )(y) = η(x) Mn (x, y) ∈ [0, 1]
x∈E

Dans ce contexte, il est aussi très utile de voir une mesure de probabilité η comme un
opérateur sur l’ensemble des fonctions f bornées sur l’espace des états E
X
η[f ] = η(x) f (x) ∈ R
x∈E

Avec ce système de notations, pour les fonctions indicatrices f = 1B d’ensembles B ⊂ R,

on retrouve la mesure η(B) des ensembles B
X
η[1B ] = η(x) = η(B)
x∈B

De même, on observe que

Mn [1B ](x) = E(1B (Xn )|Xn−1 = x) = Mn (x, B) = P(Xn ∈ B | Xn−1 = x)

Lorsqu’il n’y a pas de confusions, il est coutume de noter η(f ), et Mn (f ) plutôt que
η[f ], et Mn [f ]. Dans ce système de notations, nous avons
X
E(f (Xn )) = f (x) ηn (x) = ηn (f )
x∈E

et X
E(f (Xn ) | Xn−1 = x) = Mn (x, y) f (y) = Mn (f )(x)
x∈E
En utilisant la formule des conditionnements emboı̂tes on montre que

ηn (f ) = E(E(f (Xn )|Xn−1 ))

= E( Mn (f )(Xn−1 ) )
X
= ηn−1 (x) Mn (f )(x) = ηn−1 (Mn (f ))
x∈E

Autrement dit, les lois ηn des différents états de la chaı̂ne de Markov Xn

peuvent “se calculer” récursivement. Ces dernières sont solution d’un système
dynamique discret à valeurs dans l’espace des probabilités :

ηn = ηn−1 Mn

Notre prochain objectif est de décrire plus précisement le semigroupe d’évolution

du flot de mesures (ηn )n . On introduit pour cela la probabilité conditionnelle de Xn en
Xp
Mp,n (x, y) = PXn |Xp (y|x) = P(Xn = y | Xp = x)
avec 0 ≤ p ≤ n, Comme précédemment, on associe à toute fonction f sur E, la fonction
Mp,n (f ) sur E donnée par
X
Mp,n (f )(x) = E(f (Xn )|Xp = x) = Mp,n (x, y) f (y) (3.1)
y∈E

En particulier pour p = n, on a Mn,n (f ) = f ; et pour p = (n − 1), nous avons

Mn−1,n (f )(x) = Mn (f )(x) = E(f (Xn )|Xn−1 = x)

D’après la formule de conditionnements emboı̂tés, nous avons la formule de

récurrence

Mp,n (f )(x) = E(E(f (Xn )|Xn−1 )|Xp = x)

= E(Mn (f )(Xn−1 )|Xp = x) = Mp,n−1 (Mn (f ))(x)

Par conséquent, l’opérateur de transition Mp,n est donné par la formule (3.1) avec

Mp,n (x, y)
P
= xp+1 ∈E,...,xn−1 ∈E Mp+1 (x, xp+1 )Mp+1 (xp+1 , xp+2 ) . . . Mn (xn−1 , y)

On en déduit une formulation des transitions conditionnelles Xp Xn , en terme

de compositions d’opérateurs

Mp,n (f ) = Mp+1 Mp+2 . . . Mn (f )

En utilisant la formule

∀p ≤ n E(f (Xn )) = E(E(f (Xn )|Xp )) = E(Mp,n (f )(Xp ))

on en conclut que
X X
ηn (f ) = ηp (Mp,n f ) = ηp (x)[ Mp,n (x, y)f (y)]
x∈E y∈E
X X
= (ηp Mp,n )(f ) = ( ηp (x)Mp,n (x, y)) f (y)
y∈E x∈E

L’ordre des sommes étant sans importance, on simplifie les notations, et on

écrit tout simplement

ηn = ηp Mp,n avec Mp,n = Mp+1 Mp+2 . . . Mn (3.2)

Nous avons donc montré que les opérateurs Mp,n , ≤ p ≤ n, correspondent au

semigroupe d’évolution du flot de mesures (ηn )n≥0 .
Lorsque la chaı̂ne est homogène, ces opérateurs correspondent à la composition du
même opérateur de transition M . Dans cette situation, on utilise souvent les notations
synthétiques suivantes

Mp,p+n = M n et ηn = η0 M n avec M n = M n−1 M = M M n−1

3.2.2 Processus historique

On considère une chaı̂ne de Markov discrète Xn0 à valeurs dans un espace au plus
dénombrable E 0 . On note Mn0 (x0 , y 0 ), la probabilité de passage de Xn−1 = x0 vers
Xn = y 0 . La séquence de trajectoires aléatoires

Xn = (X00 , . . . , Xn0 ) ∈ En =déf. (E 0 )n+1

forme à nouveau une chaı̂ne de Markov à valeurs dans les espaces trajectoriels En . En
effet, le passage de Xn à Xn+1 s’effectue en deux temps. On conserve tout d’abord le
segment de trajectoire Xn = (X00 , . . . , Xn0 ), puis on lui adjoint une extension élémentaire
0
Xn+1 = x0 de loi Mn+1
0 (Xn0 , x0 ). Autrement dit, nous avons

Xn = (X00 , . . . , Xn0 ) Xn+1 = ((X00 , . . . , Xn0 ), Xn+1

0 )
| {z }
=( Xn 0
, Xn+1 ) ∈ En+1 = (En × E 0 )
 On notera que pour toute fonction fn+1 bornée sur En+1 , et pour tout segment de
trajectoire xn = (x00 , . . . , x0n ) ∈ En , nous avons

E(fn+1 (Xn+1 ) | Xn = xn )

= E(fn+1 ([X00 , . . . , Xn ], Xn+1

0 ) | (X00 , . . . , Xn0 ) = (x00 , . . . , x0n ))

0
fn+1 ([x00 , . . . , x0n ], x0n+1 ) Mn+1 (x0n , x0n+1 )
P
= x0n+1 ∈E 0

Cette équation s’exprime aussi sous la forme trajectorielle suivante

E(fn+1 (Xn+1 ) | Xn = xn )
X
= fn+1 (y00 , . . . , yn+1
0
) 1(x00 ,...,x0n ) (y00 , . . . , yn0 )Mn+1
0
(yn0 , yn+1
0
)
(y00 ,...,yn+1
0 )∈En+1

La dernière formule montre que Xn est une chaı̂ne de Markov, de probabilités de

transitions Mn+1 de En vers En+1 , données par

Mn+1 ((x00 , . . . , x0n ), (y00 , . . . , yn+1

0
)) = 1(x00 ,...,x0n ) (y00 , . . . , yn0 )Mn+1
0
(yn0 , yn+1
0
)

3.2.3 Interprétation matricielle

Lorsque l’espace d’état E est, soit fini, soit indexé de façon naturelle par N, ou encore
par Z, les semigroupes de transitions Mp,n définis dans la section 3.2.1 correspondent à
des compositions “élémentaires” de matrices. Pour préciser cette assertion, commençons
par l’exemple le plus simple, où l’espace d’état est donné par un ensemble à deux points
E = {1, 2}. Dans ce cas, la donnée des transitions de la chaı̂ne

Mn (x, y) = P(Xn = y|Xn−1 = x)

est équivalente à la donnée des matrices (2 × 2) suivantes

 
Mn (1, 1) Mn (1, 2)
Mn = (Mn (x, y))x,y∈E =
Mn (2, 1) Mn (2, 2)

Par exemple, la matrice  

6/7 1/7
M=
1/2 1/2
correspond au mouvement aléatoire entre deux états E = {1, 2}, synthétisé par le
schéma ci-dessous
 6/7 1/7

1 2 1/2

1/2

Fig. 3.1 – Chaı̂ne à deux états

Par définition de l’opérateur de transition Mn , nous avons, pour toute fonction f sur
E = {1, 2}, la formule matricielle
   
Mn (f )(1) Mn (1, 1) f (1) + Mn (1, 2) f (2)
=
Mn (f )(2) Mn (2, 1) f (1) + Mn (2, 2) f (2)
   
Mn (1, 1) Mn (1, 2) f (1)
= (3.3)
Mn (2, 1) Mn (2, 2) f (2)
D’autre part, en utilisant la formule des conditionnements emboı̂tés
E(f (Xn )) = E(E(f (Xn )|Xn−1 )) = E(Mn (f )(Xn−1 ))
on obtient une nouvelle formule matricielle
 
ηn−1 (1) Mn (1, 1) ηn−1 (2) Mn (2, 1)
[ηn (1), ηn (2)] =
ηn−1 (1) Mn (1, 2) ηn−1 (2) Mn (2, 2)
 
Mn (1, 1) Mn (1, 2)
= [ηn−1 (1), ηn−1 (2)] (3.4)
Mn (2, 1) Mn (2, 2)
Par conséquent, si l’on représente une fonction numérique f , et une mesure de
probabilité η sur E par les vecteurs colonnes et lignes
 
f (1)
f= et η = [η(1), η(2)]
f (2)
alors, les équations (3.3) et (3.4) s’expriment sous la forme d’un semigroupe matriciel
M (f ) = M f et ηn = ηn−1 Mn = η0 M1 . . . Mn
Ces interprétations matricielles restent valables pour des chaı̂nes de Markov à valeurs
dans des espaces d’états finis abstraits E = {x1 , . . . , xd }. Dans ce contexte, les matrices
de transitions sont données par
 
Mn (x1 , x1 ) . . . Mn (x1 , xd )
Mn = 
 .. .. .. 
. . . 
Mn (xd , x1 ) . . . Mn (xd , xd )
Les mesures de probabilités η, et les fonctions f sur E sont associées au vecteurs lignes
et colonnes suivants
 
f (x1 )
η = [η(x1 ), . . . , η(xd )] et f = 
 .. 
. 
f (xd )

Exemple 3.2.1 La figure suivante

1
1/2 1/2

1/4 1/3
1/3 3
2
1/3
3/4

Fig. 3.2 – Chaı̂ne à 3 états

présente un schéma d’évolution de chaı̂ne de Markov sur un espace à trois points E =

{1, 2, 3}, et de matrice de transition homogène donnée par
 
0 1/2 1/2
M =  1/4 0 3/4 
1/3 1/3 1/3

Lorsque l’espace d’états est dénombrable, et indexé par les entiers positifs E = {xi , i ∈
N}, les matrices de transition sont infini-dimensionnelles
 
Mn (x0 , x0 ) M (x0 , x1 ) M (x0 , x2 ) . . .
 Mn (x1 , x0 ) M (x1 , x1 ) M (x1 , x2 ) . . . 
Mn = (Mn (x, y))x,y∈E =  M (x , x ) M (x , x ) M (x , x ) . . . 
 
 n 2 0 2 1 2 2 
.. .. .. ..
. . . .
Dans ce cas, les probabilités η, et les fonctions f sur E sont associés au vecteurs
 
f (x0 )
 f (x1 ) 
η = [η(x0 ), η(x1 ), η(x2 ), . . .] et f =  f (x ) 
 
 2 
..
.
Exemple 3.2.2 Une marche aléatoire homogène sur N, absorbée en 0, peut être
représentée par le schéma suivant

0 1 p 2 p 3 p 4
1

1−p 1−p 1−p 1−p

Fig. 3.3 – Chaı̂ne absorbée en 0

La matrice de transition associée à cette chaı̂ne est donnée par

 
1 0 0 0 ...
 (1 − p) 0 p 0 ... 
M =
 
 0 (1 − p) 0 p ... 

.. .. .. .. ..
. . . . .
Lorsque le point 0 est réfléchissant, la matrice de transition s’écrit sous la forme
 
0 1 0 0 ...
 (1 − p) 0 p 0 ... 
M =
 
 0 (1 − p) 0 p ...  
.. .. .. .. ..
. . . . .
Cette situation correspond au schéma suivant

0 1 1 p 2 p 3 p 4

1−p 1−p 1−p 1−p

Fig. 3.4 – Chaı̂ne réfléchie en 0

Dans le dernier cas, où l’espace d’états est dénombrable, et indexé par les entiers
relatifs E = {xi , i ∈ Z}, les matrices de transition sont données par
.. .. .. .. ..
 
 . . . . . 
 . . . M (x−1 , x−1 ) Mn (x−1 , x0 ) M (x−1 , x1 ) . . . 
 
Mn =   . . . M (x0 , x−1 ) Mn (x0 , x0 ) M (x0 , x1 ) . . . 

 . . . M (x1 , x−1 ) M n (x 1 , x0 ) M (x1 , x1 ) . . . 
 
.. .. .. .. ..
. . . . .
et les probabilités η, et les fonctions f sur E sont associés au vecteurs
..
 
 . 
 f (x−1 ) 
 
 f (x0 ) 
η = [. . . , η(x−1 ), η(x0 ), η(x1 ), . . .] et f =  
 f (x1 ) 
 
..
.

Exemple 3.2.3 Une marche aléatoire simple sur Z peut être représenté par le schéma
suivant ette situation correspond au schéma suivant

p p p p

−2 −1 0 1 2
1−p 1−p 1−p 1−p

Fig. 3.5 – Marche simple sur Z

La matrice de transition de cette évolution est donnée par

.. .. .. .. ..
 
 . . . . . 
 . . . 0 (1 − p) 0 . .. 
 
M =  ... p
 0 (1 − p) . . . 

 ... p 0 ... 
 
.. .. .. .. ..
. . . . .

3.3 Quelques Exemples

3.3.1 Files d’attentes
Soit (Un )n≥1 une suite de v.a. positives et indépendantes de lois respectives (µn )n≥1 .
On considère la chaı̂ne de Markov définie de façon récursive par l’équation suivante

Xn+1 = (Xn − 1)+ + Un+1



X0 = 0

Dans la formule précédente a+ = max (a, 0) désigne le maximum entre un nombre réel
a ∈ R et 0.
 Ce processus aléatoire peut s’interpréter comme la longueur d’une file
d’attente, ou encore le temps d’attente d’un client arrivant à un guichet,
servant une personne par unité de temps. Dans ce contexte, la v.a. Un+1
représente le nombre de clients arrivant dans la file d’attente au temps (n + 1).

On peut aussi interpréter Xn comme le nombre de paquets (symboles binaires

représentant de l’information : voix, vidéo, données,...) en attente dans la
mémoire d’un canal de communication, transmettant un paquet par
unité de temps. Dans cette situation, la v.a. Un+1 représente le nombre de
paquets arrivant dans le canal à l’instant (n + 1).

On notera que les transitions de cette chaı̂ne sont données pour tout i ≥ 1, et pour
tout j ≥ 0, par la formule suivante
P(Xn+1 = (i − 1) + j | Xn = i) = µn+1 (j) = P(Xn+1 = j | Xn = 0)

3.3.2 Modèle d’urne

On considère une urne contenant initialement B0 boules blanches, et N0 boules
noires. A chaque instant n, on choisit au hasard une boule, puis on remet cette boule
dans l’urne accompagnée d’une nouvelle boule de la même couleur. On note (Bn , Nn )
le nombre de boules blanches, et noires, dans l’urne au temps n. Par construction,
le couple Xn = (Bn , Nn ) est une chaı̂ne de Markov à valeurs dans E = N2 , et de
probabilités de transitions
b m
PXn |Xn−1 (.|(b, m)) = δ(b+1,m) + δ
m+b b + m (b,m+1)

3.3.3 Marche aléatoire sur Z

La marche aléatoire simple sur Z correspond à un mouvement aléatoire d’une
particule sur les entiers relatifs, se déplaçant soit d’un pas vers la droite, avec une
probabilité p, soit d’un pas vers la gauche, avec une probabilité (1 − p), avec p ∈ (0, 1).
Ce mouvement peut être représenté schématiquement par la figure suivante :

Plus formellement, ce mouvement aléatoire est défini par une chaı̂ne de Markov
X = (Xn )n≥0 définie sur un espace de probabilités (Ω, F, P), d’origine X0 = 0 et
de probabilités de transitions homogènes
M (x, y) =déf. P(Xn = y | Xn−1 = x) = PXn |Xn−1 (y|x)
= p 1x+1 (y) + (1 − p) 1x−1 (y)
 1−p 1−p 1−p 1−p

−2 −1 0 1 2

p p p p

Fig. 3.6 – Marche aléatoire simple sur Z

Pour décrire une réalisation dynamique de cette chaı̂ne, on se donne une suite de
v.a. indépendantes U = (Un )n≥1 distribuées sur {−1, +1} selon la même loi de Bernoulli

P(Un = +1) = 1 − P(Un = −1) = p

On suppose, comme d’habitude, que cette suite est définie sur un espace de probabilités
(Ω, F, P). On associe à U , le système dynamique aléatoire donné par

Xn = Xn−1 + Un
(3.5)
X0 = 0

Cette interprétation dynamique offre un certain nombre d’avantages. Par exemple, elle
permet une représentation “explicite”de l’état Xn en terme de v.a. indépendantes, et
simplifie l’analyse des transitions de probabilités.

Exercice 3.3.1 Cet exercice a pour objectif d’analyser plus en profondeur la marche
aléatoire sur Z décrite dans l’exemple 3.3.3.
1. Vérifier que l’interprétation dynamique introduite en (3.5) correspond bien à la
donnée d’une marche aléatoire simple sur Z.
2. Montrer que la position moyenne de la particule au temps n est donnée par la
formule E(Xn ) = n × (2p − 1). En conclure que

 limn→∞ E(Xn ) = −∞ si p ∈ [0, 1/2)
E(Xn ) = 0 si p = 1/2
limn→∞ E(Xn ) = +∞ si p ∈ (1/2, 1]


3. Vérifier que les transitions de la chaı̂ne entre deux instants, m et (m + n), sont
données par la formule

P(Xm+n = x + [k − (n − k)] | Xm = x) = Cnk pk (1 − p)n−k

pour tous les k ∈ {0, . . . , n}, et

P(Xm+n 6∈ {2k − n : k = 0, . . . , n}|Xn = x) = 0

 4. En déduire que
(2k)!
P(Xm+2k = 0 | Xm = 0) = (p(1 − p))k
k!k!
√
En utilisant la formule de Stirling (k! ' 2πk k k e−k ), montrer que

(4p(1 − p))k √
P(Xm+2k = 0 | Xm = 0) ' √ (= 1/ πk si p = 1/2)
πk

3.3.4 Marche aléatoire sur Zd

On note |.| la distance l1 sur E = Zd définies par
d
X
|x| = |xi |
i=1

pour tout x = (xi )1≤i≤d ∈ Zd . On associe à une mesure de probabilité p sur l’ensemble
des 2d vecteurs unitaires directionnels

U = {u ∈ Zd : |u| = 1}

la transition homogène X
M (x, y) = p(u) 1x+u (y)
u∈U

L’évolution aléatoire de la chaı̂ne Xn associée à M est claire. A chaque étape n, la

particule choisit aléatoirement un vecteur u ∈ U avec la probabilité p(u), et se déplace
dans cette direction. Autrement dit, si (Un )n≥1 désigne une suite de v.a. indépendantes
de même loi p sur U, on a une représentation dynamique de l’évolution
n
X
Xn = Xn−1 + Un = X0 + Un
i=1

La figure suivante présente une réalisation d’une trajectoire aléatoire de la chaı̂ne Xn

sur Z2 , d’origine X0 = 0 ∈ Z2 .

3.3.5 Marche aléatoire arrétée

Supposons qu’une particule évolue sur Zd selon les principes de transitions
élémentaires décrits dans l’exemple 3.3.4, mais cette dernière ne peut se mouvoir que
s’il elle se trouve dans une région spécifique de l’espace B ⊂ Zd . Autrement dit, lorsque
la chaı̂ne Xn sort de B, elle s’immobilise. Ce modèle physique peu à nouveau être
 0

Fig. 3.7 – Marche aléatoire sur Z2

représenté par une chaı̂ne de Markov Yn à valeurs dans Zd . Dans ce contexte, les
transitions Yn−1 Yn sont données par la formule suivante

Yn−1 + Un si Yn−1 ∈ B
Yn =
Yn−1 si Yn−1 6∈ B

La figure suivante présente une réalisation d’une trajectoire s’immobilisant à la sortie

d’un segment B de Z.

3.3.6 Processus de branchements

On considère une population d’individus se développant de la façon suivante. Notons
Xn le nombre d’individus à l’instant n. A l’étape suivante, chaque individu de label
i ∈ {1, 2, . . . , Xn }, donne naissance à Nni individus. On convient que ces nombres de
branchements (Nni )i≥1 sont indépendants des configurations passés (X0 , . . . , Xn−1 ).
Dans ce cas, le processus aléatoire des tailles de population forme un processus de
Markov à valeurs entières. Son évolution peut être décrite dynamiquement par le
système
Xn
X
Xn+1 = Nn1 + Nn2 + . . . + NnXn = Nni
i=1
 Z−B Xn

Xp
B

p n temps

Fig. 3.8 – Chaı̂ne stoppée

Lorsque les v.a. de branchement (Nni )i≥1, n≥1 sont des copies indépendantes d’une même
v.a. entière N , on notera que
Xn
X
E(Xn+1 |Xn ) = E(Nni ) = Xn × E(N ) ⇒ E(Xn+1 ) = E(X0 ) E(N )n+1
i=1

Par conséquent, en supposant que E(X0 ) 6= 0, la population moyenne s’éteindra lorsque

E(N ) < 1, et elle explosera lorsque E(N ) > 1. La figure ci-dessous présente une
réalisation d’un arbre de descendances d’un individu (X0 = 1).

Examinons la situation où chaque individu se dédouble avec une probabilité p, ou

disparaı̂t avec la probabilité (1 − p). Ce modèle correspond au choix d’une v.a. N de
Bernoulli loi
PN = p δ2 + (1 − p) δ0
Dans ce contexte, on notera que E(N ) = 2p. De plus, on montre les transitions de cette
chaı̂ne sont données par la formule binomiale

P(Xn+1 = 2y | Xn = x) = Cxy py (1 − p)x−y

pour tout y ∈ {0, . . . , x}, et P(Xn+1 ∈ 2N + 1 | Xn = x) = 0, pour tout x ∈ N.

X0=1 X1=2 X2=4

Fig. 3.9 – Processus de branchement

Chapitre 4

Chaı̂nes de Markov abstraites

4.1 Description des modèles

Comme nous l’avons souligné dans l’introduction, la notion de chaı̂ne de Markov,
est loin d’être restreinte à des phénomènes aléatoires à valeurs dans des espaces discrets.
Les fluctuations de température d’un liquide par unité de temps correspondent à des
chaı̂nes de Markov à valeurs réelles. L’évolution d’une cible spatiale peut être modélisée
par une chaı̂ne de Markov à valeurs dans l’espace euclidien R3 , ou bien dans R9 si l’on
tient compte des coordonnées de position, vitesse et accélération. Enfin, si l’on tient
compte des évolutions de cette même cible entre certains paliers, on obtient une chaı̂ne
de Markov à valeurs dans des espaces de chemins ou d’excursions.
Tous ces modèles physiques, s’inscrivent dans la théorie abstraite des chaı̂nes de
Markov à valeurs dans des espaces mesurables. Leur analyse s’exprime de façon naturelle
dans la théorie de l’intégration de Lebesgue. Afin d’éviter un langage trop technique,
j’ai volontairement choisi de ne pas inscrire cette introduction à l’ingénierie stochastique
dans ce cadre trop mathématique. Nous ignorerons donc les notions d’ensembles et de
fonctions mesurables, les théor‘emes de Fubini, et d’autres propriétés fondamentales de
la théorie de l’intégration de Lebesgue. La nature est suffisament stable pour polir les
erreurs et les confusions dues à de tels écarts. J’ai plutôt essayé de coller au plus près
à la réalité scientifique et technique, tout en aiguisant la curiosité mathématique du
lecteur.
Dans la suite, on utilisera donc la terminologie “chaı̂ne de Markov abstraite” pour
désigner une chaı̂ne de Markov Xn à valeurs dans des ensembles suffisament réguliers
En , de transitions de probabilités

Mn (xn−1 , dxn ) =déf. PXn |Xn−1 (dxn |xn−1 ) = P(Xn ∈ dxn | Xn−1 = xn−1 )

et de loi initiale η0 (dx0 ) = PX0 (dx0 ) sur E0 . On note alternativement

ηn (dxn ) = PXn (dxn ) = P(Xn ∈ dxn )

79
la loi de la v.a. Xn sur En , donnée pour toute fonction fn bornée sur En par la formule
Z
E(fn (Xn )) = fn (xn ) ηn (dxn ) =déf. ηn (fn )

Dans ce système de notations abusives, une mesure donnée sur En s’interprète, soit
comme un opérateur intégral, soit comme une fonction ensembliste. Ainsi, pour tout
sous ensemble suffisament régulier Bn ⊂ En , nous avons les représentations équivalentes

ηn (Bn ) = PXn (Bn ) = ηn [1Bn ]

Mn (xn−1 , Bn ) = PXn |Xn−1 (Bn |xn−1 )
= P(Xn ∈ Bn | Xn−1 = xn−1 ) = Mn [1Bn ](xn−1 )

On dit qu’une chaı̂ne de Markov réelle X = (Xn )n≥0 est absolument continue,
lorsque la v.a. initiale X0 est absolument continue de densité p0 sur R, et lorsque
ses transitions de probabilités sont données en terme d’une famille de densités de
probabilités {pn (x, .), x ∈ R} sur R

Mn (xn−1 , dxn ) = PXn |Xn−1 (dxn |xn−1 ) = pn (xn−1 , xn ) dxn (4.1)

Dans la formule précédente dxn désigne la mesure de Lebesgue sur R. A titre illustratif,
on pourra considérer les familles de densités gaussiennes, ou exponentielles données par
la formule
1 2
e− 2 (y−x)
pn (x, y) = √ ou pn (x, y) = 1[0,∞) (y) (|x| + 1) e−(|x|+1) y
2π
Dans ce contexte, la trajectoire (X0 , . . . , Xn ) de l’origine jusqu’à l’instant n est une v.a.
absolument continue, de loi donnée par la formule

P(X0 ,...,Xn ) (d(x0 , . . . , xn )) = p0 (x0 )p1 (x0 , x1 ) . . . pn (xn−1 , xn ) dx0 dx1 . . . dxn

4.1.1 Semigroupe des transitions

Comme dans le cas des chaı̂nes discrètes, il existe un système de notations naturel
permettant de décrire les semigroupes d’évolution des lois des états Xn .
 Les probabilités de transitions Mn (xn−1 , dxn ) permettent de définir deux
opérateurs intégraux naturels :
1. Le premier agit à droite sur les fonctions bornées. A chacune de ces
fonctions f , on associe la fonction bornée Mn [f ] définie par la formule
suivante

Mn [f ] : x ∈ R 7→ Mn [f ](x) =déf. E(f (Xn )|Xn−1 = x) ∈ R

2. Le second agit à gauche sur les mesures de probabilités sur R. A chacune

de ces mesures η, on associe la mesure de probabilité (ηMn ) définie par
Z
(ηMn ) : B ⊂ R 7→ (ηMn )(B) = η(dx) Mn (x, B) ∈ [0, 1]

Dans ce contexte, une mesure de probabilité η correspond à un opérateur intégral

sur l’ensemble des fonctions f mesurables et bornées

Z
η[f ] = η(dx) f (x) ∈ R

Avec ce système de notations, pour les fonctions indicatrices f = 1B d’ensembles B ⊂ R,

nous avons η(B) = η[1B ]. Lorsqu’il n’y a pas de confusions, il est coutume de noter
η(f ) et Mn (f ) plutôt que η[f ] et Mn [f ]. Cet abus de notation évident est parfois poussé
à l’extrême. Certains auteurs notent parfois tout simplement ηf , pour insister sur le
fait que cette opération intégrale n’est autre qu’une extension naturelle du produit
matriciel.
À la différence des chaı̂nes discrètes, les lois ηn (dx) des états Xn de chaı̂nes abstraites
sont données par des équations d’évolution intégrales. Il est bien entendu hors de
question d’être tenté de calculer, ou d’estimer ces formules de transport. Il convient
néanmoins de souligner qu’il existe des stratégies numériques et probabilistes pour
le faire ! Ces techniques sont connues sous le nom de méthodes de Monte-Carlo, en
référence au fait qu’elles sont basées sur des simulations concrètes de trajectoires
aléatoires.
 En utilisant la propriété de Markov, on notera que pour toute fonction fn+1 ,
mesurable et bornée sur En+1 , on a

E(fn+1 (Xn+1 )|Xn−1 ) = E( E(fn+1 (Xn+1 )|Xn−1 , Xn ) |Xn−1 )

= E( E(fn+1 (Xn+1 )|Xn ) |Xn−1 )
= E( Mn+1 (fn+1 )(Xn ) |Xn−1 )
= Mn [Mn+1 (fn+1 ])(Xn−1 )

On note (Mn Mn+1 )(xn−1 , dxn+1 ) la collection de mesures sur En+1 , indexée par les
xn−1 ∈ En−1 et définies pour tout Bn+1 ⊂ En+1 par la formule de composition intégrale
Z
(Mn Mn+1 )(xn−1 , Bn+1 ) = Mn (xn−1 , dxn ) Mn+1 (xn , Bn+1 )
En

En terme d’indicatrices cette équation s’exprime sous la forme suivante

Z
(Mn Mn+1 )(1Bn+1 )(xn−1 ) = Mn (xn−1 , dxn ) Mn+1 (1Bn+1 )(xn )
En

Plus généralement, on a la formule de conditionnement

E(fn+p (Xn+p )|Xn−1 ) = Mn Mn+1 . . . Mn+p (fn+p )(Xn−1 )

pour toute fonction fn+p bornée sur En+p , et pour tout décalage d’indice temporel
p ≥ 1. Les opérateurs intégraux

Mn,n+p =déf. Mn Mn+1 . . . Mn+p

avec n ≥ 1 et p ≥ 1, forment un semigroupe d’opérateurs intégraux, en ce sens où

∀n ≤ m ≤ n + p Mn,n+p = Mn,m Mm,n+p

4.1.2 Équations de Chapman-Kolmogorov

D’après la formule des conditionnements emboı̂tes, nous avons

ηn (fn ) = E(E(fn (Xn )|Xn−1 )) = E(Mn (fn )(Xn−1 ))

Z
= ηn−1 (dxn−1 ) Mn (fn )(xn−1 ) = ηn−1 (Mn (fn )) (4.2)
En−1

Soit (ηn−1 Mn ) la mesure sur En définie, pour tout Bn ⊂ En , par la formule

Z
(ηn−1 Mn )(1Bn ) = ηn−1 (dxn−1 ) Mn (1Bn )(xn−1 )
En−1
Par construction, nous avons les représentations équivalentes suivantes
ηn (fn ) = ηn−1 (Mn (fn ))
Z Z 
= ηn−1 (dxn−1 ) Mn (xn−1 , dxn ) fn (xn )
En−1 En
Z "Z #
= ηn−1 (dxn−1 ) Mn (xn−1 , dxn ) fn (xn )
En En−1
= (ηn−1 Mn )(fn ) =déf. ηn−1 Mn (fn )

En utilisant (4.2), nous obtenons la formule de transport intégral des lois des
états de la chaı̂ne

ηn = ηn−1 Mn = η0 M1 M2 . . . Mn (4.3)

Cette équation intégrale, appelée la formule de Chapman-Kolmogorov, permet

de voir les lois ηn comme solution d’un système dynamique (déterministe) intégral (et
donc linéaire) sur les espaces de mesures de probabilités.

4.1.3 Processus historique

Comme nous l’avons vu pour les chaı̂nes discrètes dans la section 3.2.2, le cadre non
homogène est utile pour représenter des modèles trajectoriels, tel le processus historique
associé à une chaı̂ne de Markov élémentaire. La construction abstraite de ces modèles
trajectoriels est analogue à celle présentée à l’exemple 3.2.2. Ainsi, si Xn0 est une chaı̂ne
de Markov de transitions Mn0 sur des ensembles En0 , les séquences de trajectoires
Xn = (X00 , . . . , Xn0 )
forment un processus de Markov sur les espaces produits
En = (E00 × . . . × En0 )
Le passage de Xn à Xn+1 s’effectue en deux temps. On conserve tout d’abord le segment
de trajectoire Xn = (X00 , . . . , Xn0 ), puis on lui adjoint une extension élémentaire aléatoire
0
Xn+1 0
de loi Mn+1 (Xn0 , dx0 ). Plus formellement, nous avons

Xn = (X00 , . . . , Xn0 ) Xn+1 = ((X00 , . . . , Xn0 ), Xn+1

0 )
| {z }
=( Xn 0
, Xn+1 ) ∈ En+1 = (En × E 0 )
On notera que pour toute fonction fn+1 bornée sur En+1 , et pour tout segment de
trajectoire xn = (x00 , . . . , x0n ) ∈ En , nous avons
E(fn+1 (Xn+1 ) | Xn = xn )

= E(fn+1 ([X00 , . . . , Xn ], Xn+1

0 ) | (X00 , . . . , Xn0 ) = (x00 , . . . , x0n ))

fn+1 ([x00 , . . . , x0n ], x0n+1 ) Mn+1

0 (x0n , dx0n+1 )
R
= x0n+1 ∈E 0

Cette équation s’exprime aussi sous la forme trajectorielle suivante

E(fn+1 (Xn+1 ) | Xn = xn )
Z
= fn+1 (y00 , . . . , yn+1
0
) δ(x00 ,...,x0n ) (d(y00 , . . . , yn0 ))Mn+1
0
(yn0 , dyn+1
0
)
(y00 ,...,yn+1
0 )∈En+1

La dernière formule montre que Xn est une chaı̂ne de Markov, de probabilités de

transitions Mn+1 de En vers En+1 , données par la formule

Mn+1 ((x00 , . . . , x0n ), d(y00 , . . . , yn+1

0 ))

= δ(x00 ,...,x0n ) (d(y00 , . . . , yn0 )) Mn+1

0 (yn0 , dyn+1
0 )

Ces processus historiques interviennent de façon naturelle dans divers problèmes

issus de la physique, ou de la biologie. Ils offrent un cadre markovien naturel pour
modéliser et analyser des évolutions aléatoires complexes, liées le plus souvent à des
effets de dépendance trajectorielles. Ainsi, dans la section 7.3, ces processus historiques
nous permettrons de définir des modèles d’arbres généalogiques en terme d’algorithmes
génétiques trajectoriels. Dans la section 7.4, nous utiliserons à nouveau ces modèles
pour représenter des explorations évolutionnaires basées sur des mécanismes de mémoire
renforçant les probabilités de retours vers des sites qui ont déjà été visités.

4.2 Chaı̂nes linéaires et gaussiennes

La distribution gaussienne joue un rôle essentiel dans la théorie des probabilités.
L’importance des variables gaussiennes est en grande partie due au théorème
central de la limite. Ce dernier nous informe que toute accumulation de petites
fluctuations indépendantes, et de nature quelconque, se traduit asymptotiquement
et irrémédiablement par une variable gaussienne. Tout phénomène résultant d’une
addition d’effets aléatoires indépendants est donc nécessairement de nature gaussienne.
Dans de nombreux problèmes pratiques, il est donc naturel de considèrer comme
gaussiennes les erreurs de mesures, et les erreurs de modélisation.
4.2.1 Formulation canonique
Les systèmes linéaires et gaussiens X = (Xn )n≥0 sont définis par la donnée d’une
suite de v.a. à valeurs réelles de lois marginales
(xp −ap xp−1 )2 x2
 
n − 2 − 02
2σ 2σ
Y e p
 e 0
P(X0 ,...,Xn ) (d(x0 , . . . , xn )) =  q dxp  p 2
dx0 (4.4)
p=1 2πσ 2
p
2πσ 0

Dans le formule ci-dessus (an )n≥1 désigne une suite de nombres réels, et (σn )n≥0 une
suite de nombres strictement positifs.
En intégrant la formule (4.4) en la coordonnée xn , on montre facilement que

PXn |(X0 ,...,Xn−1 ) (dxn |x0 , . . . , xn−1 )

= PXn |Xn−1 (dxn |xn−1 )

= Mn (xn−1 , dxn ) =déf. √ 1 2

exp {− 2σ12 (xn − an xn−1 )2 } dxn
2πσn n

pour tout n ≥ 1, et x0 , . . . , xn−1 ∈ R. De même, on montre que la loi de la condition

initiale X0 est donnée par

1 x20
PX0 (dx0 ) = η0 (dx0 ) =déf. p exp {− } dx0
2πσ02 2σ02

La formule multiplicative (4.4) doit donc se lire comme une formule de Bayes
séquentielle
 
n
Y
P(X0 ,...,Xn ) (d(x0 , . . . , xn )) =  PXp |Xp−1 (dxp |xp−1 ) PX0 (dx0 )
p=1

= η0 (dx0 ) M1 (x0 , dx1 ) . . . Mn (xn−1 , dxn )

4.2.2 Formulation dynamique

Notons W = (Wn )n≥0 une suite de v.a. indépendantes et gaussiennes, et telles que
n n
!
w 2
Y Y 1 p
P(W0 ,...,Wn ) (d(w0 , . . . , wn )) = PWp (dwp ) = √ exp {− } dwp
2π 2
p=0 p=0

On associe à cette séquence de v.a. gaussiennes, le système dynamique


Xn = an Xn−1 + σn Wn
(4.5)
X0 = σ0 W0
Par construction, nous avons pour toute fonction bornée f
E(f (Xn )|Xn−1 = xn−1 )

= E(f (xn−1 + σn Wn )|Xn−1 = xn−1 )

wp2
R +∞ −
e√ 2
= −∞ f (xn−1 + σn wn ) 2π
dwp

1
R +∞ − (xn −an xn−1 )2
f (xn ) √ 1 2
R
= −∞ 2
e 2σn dxn = Mn (xn−1 , dxn ) f (xn )
2πσn

Il est donc équivalent de définir la chaı̂ne Xn soit de façon dynamique, soit directement
par la donnée de ses transitions de probabilités.

4.3 Processus de Poisson

Le processus de Poisson est souvent associé à des phénomènes de comptages dans le
temps : arrivées de clients dans une file d’attente, nombres de transactions journalières
autour d’une action boursière, arrivées et départs d’avions dans des aéroports, nombres
d’appels dans un central téléphonique, etc. Nous renvoyons le lecteur à l’ouvrage de F.A.
Haight [?] consacré aux différentes applications de ce processus.
D’un point de vue mathématique, le processus de Poisson est défini en terme d’une
suite de v.a. (Tn )n≥0 , à valeurs positives, et de lois marginales données par la formule
P(T0 ,...,Tn ) (d(t0 , . . . , tn )) (4.6)

= [1[0,∞) (t0 ) λ e−λt0 dt0 ][1[t0 ,∞) (t1 ) λ e−λ(t1 −t0 ) dt1 ]

. . . × [1[tn−1 ,∞) (tn ) λ e−λ(tn −tn−1 ) dtn ]

Dans la formule ci-dessus, λ > 0 désigne un paramètre fixé. En interprétant (4.6) comme
une formule de Bayes séquentielle, on obtient
 
n
Y
P(T0 ,...,Tn ) (d(t0 , . . . , tn )) =  PTp |Tp−1 (dtp |tp−1 ) PT0 (dt0 )
p=1

avec
PT0 (dt0 ) = 1[0,∞) (t0 ) λ e−λt0 dt0
PTn |Tn−1 (dtn |tn−1 ) = 1[tn−1 ,∞) (tn ) λ e−λ(tn −tn−1 ) dtn
Par un simple changement de variable, on montre aussi facilement que
U0 =déf. T0 et Un =déf. (Tn − Tn−1 )
forment une suite de v.a. exponentielles et indépendantes de paramètre λ. Le processus
continu de comptage des sauts
X
t ∈ R+ 7→ N (t) = 1[Tn ,∞) (t) ∈ N
n≥0

est appelé le Processus de Poisson d’intensité λ. Une réalisation de ce processus est

donnée dans la figure suivante.

N(t)=4

t temps
T0 T1 T2 T3 T4

Fig. 4.1 – Processus de Poisson

4.4 Évolutions dans des milieux absorbants

En termes physiques, la chaı̂ne de Markov suivante représente l’évolution d’une
particule physique dans un puit de potentiel absorbant. Ce modèle peut aussi
s’interpréter comme une désintégration, ou encore comme une absorption de la
radioactivité dans des containers de stockages de déchets nucléaires.
Très brièvement, une particule évolue aléatoirement dans un environnement
absorbant avec un taux de survie G(x) ∈ [0, 1], en chaque site x.
Les obstacles sont d’autant plus absorbants que les valeurs de G sont proches de 0.
Inversement, la particule évolue librement dans des régions où G vaut 1. Les régions
où le potentiel est strictement inférieur à 1 représentent des obstacles dans lesquels la
particules peut être piégée.

Ces trappes peuvent représenter des niveaux de sécurité dans tout type
d’environnements, tels des aéroports, des milieux carcéraux, des chaı̂nes de
production, des réseaux de télécommunications, des containers de stockage de
radioactivité, des tissus cellulaires, etc.
 En pharmacologie, ces modèles d’aborption sont aussi utilisé pour modéliser
les évolution des taux de leukocytes dans des traitement du cancer. Dans ce
contexte, les trappes reflètent des chutes de niveaux de leukocytes dans lesquels
un individu risque de decéder.

En biologie et en recherche médicale, ces modèles d’aborption sont aussi

utilisé pour modéliser l’évolution de photons émis par un laser sur un tissus
cellulaire. Les trappes représentent des regions cellulaires, telles des tumeurs
plus sombres absorbant les photons. L’analyse de ces modèles permet de
localiser et d’analyser ces régions de gonflement pathologiques des tissus.

Enfin, ces modèles peuvent aussi s’interpreter comme des processus

économiques tels ldes évolutions de portefeuilles, ou tout autre indicateur
économique, dans des milieux financiers ou géopolitiques-politiques.

Dans tous ces domaines d’applications, l’un des problèmes majeurs est de calculer
les probabilités de défaillances, autrement dit les probabilités pour que l’aborption de
la particule ne soit pas éffective à des instants donnés.
Pour fixer les idées, nous conviendrons par la suite que l’espace des états est donné
par le réseau Zd . Plus formellement, on adjoint à l’espace des sites Zd , un point cimetière
“c”. On considère alors une particule évoluant sur l’espace augmenté E = Zd ∪ {c},
selon une transition Yn Yn+1 se décomposant en deux étapes
absorption évolution
Yn −−−−−−−−→ Ybn −−−−−−−−→ Yn+1
Pour construire le mécanisme d’absorption, on considère à chaque instant n, une
une collection de v.a. indépendantes (n (x))x∈Zd , de Bernoulli à valeurs dans {0, 1},
de paramètre G(x). Sur l’évènement n (x) = 0, le site x devient une pure trappe
hautement absorbante. Inversement, si n (x) = 1, le site x reste viable, et la particule
peut le traverser sans encombres.
Le mécanisme d’évolution libre est analogue à celui de la marche aléatoire sur Zd
présenté dans l’exemple 3.3.4. Pour le décrire, on se donne une mesure de probabilité
p, sur l’ensemble des 2d vecteurs unitaires directionnels
U = {u ∈ Zd : |u| = 1}
ainsi qu’une suite de vecteurs indépendants (Un )n≥1 de même loi p sur U. On suppose
que les suites de v.a. d’exploration Un , et de v.a. d’absorption n (x), sont indépendantes.
On convient enfin que n (c) = 0 = g(c).
 Les transitions élémentaires de la chaı̂ne sont définies récursivement de la façon
suivante. Supposons que la particule se trouve à l’instant n sur un site Yn = x ∈ Zd (si
Yn = c, on pose Yn+1 = c).
1. Avec une probabilité (1 − G(x)), la particule est tuée, puis placée dans l’état
cimetière. Dans ce cas, on pose Ybn = c. Dans le cas contraire, la particule reste
active et l’on pose Ybn = x. Plus formellement, on a

Ybn = n (Yn ) Yn + (1 − n (Yn )) c

2. Comme il n’y a semble-t-il pas de vie après la mort, lorsque la particule a été tuée,
elle reste inactive à l’instant suivant. Dans ce cas, on pose Yn+1 = c. Dans le cas
contraire, la particule est encore active, et effectue un mouvement exploratoire
analogue à celui d’une marche aléatoire sur Zd . Plus formellement, on pose dans
ce dernier cas
Yn+1 = Ybn + Un+1

Par construction, la suite

Y0 → Yb0 → . . . → Yn → Ybn → Yn+1 → Ybn+1 → . . .

forme une chaı̂ne de Markov de transitions

P(Ybn = y | Yn = x) = G(x) 1x (y) + (1 − G(x)) 1c (y)

P(Yn+1 = z | Ybn = y) = 1c (y) 1c (z) + 1Zd (y) K(y, z)

avec la probabilité de transition

X
K(y, z) = p(u) 1y+u (z)
u∈U

La figure suivante montre deux réalisations de trajectoires absorbées dans des

“poches”d’obstacles associés à des régions où le potentiel G(x) < 1. Lorsque G(x) = 1,
la particule ne subit pas le mécanisme d’absorption, et évolue librement comme une
marche aléatoire sur Z2 . Les traits pointillés témoignent du fait qu’une particule
visitant une poche d’obstacles s’essouffle. Sa durée de vie diminue à chaque instant,
en “subissant” les v.a. trappes de Bernoulli.

Exercice 4.4.1 On note T l’instant d’absorption de la particule

T = inf {n ≥ 0 : Ybn = c} = inf {n ≥ 0 : n (Yn )n = c}

 g(x)<1
(c)
0

(c)
g(x)=1

Fig. 4.2 – Particule absorbée dans Z2

1. Montrer que  
n
Y
P(T > n) = E  Gp (Xp )
p=0

où Xn désigne une chaı̂ne de Markov sur Zd de transitions de probabilités K,

et de même loi que Y0 sur Zd . Lorsque le potentiel est uniformément majoré
G(x) ≤ e−λ , avec λ > 0, montrer que

P(T > n) ≤ e−λ(n+1)

Interpréter ce résultat.
2. Vérifier que la loi d’une particule non aborbée est donnée pour toute fonction
bornée f sur Zd par la formule renormalisée de Feynman-Kac
 
E f (Xn ) np=0 Gp (Xp )
Q
E(f (Yn ) | T > n) = Q 
n
E p=0 G p (Xp )
Chapitre 5

Chaı̂nes de Markov non linéaires

5.1 Introduction
Les modèles de chaı̂nes de Markov non linéaires présentés dans cette section sont
une extension naturelle des modèles markoviens étudiés dans la première partie de ce
chapitre.
Les premiers modèles de ce type sont semble-t-il apparus dans la littérature du
traitetement du signal, et plus particulièrement en filtrage non linéaire ([?, ?]). L’idée
de départ est la suivante. Les lois conditionnelles des états Xn du signal par rapports
aux observations reçues (Yp )0≤p<n sont données par un flot de mesures

η0 = Loi(X0 ) −→ η1 = Loi(X1 | Y0 ) −→ . . . −→ ηn = Loi(Xn | Y0 , . . . , Yn−1 ) −→ . . .

Ces mesures sont appelées les prédicteurs optimaux. Elles peuvent être calculées de
façon récursives suivant une équation de la forme

ηn = Φn (ηn−1 )

où Φn désigne une transformation plus ou moins complexe sur l’espace des distributions
sur l’espace d’état du signal. Bien évidemment, ces équations sont en général impossible
à résoudre explicitement. Pour utiliser des méthodes de simulation de type Monte Carlo,
l’idée naturelle est d’interpréter ces mesures ηn comme les lois d’une chaı̂ne de Markov
X n sur l’espace d’états du signal :

ηn = Loi(X n )

Ces chaı̂ne de Markov X n ne sont pas uniques. On peut choisir par exemple des suite
de variables indépendentes X n de lois

ηn = Φn (ηn−1 ) avec ηn−1 = Loi(X n−1 )

91
On peut aussi chercher à écrire les transformations Φn comme un transport de mesures
markovien
ηn = Φn (ηn−1 ) = ηn−1 Kn,ηn−1
Dans la formule précédente, Kn,η (xn−1 , dxn ) désigne une famille de probabilités de
transitions de En−1 vers un espace En , indexées par le paramètre temporel et par les
mesures η sur En−1 . Dans ce contexte, nous avons à nouveau

ηn = Loi(X n )

où (X n )n≥0 désigne une une chaı̂ne de Markov de loi initiale η0 = Loi(X 0 ) sur E0 et de
transitions de probabilités élémentaires d’un espace En−1 vers un espace En décrites
par :


P X n ∈ dxn | X n−1 = xn−1 = Kn,ηn−1 (xn−1 , dxn ) avec ηn−1 = Loi(X n−1 )

Nous étudierons ces modèles mathématiques avec plus de détails dans la section 5.2.
Il me semble est important de souligner que ces modèles discrets sont très proches
des équations à temps continu issues de la mécanique des fluides, telles les équation
diffusives de type McKean-Vlasov et/ou les modèles de gaz de type Boltzmann.

5.2 Description des modèles

La plus simple façon de les définir une chaı̂ne de Markov non linéaires est de
se donner un chaı̂ne de Markov (X n )n≥0 de loi initiale η0 = Loi(X 0 ) sur E0 et de
transitions de probabilités élémentaires d’un espace En−1 vers un espace En décrites
par une formule de la forme suivante :


P X n ∈ dxn | X n−1 = xn−1 = Kn,ηn−1 (xn−1 , dxn ) avec ηn−1 = Loi(X n−1 )

Dans la définition précédente, Kn,η (xn−1 , dxn ) désigne une famille de probabilités de
transitions de En−1 vers un espace En , indexées par le paramètre temporel et par les
mesures η sur En−1 .
Il est assez aisé de vérifier que l’on a


P (X 0 , . . . , X n ) ∈ d(x0 , . . . , xn ) = η0 (dx0 )K1,ηn−1 (x0 , dx1 ) . . . Kn,ηn−1 (xn−1 , dxn )

Ces mesures sont appelés les mesures de McKean associés aux transitions de probabilités
(Kn,η )n,η . Une simple intégration par rapports aux coordonnées temporelles permet de
vérifier que
Loi(X n ) = η0 K1,η0 . . . Kn,ηn−1 = ηn−1 Kn,ηn−1
Une simple récurrence permet de s’assurer que

ηn = Law(X n ) = ηn−1 Kn,ηn−1

5.3 Interprétations particulaires en champ moyen
Supposons dans un premier temps que les transitions de la chaı̂ne

X n−1 Xn

sont faciles à simuler. Dans ces conditions, pour approcher les lois ηn de la chaı̂ne à
i
chaque instant, il suffit de simuler N copies indépendantes (X n )1≤i≤N de la chaı̂ne X n .
Ce schéma de simulation numérique est illustré par la figure suivante
Kn,ηn−1 1
1 Xn
X n−1 −−−−−−−−−−→ ..
.. .
. ..
i
Kn,ηn−1 .
X n−1 −−−−−−−−−−→ i
.. Xn
. ..
Kn,ηn−1
.
N N
X n−1 −−−−−−−−−−→ X n

Une simple application de la loi des grands nombres nous donne l’approximation
suivante :
N
1 X
∀n ≥ 0 ηnN := δXni 'N ↑∞ ηn
N
i=1

Malheureusement, dans la plupart des cas les lois de ces chaı̂nes non linéaires ηn
n’ont aucune expression analytique explicite et/ou ne peuvent être simulées de façon
exacte par aucun algorithme de simulation dans un temps raisonnable.
L’idée des méthodes particulaires de type champ moyen est d’utiliser la population
courante pour approcher ces lois complexes. Sous certaines conditions de régularités
sur les transitions Kn,η , nous avons
N
ηn−1 'N ↑∞ ηn−1 −→ Kn,ηN 'N ↑∞ Kn,ηn−1
n−1

En utilisant ces formules d’approximation, on définit récursiment une chaine de

Markov
(ξn(1,N ) , ξn(2,N ) , . . . , ξn(N,N ) )n≥0
(1,N ) (2,N ) (N,N )
sur les espaces produits EnN . L’état initial (ξ0 , ξ0 , . . . , ξ0 ) est donné par N
variables aléatoires indépendantes et identiquement distribuées de loi η0 . Les transitions
élémentaires de cette chaı̂ne sont données par le schéma suivant
 Kn,ηN
n−1 (1,N )
(1,N )
ξn−1 −−−−−−−−−−→ ξn
.. ..
. .
Kn,ηN
..
(i,N ) n−1 .
ξn−1 −−−−−−−−−−→ (i,N )
ξn
.. ..
. .
Kn,ηN
n−1 (N,N )
−−−−−−−−−−→ ξn
(N,N )
ξn−1
Plus formellement, nous avons

  N
(N )
Y
1 N
P ξn+1 ∈ d(x , . . . , x ) | ξn(N ) = Kn+1,ηnN (ξn(N,i) , dxi ) (5.1)
i=1

avec les mesures d’occupation du système

N
1 X
ηnN := δξ(N,j)
N n
j=1

Autrement dit, connaissant une réalisation du système au temps n, la population des

individus au temps (n + 1)
 
(N,1) (N,2) (N,N )
ξn+1 , ξn+1 , . . . , ξn+1

est formée de N variables aléatoires indépendantes de lois respectives

Kn+1,ηnN (ξn(N,1) , dx1 ) , Kn+1,ηnN (ξn(N,2) , dx2 ) , . . . Kn+1,ηnN (ξn(N,N ) , dxN )

Sous certaines hypothèses de régularité, on peut montrer que l’on a en un certain sens
N
1 X
ηnN := δξ(N,j) 'N ↑∞ ηn
N n
j=1

5.4 Champs moyens de type gaussien

On considère l’équation d’évolution suivante sur l’ensemble P(Rd ) des mesures de
probabilités sur l’espace d-dimensionnel Rd :
Z
ηn+1 (dy) = (ηn Kn+1,ηn ) (dy) := ηn (dx) Kn+1,ηn (x, dy) (5.2)
Rd
Dans la définition précédente, les probabilités de transitions Kn,η sur En = Rd définis
par les mesures gaussiennes :
1 1
exp − (y − an (x, η))0 Q−1

Kn,η (x, dy) = d n (y − an (x, η)) dy
2
p
(2π) 2 |Qn |

où
– Qn désigne une matrice (d × d) symétrique et définie positive, |Qn | := det(Qn ).
– dy désigne la mesure de Lebesgue sur Rd , y = (y 1 , . . . , y d ) et x ∈ Rd .
– an : Rd × P(Rd ) → Rd est une fonctions suffisamment régulière et bornée.
Cet exemple est assez intéressant car les mesures de McKean associées

η0 (x0 )K1,η0 (x0 , dx1 ) . . . Kn,ηn−1 (xn−1 , dxn )

peuvent s’interpréter comme les lois des trajectoires (X 0 , . . . , X n ) d’une chaı̂ne de

Markov donnée par les équations cinétiques suivantes :

X n = an (X n−1 , ηn−1 ) + Wn , n≥1

avec
– Wn , n ≥ 1, une suite de variables aléatoires indépendantes, à valeurs dans Rd , et
de lois gaussienne
1 1
exp − y 0 Q−1

λn (dy) = d n y dy
2
p
(2π) 2 |Qn |

– X 0 une variable aléatoire de loi η0 , et indépendante de la suite Wn .

– pour chaque n ≥ 0, ηn est la loi de X n .
L’interprétation particulaire de type champ moyen (5.1) associée au modèle d’évolution
(5.2) est la chaı̂ne de Markov
 
ξn(N ) = ξn(N,1) , ξn(N,2) , . . . , ξn(N,N ) ∈ (Rd )N

définie par les équations cinétiques suivantes :

N

(N,i) N
 1 X
∀1 ≤ i ≤ N ξn(N,i) = an ξn−1 , ηn−1 + Wni avec N
ηn−1 := δξ(N,j)
N n−1
j=1

Dans cette formulation, la suite (Wni )1≤i≤N est formée de N copies indépendantes de
(N,i)
Wn . L’états initial (ξ0 )1≤i≤N est encore formé de N copies indépendantes de X 0 .
5.5 Modèles simplifiés de gaz de McKean
Dans cette section, nous présentons le modèle simplifié de gaz à deux vitesses de
McKean étudié dans l’article [?]. L’analyse qui suit peut s’étendre à des modèles plus
sophistiqués avec des vitesses multiples. Le modèle le plus élémentaire est défini sur
l’espace E = {−1, +1} en associant à chaque mesure η dans l’ensemble P({−1, +1})
des mesures de probabilités sur = {−1, +1} la transition suivante :

Kn,η (x, dy) = η(+1) δx (dy) + η(−1) δ−x (dy). (5.3)

On associe à ces transitions le processus à valeurs mesure suivant :

ηn+1 (+1) = ηn (+1) Kn+1,ηn (+1, +1) + ηn (−1) Kn+1,ηn (−1, +1)
= ηn (+1)2 + ηn (−1)2 = ηn (+1)2 + (1 − ηn (+1))2

On pourra remarquer que l’on a

ηn+1 (−1) = ηn (+1) Kn+1,ηn (+1, −1) + ηn (−1) Kn+1,ηn (−1, −1)
= 2 ηn (+1)ηn (−1)

L’interprétation particulaire (5.1) associée c̀es équations est la chaı̂ne de Markov

 
ξn(N ) = ξn(N,1) , ξn(N,2) , . . . , ξn(N,N ) ∈ {−1, +1}N

définie par la formule

(N,i)
∀1 ≤ i ≤ N ξn(N,i) = (N,i)
n ξn−1
(N,i)
Dans la définition précédente (n )1≤i≤N représente une suite de variables aléatoires
conditionnellement indépendantes sur {−1, +1} de loi :
N

(N )
 1 X
P (N,i)
n =  | ξ N N
n−1 = ηn−1 (+1) 11 () + ηn−1 (−1) 1−1 () avec N
ηn−1 := δξ(N,j)
N n−1
j=1

5.6 Flots de mesures de Feynman-Kac

5.6.1 Description des modèles
On se donne une suite de fonctions positives et bornées Gn sur En . On note PGn (En )
l’ensemble des mesures de probabilités µ sur En telles que µ(Gn ) 6= 0 On associe à
chaque fonction potentiel Gn une transformation de Boltzamnn-Gibbs ΨGn définie par
les formules suivantes

ΨGn : µ ∈ PGn (E) 7→ ΨGn (µ) ∈ PGn (E)

avec la mesure de probabilité ΨGn (µ) donnée par
1
ΨGn (µ)(dx) := Gn (x) µ(dx)
µ(Gn )
On se donne Mn (xn−1 , dxn ) une suite de probabilités de transitions de En−1 dans En ,
avec n ≥ 1. On considère enfin une mesure de probabilité η0 sur E0 telle que η0 (G0 ) > 0.
On notera par la suite Xn une chaı̂ne de Markov à valeurs dans En de loi initiale η0 et
de probabilités de transitions

Mn (xn−1 , dxn ) := P(Xn ∈ dxn | Xn−1 = xn−1 )

On convient que
∀b
ηn ∈ PGn (En ) ηbn Mn+1 ∈ PGn+1 (En+1 )
autrement dit
ηbn (Gn ) > 0 =⇒ ηbn Mn+1 (Gn ) > 0
Lorsque les fonctions potentiels sont strictement positives, les espaces PGn (En )
coincident avec l’espace de toutes les mesures de probabilités sur En . Dans ce cas,
les conditions précédentes sont trivialement satisfaites.
On considère (ηn )n≥0 le flot de mesures de probabilités définit par les équations :

ηn = Φn (ηn−1 ) := ΨGn−1 (ηn−1 )Mn (5.4)

Ces mesures ηn et leurs mesures mises à jour ou corrigées

ηbn := ΨGn (ηn )

peuvent s’exprimer sous forme d’intégrales de chemins pondérés connus sous le nom
de formules de Feynman-Kac. Ces représentations fonctionnelles sont données sur des
fonctions tests fn par les formules suivantes

ηn (fn ) = γn (fn )/γn (1) et ηbn (fn ) = γ

bn (fn )/b
γn (1) (5.5)

avec des mesures non normalisées γn et γ

bn décrites ci-dessous :
Q
γn (fn ) = E[fn (Xn ) 0≤k<n Gk (Xk )] et γ bn (fn ) = γn (fn Gn ) (5.6)

Pour vérifier cette formule, on utilise tout d’abord la propriété de Markov pour
vérifier que
 
Y
γn (fn ) = E E (fn (Xn ) | (Xp )0≤p<n ) Gp (Xp )
0≤p<n
 
Y
= E E (fn (Xn ) | Xn−1 ) Gp (Xp )
0≤p<n
Ceci entraine que
 
Y
γn (fn ) = E Mn (fn )(Xn−1 ) Gp (Xp )
0≤p<n

avec
Mn (fn )(xn−1 ) := E (fn (Xn ) | Xn−1 = xn−1 )
Dans un second temps, on constate que l’on a :
 
Y
γn (fn ) = E Gn−1 (Xn−1 )Mn (fn )(Xn−1 ) Gp (Xp )
0≤p<(n−1)

Par définition de la mesure γn−1 , nous avons

γn (fn ) = γn−1 (Gn−1 Mn (fn )) and γn (1) = γn−1 (Gn−1 )

Ceci entraı̂ne que

γn (fn ) γn−1 (Gn−1 Mn (fn ))
ηn (fn ) = =
γn (1) γn−1 (Gn−1 )
et par conséquent

γn−1 (Gn−1 Mn (fn ))/γn−1 (1) ηn−1 (Gn−1 Mn (fn ))

ηn (fn ) = =
γn−1 (Gn−1 )/γn−1 (1) ηn−1 (Gn−1 )

On en conclut que

∀fn ∈ Bb (En ) ηn (fn ) = ΨGn−1 (ηn−1 )(Mn (fn ))

m
ηn = ΨGn−1 (ηn−1 )Mn
Les mêmes arguments permettent d’analyer les mesures (b γn , ηbn ).
Ces mesures de probabilités permettent de modèliser une variété considérable de
problèmes issus de la physique, ou de la biologie : traitement du signal non linéaire,
description de macro-polymères et de chaı̂nes auto-évitantes, analyse d’évènements
rares, représentation de valeurs propres et d’états fondamentaux d’opérateurs de
Schrödinger,... Nous examinerons un certain nombre de ces questions dans le
chapitre ??, pour plus de détails nous renvoyons le lecteur aux ouvrages [?], et [?]. Dans
la section ??, nous présentons un algorithme de simulation universel de ces mesures
de Feynman-Kac. Ces modèles particulaires sont fondés sur l’évolution d’individus en
interaction explorant l’espace selon des mécanismes de mutation et sélection de type
génétique.
5.6.2 Chaı̂nes de Markov non linéaires
Commençons par remarquer qu’une transformation de Boltzmann-Gibbs
1
ΨG (µ)(dx) := G(x) µ(dx)
µ(G)

associée à une fonction potentiel G sur un espace d’état E, peut s’interpréter comme
un transport de mesure markovien. Pour être plus précis, on note (µ) une famille de
paramètre pouvant dépendre de G et µ, et telle que

(µ) G(x) ≤ 1 pour µ-presque tous les x ∈ E.

Avec ce système de notations, nous avons

 Z 
ΨG (µ) = µSµ ⇐⇒ ΨG (µ)(dy) = µ(dx)Sµ (x, dy)
E

avec la famille de probabilités de transition sur E données par la formule suivante :

Sµ (x, dy) = (µ) G(x) δx (dy) + (1 − (µ) G(x)) ΨG (µ)(dy).

Pour vérifier cette assertion, on se donne une fonction test bornée f sur E et l’on
observe que

Sµ (f )(x) = (µ) G(x) f (x) + (1 − (µ) G(x)) ΨG (µ)(f )

On en déduit les formules suivantes :

µ (Sµ (f )) = (µ) µ(Gf ) + (1 − (µ) µ(G)) ΨG (µ)(f )

µ(Gf )
= (µ) µ(Gf ) + ΨG (µ)(f ) − (µ) µ(G) = ΨG (µ)(f ).
µ(G)

Reprenons les équations non linéaires (5.4) vérifiées par les flots de mesures de
Feynman-Kac (ηn )n≥0 introduites en (5.5)

ηn+1 = Φn+1 (ηn ) := ΨGn (ηn )Mn+1

D’après les calculs précédents, nous avons l’équation de transport

ΨGn (ηn ) = ηn Sn,ηn

avec les transitions Sn,ηn (xn , dyn ) sur En définies par les formules suivantes

Sn,ηn (xn , dyn ) = n (ηn ) Gn (xn ) δxn (dyn ) + (1 − n (ηn ) Gn (xn )) Ψn (ηn )(dyn )
avec des paramètres n (ηn ) ≥ 0 tels que n (ηn )Gn (xn ) ≤ 1, pour tous les xn ∈ En . On
en conclut que
ηn+1 = ηn Kn+1,ηn
avec la famille de transitions de probabilités composées

Kn+1,ηn (xn−1 , dxn ) = Sn,ηn Mn+1 (xn−1 , dxn )

Z
= Sn,ηn (xn , dyn ) Mn+1 (yn , xn+1 )dxn+1 (5.7)
En

5.6.3 Champs moyens de type évolutionnaire

L’interprétation particulaire (5.1) associée aux équations (5.2) est la chaı̂ne de
Markov  
ξn(N ) = ξn(N,1) , ξn(N,2) , . . . , ξn(N,N ) ∈ (En )N

de transitions élémentaires données par la formule suivante :

  N   
(N )
Y
P ξn+1 ∈ dxn+1 | ξn(N ) = Sn,ηnN Mn+1 ξn(N,i) , dxin+1 (5.8)
i=1

avec les mesures d’occupation ηnN données par

N
1 X
ηnN := δξ(N,j)
N n
j=1

On notera que l’on a

  Z
(N )
P ξn+1 ∈ dxn+1 | ξn(N ) = Sn (ξn(N ) , dxn ) Mn+1 (xn , dxn+1 )
N
En

avec les transitions de Boltzmann-Gibbs Sn de EnN dans lui même et les transitions
d’exploration Mn+1 de EnN dans En+1
N définies par

N
Y
Sn (ξn(N ) , dxn ) = Sn,ηnN (ξn(N,i) , dxin )
i=1
YN
Mn+1 (xn , dxn+1 ) = Mn+1 (xin , dxin+1 )
i=1

Ces décompositions soulignent le fact que la chaı̂ne de Markov suit les deux mêmes
étapes de mise à jour/correction et d’exploration/mutation que le flot des mesures de
Feynman-Kac limites. Plus formellement, les deux étapes de correction/prediction dans
l’espace des distributions :
Sn,ηn Mn+1
ηn ∈ P(En ) −−−−−−−−→ ηbn = ηn Sn,ηn = ΨGn (ηn ) ∈ P(En ) −−−−−
−−→ ηn+1 = ηbn Mn+1
(5.9)
sont approchées par deux étapes évolutionnaires de type sélection/mutation dans
l’espace des mesures empiriques :
selection (N ) mutation
ξn(N ) ∈ EnN −−−−−−−−→ ξbn(N ) ∈ EnN −−−−− N
−−→ ξn+1 ∈ En+1 (5.10)

L’équation d’évolution complète du système peut se résumer par le diagramme

synthétique suivant :

(N,1)
 
(N,1) Mn+1 (N,1)

ξn ξb −−−−−−−−−−→ ξn+1
..   n. .. 
. 
 Sn,ηN

 .. . 

(N,i) n (N,1)
 −−−−−−−−−−→  ξbn
 (N,i)
ξn ξn+1
 
−−−−−−−−−−→ 
..  
.. .. 
. .
  
  . 
(N,N ) (N,N ) (N,N )
ξn ξn
b −−−−−−−−−−→ ξn+1

avec les probabilités d’acceptation-rejets :

(N,i)
Sn,ηnN (ξn , dx)

(N,i)
:= n (ηnN ) Gn (ξn ) δξ(N,i) (dx)+
n

  P (N,j)
(N,i) N Gn (ξn )
1 − n (ηnN )Gn (ξn ) j=1 PN (N,k) δξ (N,j) (dx)
k=1 Gn (ξn ) n

Le choix du paramètre d’acceptation n’est pas unique. L’algorithme genetique simple

de type mutation/selection correspond au choix n (ηnN ) = 0. Si on pose n (ηnN ) =
(N,i)
1/ supi Gn (ξn ), alors les élites de la population qui maximisent le potentiel sont
toujours acceptées. Lorsque kGn k est explicitement connu, on peut aussi poser n (ηnN ) =
1/kGn k.
Pour plus de détails concernant ces modèles avec notamment la description des
arbres généalogiques associées nous renvoyons le lecteur au chapitre 7, ainsi qu’à
l’ouvrage [?].
Chapitre 6

L’équation de la chaleur

6.1 Les fluctuations browniennes

Les physiciens se sont longtemps penché sur la définition de la chaleur. Depuis la
nuit des temps, ou presque, nous savons que la chaleur peut être générée par frottements
vigoureux entre objets. En est-il de même au niveau atomique ? La répartition de la
chaleur dans un corps serait elle dictée par les excitations et les évolutions chaotiques
des différents atomes constituant la matière ? Est il envisageable simuler l’agitation
thermique induite par des collisions de molécules monoatomiques dans des gaz, ou
dans des liquides ? Peut on tenir compte des différents degrés de de frottement induits
par la matière ?
Pour comprendre l’origine de ces questions, il faut remonter au début du 19ième
siècle, et plus précisément en 1827. Cette année là, le biologiste Robert Brown s’amuse à
observer les mouvements erratiques de différents grains de pollen dans une goutte d’eau.
Ces grains de matière subissent des chocs insaissants et imprévisibles avec les molécules
du liquide en agitation permanente. Leur trajectoires aléatoires totalement erratiques
et chaotiques, semblent néanmoins dictées par des lois physiques bien précises.
D’après la théorie stochastique moderne, la répartition de la chaleur est elle aussi
l’expression de la répartition spatiale de “grains de chaleur” animés par ces mouvements
browniens. Si l’on chauffe une barre de fer rectiligne en un point donné, les molécules de
matière en ce point prennent de l’énergie, et se mettent à vibrer, et à se cogner en elles.
Ces successions de collisions de chaque cotés confèrent aux particules des mouvements
imprévisibles et chaotiques, allant vers la droite ou vers la gauche. Lorsque le temps
s’écoule, les grains de chaleur se trouvent distribués sur la barre selon une distribution
gaussienne centrée en la source de chaleur, et de plus en plus étalée sur les bords.
Le mouvement brownien uni-dimensionnel peut être peut être construit à partir
d’une simple marche aléatoire sur la droite des réels R. Un marcheur virtuel partant
de l’origine choisit à chaque instant d’évoluer aléatoirement vers la droite ou vers la
gauche. Ses déplacements vers la gauche −∆X ou vers la droite +∆X, par unité de

103
temps ∆t sont extrèment amples et erratiques, en ce sens où

[∆X]2 = ∆t
√
Autrement dit, le marcheur éffectue des amplitudes de déplacement du type ± ∆t,
par unité de temps ∆t. Par exemple, pour des fréquences temporelles très rapides de
∆t = 10−10 secondes, les amplitudes spatiales sont données par, disons ∆X = 10−5
mètres. La vitesse de ce marcheur est extrèmement rapide. Il court tout simplement à
∆X/∆t = 105 mètres par secondes !
Supposons que notre marcheur évolue dans des échelles sub-atomiques. Pour un
observateur macroscopique cela revient à s’éloigner de plus en plus du mouvement
décrit plus haut en faisant tendre les deux pas de la grille spatio-temporelle ∆X et
∆t vers 0. Pour cet observateur, ces particules browniennes ont tout simplement une
amplitude de vitesse infinie
∆X ±1
=√ −→ ±∞ lorsque ∆X et ∆t ↓ 0
∆t ∆t
La construction probabiliste du mouvement brownien développée en 1923 par
l’américain Nobert Wiener est assez sophistiquée. Dans ce qui suit, nous allons essayer
d’en donner les grandes lignes. Par soucis de clarté, nous conviendrons que l’intervalle
de temps est simplement donné par le segment unité [0, 1]. On choisit une séquence
d’instants ti équi-distribués

t0 = 0 < t1 = 1/n < . . . < ti = i/n < . . . < tn−1 = (n − 1)/n < tn = 1

avec
ti − ti−1 = ∆t = 1/n
On construit sur cette subdivision, une marche aléatoire simple (βtni )i=0,...,n centrée en
l’origine : √
β0n = 0 et ∆βtni = βtni − βtni−1 = ti ∆t
Les variables aléatoires ti représentent les mouvements aléatoires indépedants, vers la
gauche ou bien vers la droite, de notre marcheur indécis
1
P(ti = +1) = P(ti = −1) =
2
Par construction, cette marche aléatoire possède les deux propriétés essentielles
suivantes

∆βtni × ∆βtni = ∆t
E(∆βtni | βtn0 , . . . , βtni−1 ) = 0
Pour vérifier la seconde assertion, on pourra noter que l’on a
√ √
 
n n n 1 1
E(∆βti | βt0 , . . . , βti−1 ) = E(ti ) ∆t = − ∆t = 0
2 2
La première propriété concerne la vitesse du marcheur, la seconde souligne le fait que
ses déplacements aléatoire et locaux sont nuls, en moyenne.
Entre chaque instant ti−1 et ti , on convient que le marcheur évolue de façon
rectiligne, sur la pente entre βtni−1 et βtni . Autrement dit, sa position aux temps
t ∈ [ti−1 , ti [ est donnée par l’interpolation linéaire des deux extrémités βtni−1 et βtni .
Plus formellement, nous avons
βtn −βtn
βtn = βtni−1 + i i−1
ti −ti−1 (t − ti )

Par définition de la séquence d’instants ti , nous avons ti − ti−1 = 1/n, et donc

√
βtn = βtni−1 + (nt − nti−1 ) ∆βtni = βtni−1 + n (t − [nt]/n) ti

Après avoir noter que

 
i−1 i [nt] (i − 1)
t ∈ [ti−1 , ti [= , ⇐⇒ nt ∈ [i − 1, i[ ⇐⇒ = = ti−1
n n n n
on obtient pour chaque t ∈ [ti−1 , ti [, la décomposition suivante

βtn = = βtni−1 + (βtn − βtni−1 )

i−1
X √
= ∆βtnj + n (t − [nt]/n) t([nt]+1)/n
j=1
[nt]
X √
= ∆βtnj + n (t − [nt]/n) t([nt]+1)/n
j=1

Pour poursuivre notre discussion, il convient de rapeller que les accroissements de la

marche aléatoire ∆βtnj forment une suite de variables aléatoires donnés par la formule
√ 1
∆βtnj = ti ∆t = √ ti
n
On a de plus les majorations
√
2
0 ≤ (t − [nt]/n) ≤ 1/n ↓ 0 et [βtn − βtni−1 ]2 = n (t − [nt]/n) ≤ 1/n ↓ 0

On en conclut que
 
[nt] r [nt]
1 X [nt]  1 X
βtn = √ × ti + o(1) = p × ti  + o(1)
n n [nt]
i=1 i=1
où o(1) désigne une fonction aléatoire du paramètre n qui tend presque surement vers
0, lorsque n tends vers l’infini.
P[nt]
Le terme somme i=1 ti est formé de [nt] variables aléatoires (epsilonti )i=1,...[nt]
indépendantes et de même loi. Le théorème central de la limite nous assure que la suite
[nt]
1 X
p × ti
[nt] i=1

converge faiblement, lorsque n tend vers l’infini, vers une variable aléatoire gaussienne,
de moyenne nulle et de variance unité. q √
Il reste à noter que le facteur déterministe [nt]
n converge vers t. Par des arguments
probabilistes élémentaires (le lemme de Slutky), on en conclut que

βtn −→ βt
√
où βt désigne une variable aléatoire gaussienne de moyenne nulle et de variance t.
Autrement dit, la probabilité pour que notre marcheur βt se trouve au temps t dans
un intervalle [a, b] ⊂ R, est donnée par
Z b
1 x2
P(a ≤ βt ≤ b) = √ e− 2t dx
a 2πt
Si l’on traduit ces probabilités de présence, par des répartitions de chaleur sur une barre
rectiligne chauffée en l’origine, il s’en suit que la chaleur est plus élévée autour de la
source, et de plus en plus faible lorsque l’on s’en éloigne.
Par des raisonnements analogues, il est possible vérifier que pour tout couple
d’instants s, t ∈ [0, 1], avec s < t, on a la convergence faible

(βtn − βsn ) −→ (βt − βs )

où (βt − βs ) désigne une variable aléatoire gaussienne de moyenne nulle et de variance
√
t − s. On peut approfondir cette étude, et montrer que le processus gaussien (βt )t∈[0,1]
que nous venons de construire est un processus à trajectoires continues, nulle part
dérivables, à accroissements indépendants et de nature gaussienne.

6.2 La loi des grands nombres

Dans la section précédente, nous avons interprété la répartition de chaleur sur une
barre rectiligne et chauffée en l’origine, en terme de la probabilité de présence d’un
mouvement brownien sur la droite réelle centré en l’origine. Avec cette interprétation,
la répartition de la chaleur sur un segment [a, b] est donnée par la probabilité pour que
βt appartiennent à cet intervalle. On peut écrire cette quantité en terme de moyenne
sur les trajectoires aléatoires conduisant le mouvement brownien à l’instant t dans
l’intervalle [a, b]
P(βt ∈ [a, b]) = E(1[a,b] (βt ))
La loi des grands nombres nous permet d’approcher ces moyennes par des moyennes
empiriques fondées sur la simulation d’un grand nombre N de copies indépendantes βti
de βt . Plus précisément, la probabilité P(βt ∈ [a, b]) est approximativement égale à la
proportion de trajectoires simulées ayant réussi à atteindre cet intervalle au temps t
N
1 X
E(1[a,b] (βt )) ' 1[a,b] (βti )
N
i=1

Dans la théorie de probabilités, ces estimateurs empiriques fondés sur la simulation

de variables aléatoires font partie des algorithmes d’estimation dits de Monte-Carlo.
Leur convergence, lorsque le nombre de simulation augmente, est une conséquence de
la loi des grands nombres. Il existe une variété de théorèmes limites, et d’estimations
de probabilités d’erreurs.
Sans rentrer en profondeur dans l’étude des convergences de ces schémas, on peut se
convaincre de leur qualité numérique par des arguments probabilistes très élémentaires.
Commençons par noter que pour toute fonction fPsuffisament régulière et bornée, les
variances d’erreurs entre l’estimateur empirique N1 N i
i=1 f (βt ), et sa moyenne E(f (βt )),
sont données par la formule suivante
h i2 
1 PN i
E N i=1 f (βt ) − E(f (βt ))

h i2 
1 PN i
=E N i=1 [f (βt ) − E(f (βt ))]

 
E [f (βti ) − E(f (βt ))] [f (βtj ) − E(f (βt ))]
1 PN 1
E [f (βti ) − E(f (βt ))]2 +

P
= N2 i=1 N2 i6=j

1
E [f (βt ) − E(f (βt ))]2 −→ 0


= lorsque N ↑ ∞
N
La dernière assertion provient du fait que

E [f (βti ) − E(f (βt ))]2 = E [f (βt ) − E(f (βt ))]2

et
 
 
E [f (βti ) − E(f (βt ))] [f (βtj ) − E(f (βt ))] = E [f (βti ) − E(f (βt ))] E [f (βtj ) − E(f (βt ))]
= 0
Pour traduire ces estimations de moyennes en terme de probabilités de défauts ou
d’erreurs, on utilise l’inégalité de Markov. Cette inégalité affirme que pour toute variable
aléatoire positive X, et pour tout nombre  > 0, on a

P(X > ) ≤ −1 E(X) (6.1)

Ce résultat résulte simplement des majorations élémentaires suivantes

E(X) = E(X 1X≤ ) + E(X 1X> ) ≥ E(X 1X> ) ≥  E(1X> ) = P(X > )

Notons que l’on a

P(X > ) = P(X 2 > 2 ) ≤ −2 E(X 2 )
Ainsi, si l’on pose
N
1 X
X= f (βti ) − E(f (βt ))
N
i=1
dans l’inégalité précédente, on obtient une estimation de la probabilité de faire une
erreur supérieure à 
N
!
1 X 1
f (βti ) − E(f (βt )) >  ≤ E [f (βt ) − E(f (βt ))]2


P 2
N N 
i=1

Notre situation correspond à la fonction indicatrice f = 1[a,b] . On obtient dans ce cas

E [f (βt ) − E(f (βt ))]2 = E(f (βt )2 ) − E(f (ηt ))2

= E(1[a,b] (βt )) − E(1[a,b] (βt ))

1
= P(βt ∈ [a, b]) (1 − P(βt ∈ [a, b])) ≤
4
En injectant cette majoration dans l’estimation précédente, on en conclut que
N
!
1 X i 1
P 1[a,b] (βt ) − P(βt ∈ [a, b]) >  ≤
N 4N 2
i=1

6.3 Marches aléatoires

D’après les constructions probabilistes précédentes, on peut simuler des mouvements
browniens de deux façons distinctes. La première revient à considérer la marche aléatoire
simple
(βtn0 , βtn1 , . . . , βtnn )
décrite au début de la section 6.1, sur une subdividion suffisamment fine (ti )i=0,...,n de
l’intervalle [0, 1].
 La seconde, est fondé sur les théorèmes limites gaussiens décrits à la fin de la
section 6.1. Supposons que l’on souhaite simuler un mouvement brownien, centré en
l’origine βt0 = 0, sur une suite de pas de temps (ti )i=1,...,n

(βt1 , . . . , βtn )

En terme d’accroissements, cette séquence peut s’exprimer sous la forme suivante :

(βt1 , . . . , βtn ) = [βt1 − βt0 ], βt1 + [βt2 − βt1 ] . . . , βtn−1 + [βtn − βtn−1 ]
2 n
!
X X
= [βt1 − βt0 ], [βti − βti−1 ], . . . , [βti − βti−1 ]
i=1 i=1

On sait de plus que les accroissements

([βt1 − βt0 ], . . . , [βtn − βtn−1 ])

forment une suite de n variables aléatoires gaussiennes de moyenne nulle, et de variances

respectives données par les composantes du vecteur
√ p
( t1 − t0 , . . . , tn − tn−1 )

Ces variables indépendantes peuvent être générées par un algorithme de simulation de

gaussiennes, tel l’agorithme de Box-Muller décrit à la page 46.
Sur un nombre pair de pas de temps, n = 2m, l’algorithme de Box-Muller est fondé
sur la simulation de 2m variables (Uk , Vk )k=1,...,m uniformes, et indépendantes sur [0, 1].
Pour simuler les 2m accroissements gaussiens, ils reste à poser pour indice k = 1, . . . , m
p
(βt2k−1 − βt2k−2 ) = −2(t2k−1 − t2k ) log Uk cos (2πVk )
p
(βt2k − βt2k−1 ) = −2(t2k − t2k−1 ) log Uk sin (2πVk )

6.4 L’équation de la chaleur

Dans les précédentes questions, nous avons montré que la répartition de la chaleur
sur une barre rectiligne était donnée par la distribution gaussienne d’un mouvement
brownien βt sur R. Plus formellement, nous avons pour toute fonction régulière et
bornée Z
E(f (βt )) = f (x) pt (x) dx
R | {z }
P(βt ∈dx)

avec la densité gaussienne

1 x2
pt (x) = √ e− 2t
2πt
Cette fonction spatio-temporelle p : (t, x) ∈]∞[×R 7→ pt (x) ∈]0, ∞[ est la solution
d’une équation aux dérivées partielles, appelée l’équation de la chaleur
∂pt 1 ∂ 2 pt
=
∂t 2 ∂x2
Pour vérifier cette formule, il suffit simplement de dériver la fonction par rapport aux
différentes variables. La dérivée par rapport au paramètre temporel est clairement
donnée par
x2 1 x2 1
   
∂pt 1 1 x2
(x) = √ − + 2 e− 2t = − pt (x)
∂t 2πt 2t 2t 2 t2 t
La dérivée seconde par rapport à la coordonnée spatiale est aussi donnée par la même
formule :
1 ∂ 2 pt 1 ∂ x  1  x2 1 
(x) = − pt (x) = − pt (x)
2 ∂x2 2 ∂x t 2 t2 t

6.5 Une formulation faible

La dérivation de l’équation de la chaleur décrite ci-dessus est simplement basée sur
un jeu de dérivations de fonctions, sans aucun fondement physique. Revenons quelques
secondes à notre interprétation de la répartition de la chaleur en terme de probabilité de
présence d’une particule brownienne. Dans cette interprétation physique, les variations
temporelles des trajectoires moyennes

t ∈]0, ∞[7→ E(f (βt ))

sont données par les variations temporelles de la densité de probabilité gaussienne

Z Z
∂ ∂ ∂pt
E(f (βt )) = f (x) pt (x) dx = f (x) (x) dx
∂t ∂t R R ∂t
Essayons de montrer que pour toute fonction continue et bornée f , deux fois dérivable,
à dérivés continues et bornées, nous avons
∂ 2 pt
Z Z
∂pt
f (x) (x) dx = f (x) (x) dx (6.2)
R ∂t R ∂x2
Cette classe de fonctions est suffisament vaste pour remplir en un certain sens l’ensemble
des fonctions continues et bornées, dual des mesures signées sur R. Cette propriété
topologique, nous permet d’identifier deux mesures signées dès que leurs intégrales sur
ces fonctions test f continues et bornées coincident. La propriété recherchée (6.2) nous
permettra donc de conclure à l’équalité au sens faible des densités
∂pt ∂ 2 pt
(x) = (x)
∂t ∂x2
On parle dans de cas de solution faible de l’équation de la chaleur.
Après ce léger aparté topologique, examinons de plus près les variations temporelles
des moyennes m(t) = E(f (βt )). Par définition de la dérivé temporelle d’une fonction,
et avec quelques abus de notations, nous avons
∂ m(t + ∆t) − m(t) E[f (βt + ∆βt ) − f (βt )]
E(f (βt )) ' =
∂t ∆t ∆t
avec
∆βt = (βt+∆t − βt ) et un pas de temps ∆t ' 0
Le développement de Taylor au second ordre de la fonction f , autour du point βt nous
conduit à la formule
∂f 1 ∂2f
f (βt + ∆βt ) − f (βt ) = (βt ) ∆βt + (βt ) [∆βt ]2 + O((∆βt )3 )
∂x 2 ∂x2
où O((∆βt )3 ) désigne une fonction aléatoire telle que
p √
E(O((∆βt )3 )) ≤ Cte E(( ∆βt )3 ) ≤ Cte ∆t ∆t
pour une certaine constante finie Cte < ∞ ne dépendant que de la fonction f . Après
avoir observé que
 2
∂2f

∂f ∂ f 2
E( (βt ) ∆βt | βt ) = 0 et E (βt ) [∆βt ] | βt = (βt ) × ∆t
∂x ∂x2 ∂x2
on en conclut que
1 ∂2f
 
1
E[f (βt + ∆t ) − f (βt )] = E (βt ) + o(1)
∆t 2 ∂x2
où o(1) désigne une fonction déterministe convergeant vers 0, lorsque ∆t tend vers 0.
En faisant tendre le pas de temps vers 0, nous avons montré l’identité suivante
∂2f
Z Z
∂ ∂pt 1
E(f (βt )) = f (x) (x) dx = (x) pt (x) dx
∂t R ∂t 2 R ∂x2
Il nous reste donc à vérifier que l’on a
∂2f ∂ 2 pt
Z Z
2
(x) p t (x) dx = f (x) (x) dx
R ∂x R ∂x2
On démontre cette ultime propriété à l’aide de deux intégrations par parties successives
∂2f ∂f ∞
Z   Z Z
∂f ∂pt ∂f ∂pt
2
(x) pt (x) dx = pt − (x) (x) dx = − (x) (x) dx
R ∂x ∂x −∞ R ∂x ∂x R ∂x ∂x
et
∞
∂ 2 pt ∂ 2 pt
Z  Z Z
∂f ∂pt ∂pt
(x) (x) dx = f − f (x) (x) dx = − f (x) (x) dx
R ∂x ∂x ∂x −∞ R ∂x2 R ∂x2
Chapitre 7

Dynamiques de population avec

branchements

7.1 Processus de branchements spatio-temporels

Dans la section 3.3.6 nous avons étudié des processus de branchement élémentaires
dans lesquels chaque individu donne naissance à un certain nombre d’enfants. Ces
modèles simplifiés sont assez éloignés de la réalité. Tout d’abord ils ne tiennent pas
compte des explorations des individus dans l’espace. De plus, les taux de branchement
sont totalement indépendant des régions plus ou moins acceuillantes dans lesquelles se
trouvent les populations.
Dans ce qui suit, nous allons essayer de rafiner ces modèles stochastiques pour
rendre compte de ces deux paramètres, et essayer de coller au mieux à la réalité
scientifique. Nous conviendrons que les individus évoluent à chaque instant n dans
un espace d’état En . Les populations d’individus seront représentées par des vecteurs
de Enp . Le paramètre entier p ≥ 0 correspond à la taille des populations. Lorsque p = 0,
on conviendra que l’espace d’état se réduit à un état cerceuil, ou cimetière En0 = {c}.
L’espace d’état du système est donc donné par l’ensemble

Sn = ∪p≥0 Enp

La dynamique d’exploration de chaque individu est associée à des transitions

markoviennes Mn de En−1 vers En . Les mécanismes de branchements dépendent de
fonctions potentiel Gn : EN → [0, ∞) représentant les différents degrés de “fertilité”
du milieu. Plus précisément, un individu sur un site xn ∈ En donnera naissance à un
nombre aléatoire d’enfants gn (xn ) avec

E(gn (xn )) = Gn (xn )

Pour illustrer ce modèle, on peut supposer que les populations évoluent dans le réseau
du plan En = Z2 selon une marche aléatoire simple, et les variables aléatoires de

113
branchement sont de Bernoulli
P (gn (xn ) = [Gn (x)] + 1) = Gn (xn ) − [Gn (xn )]
= 1 − P (gn (xn ) = [Gn (x)])
où [a] désigne la partie entière d’un nombre a ∈ [0, ∞). On peut aussi choisir des
branchements poissonniens
Gn (xn )m
∀m ∈ N P (gn (xn ) = m) = exp (−Gn (xn ) m)
m!
Pour définir plus formellement notre dynamique de population, il convient d’introduire
une suite de variable (gni (xn ))i≥1 indépendantes, et de même loi que gn (xn ). On
supposera de plus que les variables sur des sites distincts sont indépendantes.
Nous sommes enfin en mesure de construire récursivement le modèle. Initialement,
la population X0 est formé d’un seul individu X01 dans l’état x0 ∈ E0
X0 = X01 = x0 ∈ E0p0 = E01
Cet individu donne naissance à
pb0 = g01 (X01 )
enfants que l’on note
Xb0 = (Xb01 , . . . , Xb0pb0 ) ∈ E0pb0
Chacun d’entre eux explore aléatoirement l’espace Xb0i X1i selon la transition M1 .
Cette transition revient à simuler pb0 variables aléatoires X1 de loi M1 (Xb0i , dx1 ). Lorsque
i

cette étape d’exploration est terminée, nous avons une population formée de
p1 = pb0
individus. Cette transition peut s’exprimer de façon synthétique par la formule
suivante :
Xb0 = (Xb01 , . . . , Xbpb0 )
0 X1 = (X11 , . . . , X p1 )
1
Durant l’étape de branchement suivante chaque individu X1i donne naissance à
g11 (X1i )
enfants. À la fin de ce processus nous avons une population formée de pb1
individus
p1
pb1 pb1
X
1
X1 = (X1 , . . . , X1 ) ∈ E1 avec pb1 =
b b b g11 (X1i )
i=1
Chacun de ces individus explore aléatoirement l’espace E2 selon la transition M2 , etc.
branchement exploration
Xn = (Xni )1≤i≤pn −−−−−−→ Xbn = (Xbni )1≤i≤bpn −−−−−−→ Xn+1
Si le système meurt au bout d’un certain temps n, nous avons pbn = 0. Dans ce cas, on
pose
Xbp = Xp+1 = 0
pour tout les instants suivants p ≥ n.
Exercice 7.1.1P 1. Pour toutes fonctions bornées fn sur En , exprimer les variables
pbn bi i
i=1 fn (Xn ) en fonction de pn , Xn , et des variables gn (on utilisera la
aléatoires P
convention ∅ = 0, lorsque la population est éteinte).
2. On considère les mesures empiriques aléatoires
p
X
s(xn ) =def. δxin pour chaque xn = (xin )1≤i≤d ∈ Enp
i=1

Pour toute fonction bornée fn sur En , calculer les moyennes conditionnelles

E(s(Xbn )(fn ) | Xn ) et E(s(Xn+1 )(fn+1 ) | Xbn )

En déduire que

E(s(Xn+1 )(fn+1 ) | Xn ) = s(Xn )(Gn Mn+1 (fn+1 ))

3. Vérifier que les premiers moments des dynamiques de population sont donnés par
la formule suivante
n
Y
E(s(Xn+1 )(fn+1 )) = Ex0 (fn+1 (Xn+1 ) Gk (Xk ))
k=0

7.2 Algorithme génétique

En termes biologiques, les algorithmes génétiques représentent les dynamiques de
population à nombre constant d’individus. Ces individus explorent des régions plus ou
moins acceuillantes, selon des mécanismes de mutation et de sélection. Les individus
meurent ou donnent naissance à des enfants suivant la qualité de leur milieu, mesurée
en terme d’une fonction potentiel.
Dans la section ??, nous montrerons que ces modèles génétiques sont une
simple expression “microscopique” des modèles d’évolution-absorption étudiés dans la
section 4.4.
Leurs interprétations sont à l’image des différents domaines d’applications décrits
à la page 87.

Lorsque la fonction potentiel représente des niveaux de sécurité, tels

des risques de collisions dans des aéroports, l’étape de sélection permet
par exemple de choisir les configurations les moins sécurisées. L’étape de
mutation consiste alors à explorer plus en profondeur les risques de collisions
suivantes. Dans ce contexte, l’agorithme génétique peut s’interpréter comme
des séquences d’évolutions probables conduisant à des collisions d’avion.
 En pharmacologie, la fonction potentiel peut représenter les différents
niveaux de leukocytes dans un organisme vivant. Dans ce contexte, l’étape de
mutation représente les différentes possibilités d’évolution de ces niveaux. La
sélection permet d’évaluer les chutes possibles, et les entrées dans des niveaux
mortels. L’algorithme génétique correspondant permet de décrire les histoires
possibles conduisant aux déces d’un organisme.

Enfin, dans des modèles économiques, les potentiels peuvent représenter

des valeurs de portefeuilles. L’étape de mutation consiste à décrire les
évolutions envisageables de ces quantités partir de niveaux donnés. L’étape de
sélection permet de choisir ces propositions suivant leurs tendances à la hausse
ou à la baisse. Le modèles génétiques correspondent alors à des évolutions de
portefeuilles vers la ruine ou la fortune dans des milieux financiers.

Enfin, dans les modèles d’apprentissage développés en intelligence

artificielle, les potentiels représentent la qualité d’une proposition émise.
L’étape de mutation consiste tout d’abord à élaborer des séquences de
propositions envisageables, plus ou moins bonnes. L’étape de sélection permet
ensuite d’affiner le raisonnement, en choisissant les propositions les mieux
adaptées pour compléter une action donnée, ou pour reconstruire une
information partiellement observée.

En termes mathématiques, un algorithme génétique est une chaı̂ne de Markov Xn =

(Xni )1≤i≤N sur un espace produit E N . Pour fixer les idées on pourra supposer que
l’espace d’état est donné par E = Zd ou E = Rd , avec d ∈ N − {0}. Chacune des
composantes Xni représente la position de l’individu de label i, avec i = 1, . . . , N .
L’évolution de cette chaı̂ne se décompose en deux mécanismes bien distincts. Le
premier correspond à une sélection des individus, selon un certain critère de qualité.
Le second est une exploration pure de l’espace des états. En biologie, cette étape est
souvent appelé mutation, par référence au fait que les codes génétiques des individus
changent au cours du temps.

sélection mutation
Xn = (Xni )1≤i≤N −−−−−−→ X bni )1≤i≤N −−−−−−→ Xn+1 = (Xn+1
bn = (X i
)1≤i≤N
7.2.1 Sélection/Adaptation
Pour décrire l’étape de sélection, on se fixe une fonction potentiel, strictement
positive, G sur l’espace E. Une fois connue la configuration de la chaı̂ne Xn =
(Xni )1≤i≤N au temps n, l’étape de sélection consiste à simuler N v.a. (X b i )1≤i≤N
n
indépendantes de même loi
N
X G(Xni )
PN j
δXni
i=1 j=1 G(Xn )

Autrement dit chaque v.a. X b k choisit l’une des valeurs X i , avec la probabilité
n n
i
G(Xn )
PN j . On remarquera que ce procédé de selection peut aussi s’interpréter comme
j=1 G(Xn )
un mécanisme de naissances et morts. Dans cette interprétation, les individus Xni
disparaissent, ou donnent naissance à un certain nombre de copies.
Il existe divers variantes pour sélectionner les individus les mieux adaptés au
potentiel G. Dans le cas où le potentiel G est à valeurs dans [0, 1], il est bien plus naturel
“d’accepter” chaque individu Xni avec une probabilité G(Xni ), et de le remplacer (avec
une probabilité [1 − G(Xni )]) par un individu choisi avec la loi discrète
N
X G(Xni )
PN j
δXni
i=1 j=1 G(Xn )

Plus formellement, ce mécanisme de sélection est équivalent à poser pour

chaque i = 1, . . . , N
 i
i Xn avec probabilité G(Xni )
Xn =
b
X̃ni avec probabilité 1 − G(Xni )
PN G(Xnj )
où X̃ni désigne une v.a. de loi j=1 P N k
δXnj .
k=1 G(Xn )

Pour des fonctions potentiel pouvant s’annuler sur certaines régions de l’espace,
il est possible que tous les individus aient des potentiels nul. Dans cette situation,
l’algorithme est stoppé.

7.2.2 Mutation/Exploration
Durant la phase de mutation, les individus sélectionnés explorent l’espace
indépendamment les uns des autres, selon des transitions de probabilités élémentaires
M (x, y). Autrement dit, nous avons
bni
X i
Xn+1
i
où Xn+1 bni , .).
désigne une v.a. de loi M (X
 Plus formellement, nous avons
1
P(Xn+1 ∈ dx1 , . . . , Xn+1
N
∈ dxN | X
b1 . . . , X
n
b N ) = M (X
n
b 1 , dx1 ) . . . M (X
n
b N , dxN )
n

À titre d’exemple, si M (x, y) désigne la matrice de transition d’une marche aléatoire

sur Zd , on peut réaliser dynamiquement ces explorations locales en posant
i bi + Ui
Xn+1 =Xn−1 n+1

i
où Un+1 désigne une suite de v.a. indépendantes de même loi p, sur l’ensemble des
vecteurs unitaires directionnels U = {u ∈ Zd : |u| = 1}.
Un exemple schématique d’évolution de N = 4 individus est représenté dans la figure
suivante. Les nombres entiers entre parenthèse correspondent au nombre d’individus sur
le site en question, après l’étape de sélection.

Mutation Selection Mutation Selection Mutation

(0)
(0)

(2) (0)

(1)
(2) (3)

(0)

Fig. 7.1 – Algorithme génétique (N = 4)

Exercice 7.2.1 Décrire mathématiquement, et schématiquement, l’algorithme

génétique sur Z associé à la fonction de potentiel indicatrice G(x) = 1[−L,L] , avec
L ≥ 1. On conviendra que les mutations sont données par les transitions d’une marche
aléatoire sur Z, et l’on initialisera les individus en l’origine.

Exercice 7.2.2 Décrire l’algorithme génétique sur R associé a des mutations

gaussiennes  
1 1 2
M (x, dy) = √ exp − (y − x) dy
2π 2
et un potentiel quadratique centré autour d’un point a ∈ R
 
1 2
G(x) = exp − (x − a)
2

7.3 Modèles d’arbres généalogiques

7.3.1 Modèles non homogènes
L’algorithme génétique décrit dans la section précédente peut être étendu de façon
naturelle à des espaces d’états En dépendants du paramètre temporel n ∈ N. Dans ce
contexte, les individus Xni vivent à chaque instant n dans l’espace En . Les sélections
s’effectuent dans ces mêmes espaces, tandis que les mutations s’expriment comme des
passages aléatoires d’un état de En vers un nouvel état dans En+1 .

Plus formellement, les populations d’individus sont données par des N -uplets

Xn = (Xni )1≤i≤N ∈ EnN et X bni )1≤i≤N ∈ EnN

bn = (X

où EnN désigne l’espace produit (En × . . . × En ), avec N termes.

Dans ce contexte, les transitions de sélection/mutations s’expriment entre des états

non homogènes :

sélection mutation
Xn ∈ EnN −−−−−−→ X bn = E N −−−−−−→ Xn+1 = E N
n n+1

Supposons que les sélections à chaque instant n soient aussi dictées par des potentiels
non homogènes
Gn : xn ∈ En 7→ Gn (xn ) ∈ [0, 1]
Dans ce cas, le mécanisme de sélection dans En , s’exprime sous la forme suivante :
 i i
bni = Xn avec proba Gn (Xn )
X
X̃n avec proba 1 − Gn (Xni )
i

où X̃ni désigne une v.a. de loi

N
X Gn (Xnj )
PN δXnj
k
j=1 k=1 Gn (Xn )
Les individus ainsi sélectionnés Xb i vivent dans l’espace En .
n
Durant la mutation, ces individus passent de l’état En vers un nouvel état En+1 ,
selon des transitions de probabilités Mn+1 (xn , dxn+1 ) de En vers En+1 . Autrement dit,
nous avons
Xb i (∈ En ) i
Xn+1 (∈ En+1 )
n

où Xn+1i désigne une v.a. de loi Mn+1 (Xb i , .). Ces nouveaux individus sont alors
n−1
séléctionnés en fonction d’un potentiel Gn+1 sur En+1 , puis ils mutent de En+1 vers
En+2 selon une transition markovienne de En+1 vers En+2 , etc.

7.3.2 Modèles trajectoriels

Dans ce qui précéde, nous n’avons pas précisé la nature des espaces En , ni a fortiori
celle des transitions Mn (xn−1 , dxn ). Tout ceci semble donc bien abstrait ! Revenons
donc sur terre en supposant que les espaces En sont donnés par des espaces produits
0
. . × E}0
| × .{z
En = E
(n+1)-fois

où E 0 désigne un ensemble quelconque, suffisament réguliers. Pour fixer les idées, on
pourra supposer que E = Zd , ou E = Rd . On conviendra que les points de En sont
donnés par des (n + 1)-uplets représentant des trajectoires de longeur n dans l’espace
E0
xn = (x00 , . . . , x0n ) ∈ En = (E 0 )n+1
L’algorithme génétique précédent est, à chaque étape n ∈ N, formé de N variables
aléatoires trajectorielles, à valeurs dans En , que l’on notera

Xni = X0,n
i i i
∈ En = (E 0 )n+1 ,


, X1,n , . . . , Xn,n 1 ≤ i ≤ N.
et  
bni = X
X i
b0,n ,X i
b1,n ,...,X i
bn,n ∈ En = (E 0 )n+1 , 1 ≤ i ≤ N.

b i représente un chemin dans E 0 de l’origine

Chacun des individus Xni , et Xn
jusqu’au temps n.

Il est important de souligner que la population initiale

X0i = X00 i ∈ E0 = E 0 avec 1≤i≤N

est tout simplement formée de N variables aléatoires à valeurs dans E 0 .

Sélection trajectorielle :
Dans notre cadre trajectoriel, le mécanisme de sélection, peut s’interpréter comme une
sélection de trajectoires, en fonction des différents potentiels

Gn (Xni ) = Gn X0,n
i i i


, X1,n , . . . , Xn,n

L’expression de cette transition reste inchangée. Nous avons à nouveau

 i
i Xn avec proba Gn (Xni )
Xn =
b
X̃ni avec proba 1 − Gn (Xni )
 
où X̃ni = Xbi , X
0,n
bi , . . . , X
1,n
bi
n,n désigne une v.a. de loi

N N j j
X Gn (Xnj ) X Gn (X0,n , . . . , Xn,n )
δXnj = δ(X j j
0,n ,...,Xn,n )
PN k
PN k k
j=1 k=1 Gn (Xn ) j=1 k=1 Gn (X0,n , , . . . , Xn,n )

b i vivent désormais dans l’espace de chemins

Les trajectoires ainsi sélectionnés Xn
En = (E 0 )n+1 .

Mutation trajectorielle :

L’étape de mutation dépend uniquement de la nature des transitions Mn+1 (xn , dxn+1 )
de En dans En+1 . Supposons que ces dernières correspondent aux transitions de
probabilités d’un processus historique associé à une évolution markovienne sur E 0 .
Dans cette situation, on rappelle que ces transitions Mn+1 (xn , dyn+1 ) s’expriment
sous la forme suivante :

Mn+1 ((x00 , . . . , x0n ), d(y00 , . . . , yn+1

0
)) = δ(x00 ,...,x0n ) (d(y00 , . . . , yn0 ))Mn+1
0
(yn0 , dyn+1
0
)
0
Mn+1 est une transition de Markov de E 0 vers lui même. Autrement dit, simuler une
variable aléatoire trajectorielle de loi

Mn+1 ((x00 , . . . , x0n ), d(y00 , . . . , yn+1

0
))

revient à conserver tout d’abord le segment de trajectoire

xn = (x00 , . . . , x0n )
0
On lui adjoint ensuite une extension élémentaire aléatoire de loi Mn+1 (x0n , dyn+1
0 )

x0n 0
yn+1
 Durant l’étape de mutation correspondante
 
bni = X 0 i 0 i
X b0,n ,...,X
bn,n ∈ EnN
| {z }
↓
z }| { 
i 0 i 0 i 0 i N
Xn+1 = (X0,n+1 , . . . , Xn,n+1 ), Xn+1,n+1 ∈ En+1

0
chaque chemin s’étend selon un déplacement élémentaire de loi Mn+1 , c’est à dire

 
i
Xn+1 = [X0,n+1 , . . . , Xn,n+1 ], Xn+1,n+1 ∈ En+1 = (E 0 )n+2
0 i 0 i 0 i
| {z }
||
z }| { !
 
0 i 0 i 0 i bi , X0 i
= [Xb , ... , X
0,n
b
n,n ], X n+1,n+1 = Xn n+1,n+1

0 i
avec une variable aléatoire Xn+1,n+1 0
de loi Mn+1 b 0 i , .).
(Xn,n

0
À titre d’exemple, dans le cas où Mn+1 désigne la matrice de transition d’une marche
0 d
aléatoire sur E = Z , on peut réaliser dynamiquement ces explorations locales en
posant
0 i 0 i i
Xn+1,n+1 =X bn,n + Un+1
i
où Un+1 désigne une suite de v.a. indépendantes de même loi p, sur l’ensemble des
vecteurs unitaires directionnels

U = {u ∈ Zd : |u| = 1}

Lignes ancestrales :

L’algorithme génétique trajectoriel décrit ci-dessus correspond bien à une évolution

d’arbres généalogiques.

Dans cette interprétation, chaque trajectoire

Xni = X0,n i i i
∈ En = (E 0 )n+1


, X1,n , . . . , Xn,n
i
représente la ligne ancestrale de l’individu courant Xn,n à la date n. Les
i
coordonnées Xp,n , avec 0 ≤ p ≤ n, correspondent aux différents ancêtres
de cet individu, à chaque niveau temporel 0 ≤ p ≤ n.
Dans ce contexte, X0,n i i
représente l’ancêtre initial, X1,n i
sa première descendance, X2,n
sa seconde, etc. Dans le modèle général que nous avons développé, les disparitions et les
sélections de lignes ancestrales dépendent de la qualité de la trajectoire de descendance
complète. Les différents degrés d’adaptation d’une lignée complète

xn = (x00 , . . . , x0n ) ∈ En

sont mesurés par une fonction potentiel Gn (x00 , . . . , x0n ). Lorsque ces fonctions ne
dépendent que des composantes terminales, c’est à dire lorsque l’on a

Gn (x00 , . . . , x0n ) = G0n (x0n ) ,

la sélection des lignes ancestrale ne dépend que de la qualité d’adaptation du dernier

descendant. Dans cette situation, nous laissons le soin au lecteur de se convaincre que
ce modèle trajectoriel correspond à l’arbre généalogique d’une population d’individus
explorant l’espace E 0 selon Mn0 , et s’adaptant en fonction des potentiels de selection
G0n .

7.4 Chaı̂nes renforcées

Ces modèles de renforcement sont associés à des évolutions aléatoires où chaque
transition dépend de la mesure d’occupation des sites visités dans le passé. Ces processus
permettent de représenter des stratégies humaines d’exploration de sites comme des
rues, des magasins, des restaurants, ou tout autre endroits d’une ville. Dans ce contexte,
un site donné est d’autant plus attractif s’il répond à un certain critère de qualité, ou
lorsqu’il a été déjà visité de nombreuses fois par le passé !
Le critère de qualité est représenté par une fonction potentiel G : E → [0, 1], sur un
espace d’états E. Les valeurs de G en un site x sont d’autant plus grandes que le site
est attrayant. Ainsi, notre explorateur se trouvant en Xn = x, au temps n, choisit d’y
rester avec une probabilité G(x), soit préfère retourner vers l’un des sites précédemment
visité X0 , . . . , Xn−1 . Dans cette situation, il choisit un nouveau site avec une probabilité
n−1
X G(Xp )
Pn−1 δXp
p=0 q=0 G(Xq )
 Cette sélection aléatoire peut être vue comme une transition élémentaire

Xn X
bn

de probabilités de transitions
n−1
X G(Xp )
G(Xn ) δXn + (1 − G(Xn )) Pn−1 δXp
p=0 q=0 G(Xq )

Partant du site sélectionné X

bn , l’individu effectue une nouvelle exploration de
la région Xn
b Xn+1 , selon une transition de probabilité Mn+1 . Autrement
dit, nous avons

P(Xn+1 ∈ dy | X
bn = x) = Mn+1 (x, dy)

Cette transition Mn peut, par exemple, représenter une exploration uniforme des
sites voisins à X
bn .
Le lecteur aura certainement noté que Xn n’est pas une chaı̂ne de Markov, mais le
processus historique possède toujours la propriété markovienne.

Exercice 7.4.1 Examiner la situation où le potentiel est constant G(x) = , avec
 ∈ [0, 1]. Établir une analogie avec la loi des grands nombres.