Modèles et Simulation stochastiques

N . Bessah

13 mars 2022
Table des matières

1 Simulation 1
1.1 Génération de v.a. uniformes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Générateurs de nombres aléatoires . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 Qualité des nombres aléatoires . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Génération de v.a. non uniformes . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Méthode de l’inverse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 Méthodes particulières pour les variables aléatoires non uniformes . . . . 4
1.3 Méthodes de Monté Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.1 Calcul d’une intégrale définie par la méthode de Monté Carlo . . . . . . 7
1.4 Réduction de la variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

i
Chapitre 1

Simulation

Introduction La simulation est un ensemble de techniques permettant de reproduire approxi-

mativement le comportement d’un système. On peut simuler un système industriel, économique,
social, politique ou biologique.

Exemples Un système d’attente : Son état peut être décrit par le nombre de clients dans
la file, le temps d’attente, etc. Les événements sont les arrivées et les départs de clients, qui
peuvent se produire aléatoirement.
On a aussi recours à la simulation pour le calcul d’intégrales ( difficile à faire analytiquement),
en optimisation, etc...

1.1 Génération de v.a. uniformes

Un nombre au hasard x dans l’intervalle [0, 1] (resp. dans l’ensemble {0, ..., 9}) est une
réalisation de la v.a. X, de loi U [0, 1] (respectivement U {0, . . . , 9}).
On peut produire un nombre au hasard par des moyens mécaniques ou par ordinateur

1.1.1 Générateurs de nombres aléatoires

La génération de nombres aléatoires a connu plusieurs étapes :
1. Génération suivant un processus physique (lancer d’un dé, une pièce..) mais : La séquence
n’est pas reproductible donc impossible à tester ;
2. Etablissement de tables de nombres aléatoires : Les chiffres contenus sont équiprobables,
mais la taille de l’échantillon peut dépasser la capacité de la table, donc risque d’utiliser
le même échantillon ;
3. Mise au point de méthodes algorithmiques : Elles sont faciles à programmer, les séquences
sont reproductibles (i.e. en relançant le générateur avec les mêmes conditions initiales on
obtient la même séquence de nombres.
Les algorithmes les plus employés sont basés sur des suites récurrentes d’entiers.

Générateurs congruentiels linéaires

Méthodes congruentielles multiplicatives

xi+1 = axi mod [m]

1
CHAPITRE 1. SIMULATION 2

(reste de la division de axi par m) m et a sont des paramètres à choisir convenablement de sorte
que la longueur du cycle soit égale à m − 1. La racine du générateur x0 appelée seed ou graine
doit être fixée.
Pour m = 2M , la période maximale est 2M −2 , obtenue si a ≡ 3 mod [8] ou a ≡ 5 mod [8]
et le seed est impair.

Exemple 1 :
xn+1 = 5xn mod 25

pour

x0 = 1

on obtient la suite

5, 25, 29, 17, . . . 28, 9, 13, 1, 5, . . . (période maximale)

et on a pour

x0 = 2
10, 18, 26, 2

Méthodes congruentielles mixtes

xi+1 = (axi + b) mod [m]
On prend m le plus grand possible pour avoir une grande période. Si m = 2M , et b > 0
la période maximale est 2M obtenue si a ≡ 1 mod [4] et b impair. La séquence peut être
initialisée par n’importe quel nombre  entre 1 et m − 1
compris
xi xi

Les nombres ∈ [0, 1] ∈ [0, 1[ sont des nombres aléatoires. (Résultat dû a Hull et
m−1 m
Dobell (1962)).

Exemple
Xi = (9Xi−1 + 3) mod 25 i = 1, 2,
X0 = 1
on obtient les valeurs suivantes : 12, 5, 16, . . .

Remarque Les nombres obtenus devraient être uniformément distribués sur [0, 1]. Une
représentation graphique des pairs (x1 , x2 ), (x2 , x3 ), . . . (xi , xi+1 ), . . . devrait montrer une ré-
partition uniforme de ces points dans le carré. Ceci n’est pas toujours vrai, car la relation :
xi+1 = (axi + b)mod[m] signifie que xi+1 = (axi + (b − km) où k ∈ {0, 1, ..., a − 1} qui est
l’équation d’une droite. D’où la nécessité de choisir m grand ( on recommande un m > 230 ).

Définition (Lecuyer 1994 ) D’une manière générale, on définit un générateur de nombres

aléatoires par la donnée d’un ensembles d’états S (d’ entiers, vecteurs,...), d’une fonction de
transition f : S 7−→ S, tel que si+1 = f (si ), ( s0 valeur initiale ou Seed ) et une fonction
de sortie g : S 7−→ [0, 1] (ensemble de sortie ) telle que ui = g (si ). Ces valeurs sont les
nombres aléatoires produits par le générateur. La période du générateur est min (p) vérifiant :
∃ n0 , ∀n ≥ n0 , sn+p = sn (Voir site web http ://www.iro.umontreal. ca/ lecuyer).
CHAPITRE 1. SIMULATION 3

1.1.2 Qualité des nombres aléatoires

Un générateur de nombres aléatoires est dit valide s’il vérifie les propriétés suivantes :
1. les nombres générés sont uniformes dans [0, 1] ;
2. non corrélés ;
3. longueur du cycle élevée.
Remarques
— Le choix de m sous la forme 2M a pour but de rendre facile et rapide le calcul par
ordinateur. M est le nombre de bits utilisés pour représenter un entier positif dans la
machine.
— Le plus grand entier est m = 216 resp.( 232 , 264 ) pour les entiers à 16 bits resp(32ou64
bits) la suite obtenue varie dans {0, . . . , m − 1}.
— Il existe des instructions telles que Randomize, srand (en C) permettant de régler la
valeur initiale à l’horloge, ce qui évite de re-générer le même échantillon.

1.2 Génération de v.a. non uniformes

1.2.1 Méthode de l’inverse
v.a. continue
Soit X une v.a. continue de fonction de répartition F ; la variable aléatoire U = F (X) est
de loi uniforme sur [0, 1].
En effet P (U < x) = P (F (X) < x) = P (X < F −1 (x)) = F (F −1 (x)) = x (car F est
monotone et continue).
Inversement si U une v.a. uniforme sur [0, 1], alors la variable aléatoire F −1 (U ) a la même
loi que X
En effet P (F −1 (U ) < x) = P (U < F (x)) = F (x)
On peut donc simuler une variable aléatoire continue X de fonction de répartition F sachant
qu’on sait simuler U ; la difficulté réside parfois dans le calcul d’une fonction de répartition et
de son inverse.

Remarque dans le cas où F n’est pas strictement croissante donc non inversible, on utilise
l’inverse généralisé de la fonction de répartition défini par

F −1 (s) = inf{x, tel que F (x) ≥ s} ou bien F −1 (s) = sup{x, tel que F (x) ≤ s} (1.1)

Exemple Soit X une variable aléatoire suivant la loi E(λ)

La fonction de répartition est F (x) = 1 − exp(−λx) pour x ≥ 0
Soit U une variable aléatoire uniforme sur [0,1], U = F (X) ⇔ U = 1 − exp (−λX) ⇔
1
− ln (1 − U) = X.
λ
D’où pour simuler une réalisation de la variable aléatoire X, il suffit de simuler une réalisation
uniforme.

Remarque Puisque la variable aléatoire 1−U est également uniforme sur [0, 1] on la remplacera
souvent par la variable V uniforme sur [0, 1] et ce afin de réduire le nombre d’opérations.
CHAPITRE 1. SIMULATION 4

v.a. discrète
Soit X une variable aléatoire discrète dont les valeurs sont x0 , x1 , x2 , . . . , xn et de probabilités
respectives p0 , p1 , ..., pn . L’algorithme pour générer une réalisations de X est décrit comme suit :

Algorithm 1 Générer une réalisation d’une v.a. discrète

1: Générer U ∼ U [0, 1]
2: if 0 ≤ U < p0 then X = x0
3: else if p0 ≤ U < p0 + p1 then X = x1 .
4: etc..
5: end if

D’une manière générale, soit U une variable aléatoire uniforme sur [0, 1].
i−1
X i
X
Si pj ≤ U < pj alors X = xi (1.2)
j=0 j=0

Justification
 
i−1
X i
X i
X i−1
X
P pj ≤ U ≤ pj  = pj − pj
j=0 j=0 j=0 j=0

= pi = P (X = xi )

1.2.2 Méthodes particulières pour les variables aléatoires non uni-

formes
Génération d’une variable aléatoire de loi N(0,1)
Rappel : Théorème central limite Soient X1 , X2 , ..., Xn des variables aléatoires indépen-
dantes de même loi de moyenne µ et d’écart type σ(f ini) alors la variable aléatoire
n
Xi − nµ
P
i=1
√ (1.3)
nσ
tend asymptotiquement vers une loi normale N (0, 1)

Méthode par approximation Pour générer une v.a. selon une loi N (0, 1) on utilisera le
TCL pour l’approcher par des v.a. de loi uniforme.
On simule pour cela n v.a.i.i.d. suivant une loi uniforme sur [0, 1],
n n
Ui −
P
i=1 2
X= q
n/12

est asymptotiquement normale. L’expérience a montré que la convergence est rapide (à partir
1 1
de n = 3) comme E (U ) = et var (U ) = et afin de simplifier les calculs on prendra n = 12
2 12
et l’on a
12
X
Ui − 6 suit une loi N (0, 1)
i=1
CHAPITRE 1. SIMULATION 5

Méthode exacte, ( BOX et MULLER (1958) ) Si U1 , U2 sont deux variables aléatoires

indépendantes de loi uniforme sur [0, 1] alors
q q
N1 = −2 log U1 cos(2πU2 ) et N2 = −2 log U1 sin(2πU2 )
sont deux variables aléatoires indépendantes de loi N (0, 1).

Preuve Soit N1 et N2 deux v.a.i de même loi N (0, 1). Considérons N1 et N2 comme coordonnées
cartésiennes d’un point. Par transformation en coordonnées polaires, on obtient
N1 = R cos θ

et

N2 = R sin θ
où R et θ sont deux variables aléatoires.
Calculons la loi de R et θ.
1
N12 + N22 = R2 =⇒ R2 ∼ χ21 = E( )
2
d’où
r2
−
P (R < r) = P (R2 < r2 ) = 1 − e 2

et la densité de probabilité est

r2
−
fR (r) = re 2 r>0

Cherchons la loi de θ et montrons que les v.a. R et θ sont indépendantes.

La densité conjointe de (N1 , N2 ) est le produit des deux densités de la loi normale centrée
réduite.

f (n1, n2 ) = f1 (n1 )f2 (n2 )

n21 n22
1 − 1 −
=√ e 2 √ e 2
2π 2π
(n1 + n22 )
2
1 −
= e 2
2π
La densité du couple (R, θ) est

g(r, θ) = f (n1, n2 ) |J|

où J est le jacobien

∂n1 ∂n1
J= ∂r ∂θ =r
∂n2 ∂n2
∂r ∂θ
CHAPITRE 1. SIMULATION 6

d’où

((r cos θ)2 + (r sin θ)2 )

1 −
g(r, θ) = e 2 r
2π
r2
1 −
= e 2r
2π
r2
− 1
re 2 est la densité de R ce qui implique que est la densité de θ (loi uniforme sur [0, 2π]
2π
de plus les deux v.a. sont indépendantes ( puisque leur densité conjointe est le produit des deux
densités)

Méthode de rejet simplifiée

Cas continu Soit X une v.a. continue de densité de probabilité f définie sur un intervalle
compact [a, b]
— f est bornée i.e. il existe une constante M telle que f (x) ≤ M pour tout x ∈ [a, b] ;
— on insère cette fonction dans un rectangle qui a pour cotés (0, M ) et (a, b) et on répartit
des points uniformément à l’intérieur du rectangle ;
• Soit P (x, y) un de ces points, ses coordonnées sont : x = (b − a)u + a et y = M v
où u et v sont des réalisations de variables aléatoires uniformes sur (0, 1).
— Les points sous la courbe sont des réalisations de la va X. En effet, soit

C = {(x, y) : α ≤ x ≤ β, y ≤ f (x)}
Rβ
f (t)dt
P robabilité(p ∈ C) = Rαb
a f (t)dt
Rβ
αf (t)dt
=
1
= P (α ≤ X ≤ β)

x est bien une réalisation de X.

Algorithm 2 Méthode de rejet simplifiée pour une réalisation de X

1: Générer x ∼ U [a, b]
2: Générer y ∼ U [0, M ]
3: if y ≤ f (x) then sortir x
4: else aller à 1
5: end if

Si le point P est sous la courbe de f , x est une réalisation de X sinon le point est rejeté.

La probabilité d’acceptation (ou de rejet) Un point est accepté s’il est sous la courbe
Rb
af (t)dt
P (acceptation) = Pa =
M (b − a)
1
= .
M (b − a)
La probabilité de rejet est 1 − Pa , quand M décroit Pa croit
CHAPITRE 1. SIMULATION 7

Inconvénient de la méthode
— La méthode est limitée puisque les densités doivent être définies uniquement sur un
intervalle borné (les variables uniformes ne sont simulées que sur ce type d’intervalle.
— La probabilité de rejet risque d’être importante quand l’allure de la densité présente des
pics.
La méthode suivante permet d’éviter ce type d’inconvénients.

Méthode de rejet généralisée

Supposons qu’il existe une fonction de densité g (dont nous savons simuler la variable
aléatoire associée) et il existe une constante a telle que quel que soit x, f (x) ≤ ag(x)

Algorithm 3 Méthode de rejet généralisée pour une réalisation de X

1: Générer x suivant la densité de probabilité g
2: Générer y ∼ U [0, ag(x)]
3: if y ≤ f (x) then sortir x
4: else aller à 1
5: end if

Probabilité de rejet

(ag(x) − f (x))dx
R
Prej = R
ag(x)dx
1
=1−
a
quand a augmente la probabilité de rejet croit d’où l’intérêt de choisir a petit.

1.3 Méthodes de Monté Carlo

Les méthodes de Monté Carlo sont l’application de la simulation à l’estimation.

1.3.1 Calcul d’une intégrale définie par la méthode de Monté Carlo

Méthode 1, succès échec
Soit f une fonction définie sur (a, b) telle que sup {f (x), x ∈ (a, b)} = c < ∞
Rb
En supposant que f > 0, il est question d’estimer la valeur de l’intégrale I = f (x)dx
a
1. On génère n valeurs (x1 , . . . , xn ) de loi uniforme sur [a, b]
2. en considérant les valeurs xi comme abscisses des points, on génère donc leurs ordonnées
yi ∼ U [0, c]
si yi < f (xi ) alors le point est sous la courbe.
Soit nf le nombre de points appartenant à la courbe alors

Rb
nf f (x)dx Zb
nf
a
= ⇒ f (x)dx = c(b − a) (1.4)
n c(b − a) a
n
quand n tend vers l’infini l’estimation s’améliore.
CHAPITRE 1. SIMULATION 8

Méthode 2 (moyenne)
Elle consiste à exprimer une quantité à estimer sous la forme
R
d’une espérance.
Soit f une fonction définie sur un domaine D et I = f (x)dx la valeur à estimer. On
D
cherche à reconnaitre une densité de probabilité définie sur D i.e. voir si f (x) = g(x).h(x) où
g(x) est une densité de probabilité définie sur D. Dans ce cas,
Z Z
I= f (x)dx = h(x)g(x)dx = E(h(X)) (1.5)
D D

espérance mathématique de la variable aléatoire h(X) de densité g dont l’ estimateur est la

moyenne empirique.
N
1 X
h(Xi ) où X1 , . . . , XN est N échantillon de X (1.6)
N i=1

Dans le cas particulier où le domaine D est de la forme (a, b), a et b finis il suffit d’introduire la
1
constante qui est la densité d’une loi uniforme sur (a, b)
b−a
Zb
I= f (x)dx
a
Zb
1
= (b − a) f (x) dx
a
b−a
= (b − a)E(f (U )) U ∼ U [a, b]

l’ estimateur de I est
N
(b − a) X
f (Ui )
N i=1

R1 √ R∞
Exemples I1 = 1 − x2 dx I2 = exp(−x)(1 − exp(−x))2 dx
0 0

Théorème : Loi forte des grands nombres soit (Xn )n une suite de v.a.i.i.d. On suppose
E (|X|) < ∞ alors
n
1X
E(X) = n→∞
lim Xi P.s (1.7)
n i=1
La loi forte des grands nombres nous garanti la convergence des méthodes de Monté Carlo ; la
précision de l’approximation est mesurée en utilisant le théorème central limite. Soit
n
1X
n = E(X) − Xi
n i=1

d’après le TCL
√
nn
→ N (0, 1)
σ
CHAPITRE 1. SIMULATION 9

d’où
" √ # c2
nn Z
P c1 < < c2 = f (x)dx
σ c 1

où f est la densité de la loi normale i.e.

x2
1 −
f (x) = √ e 2 x ∈ R
2π

Erreur de Monté Carlo Nous pouvons présenter

n
1X σ
E(X) = Xi ± √ (1.8)
n i=1 n

ou bien on donne un intervalle de confiance

n n
" #
1X σ 1X σ
Ic = xi − tα √ , xi + tα √ , tα
n i=1 n n i=1 n

est tel que

" √ #
nn
P < tα = 1 − α
σ

exemple

α = 5% tα = 1.96.

Remarque Quand n augmente la précision s’améliore, or dans Ic la longueur de l’intervalle

σ
est h = 2tα √ si l’on souhaite réduire h de moitié cela nécessite une taille d’échantillon N
n
σ σ
vérifiant 2tα √ = tα √ i.e. N = 4n.
N n

1.4 Réduction de la variance

L’estimation d’un caractère nécessite de donner une précision, cela est possible en augmentant
la taille de l’échantillon. Cependant, cette augmentation est parfois couteuse d’où l’idée de
réduire la variance de l’estimation en corrigeant l’échantillon de valeurs tiré de manière à être
plus proche d’une distribution uniforme.

Théorème soient ui des v.a.i.i.d. suivant la loi uniforme sur (a, b) et ϕ une fonction monotone
de chacune de ses variables.
1 1
X = (ϕ(u1 , . . . , un ) + ϕ(un+1 , . . . , u2n )) et Y = (ϕ(u1 , . . . , un ) + ϕ(1 − u1 , . . . , 1 − un ))
2 2
vérifient var(Y ) ≤ var(X)
Bibliographie

[1] A.H.S. Ang and W.H. Tang. Probability Concepts in Engineering Planning and Design :
Decision, risk and reliability. Probability Concepts in Engineering Planning and Design.
Wiley, 1975.
[2] T Bodineau. Modélisation de phénomènes aléatoires : introduction aux chaînes de markov
et aux martingales. Ecole Polytechnique, 2015.
[3] Pierre Brémaud. Initiation aux Probabilités : et aux chaînes de Markov. Springer Science
& Business Media, 2009.
[4] JT Byron. Morgan, elements of simulation, 1984.
[5] J S Dagpunar. Simulation and Monte Carlo With applications in finance and MCMC. 2007.
[6] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from
incomplete data via the em algorithm. Journal of the Royal Statistical Society : Series B
(Methodological), 39(1) :1–22, 1977.
[7] Luc Devroye. Non-uniform random variate generation (1986), 2008.
[8] Yadolah Dodge and Giuseppe Melfi. Premiers pas en simulation. Springer edition, 2008.
[9] D.Petris. simulation. 1999.
[10] James E Gentle. Random number generation and Monte Carlo methods. Springer Science
& Business Media, 2006.
[11] WR Gilks, S Richardson, and DJ Spiegelhalter. Markov chain monte carlo in practice
chapman & hall : London. Markov Chain Monte Carlo in practice. Chapman and Hall,
London., 1996.
[12] Lena Gorelick, Meirav Galun, Eitan Sharon, Ronen Basri, and Achi Brandt. Shape
representation and classification using the poisson equation. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 28(12) :1991–2005, 2006.
[13] Marco Gori and Augusto Pucci. Research paper recommender systems : A random-walk
based approach. In 2006 IEEE/WIC/ACM International Conference on Web Intelligence
(WI 2006 Main Conference Proceedings)(WI’06), pages 778–781. IEEE, 2006.
[14] Geoffrey Grimmett, David Stirzaker, et al. Probability and random processes. Oxford
university press, 2001.
[15] Jean-François Hêche, Thomas M Liebling, and Dominique De Werra. Recherche opération-
nelle pour ingénieurs, volume 2. PPUR presses polytechniques, 2003.
[16] Dirk P Kroese, Thomas Taimre, and Zdravko I Botev. Handbook of monte carlo methods,
volume 706. John Wiley & Sons, 2013.
[17] Linyuan Lü and Tao Zhou. Link prediction in complex networks : A survey. Physica A :
statistical mechanics and its applications, 390(6) :1150–1170, 2011.
[18] Marina Meilă and Jianbo Shi. A random walks view of spectral segmentation. In Interna-
tional Workshop on Artificial Intelligence and Statistics, pages 203–208. PMLR, 2001.

10
BIBLIOGRAPHIE 11

[19] Adam Prügel-Bennett. The Probability Companion for Engineering and Computer Science.
Cambridge University Press, 2020.
[20] Lawrence R Rabiner. A tutorial on hidden markov models and selected applications in
speech recognition. Proceedings of the IEEE, 77(2) :257–286, 1989.
[21] Christian Robert and George Casella. Monte Carlo statistical methods. Springer Science &
Business Media, 2013.
[22] Sheldon M. Ross. Simulation. USA, fourth edition, 2006.
[23] Reuven Y Rubinstein and Dirk P Kroese. Simulation and the Monte Carlo method,
volume 10. John Wiley & Sons, 2016.
[24] Alan Ruegg. Processus stochastiques : avec applications aux phénomènes d’attente et de
fiabilité, volume 6. PPUR presses polytechniques, 1989.
[25] Purnamrita Sarkar and Andrew W Moore. Random walks in social networks and their
applications : a survey. In Social Network Data Analytics, pages 43–77. Springer, 2011.
[26] Feng Xia, Jiaying Liu, Hansong Nie, Yonghao Fu, Liangtian Wan, and Xiangjie Kong.
Random walks : A review of algorithms and applications. IEEE Transactions on Emerging
Topics in Computational Intelligence, 4(2) :95–107, 2019.