Simulation Agreg 1718

Modélisation option A Simulation
Simulation de variables aléatoires

4 avril 2018
Pierre Rousselin
La dernière version de ce document est disponible sur le site :
https://www.math.univ-paris13.fr/~rousselin/
N’hésitez pas à me faire part de vos commentaires et à me signaler les coquilles, erreurs et
imprécisions à l’adresse indiquée sur ce site.
1 Introduction
Bien que le sujet soit intéressant, nous ne nous attarderons pas trop sur les les questions
que posent la simulation du hasard par un ordinateur. Pour en savoir plus sur ce sujet, voir par
exemple le début du chapitre 1 de [BC07].
On suppose donc disposer d’un « bon » programme renvoyant lors de ses appels successifs des
nombres réels aléatoires indépendants suivant la loi uniforme sur [0 ; 1]. Les problèmes surgissent
dès la définition : d’une part les nombres réels ne sont pas représentables en machine. On travaille
le plus souvent avec des nombres à virgule flottante en double précision qui sont en nombre fini,
bien qu’important. D’autre part, l’ordinateur est une machine déterministe.
On dispose donc d’algorithmes permettant de générer des suites de nombres, de façon déter-
ministe, qui « ressemblent » aux premières décimales d’une suite de réalisations indépendantes
de la loi uniforme sur [0 ; 1]. Le choix de tel ou tel algorithme dépend de l’utilisation que l’on en
fait.
Pour la simulation, l’algorithme le plus utilisé aujourd’hui est sans doute « Mersenne Twister
19937 », mais son utilisation est fortement déconseillée lorsqu’on souhaite, par exemple, fabriquer
une clé privée aléatoirement car il est jugé trop « prévisible ». En revanche, il passe un grand
nombre de tests statistiques et est assez performant.
Le logiciel Scilab est assez généreux en ce qui concerne la simulation. Les deux principales
fonctions que nous utiliserons seront rand, le générateur basique de Scilab, qui utilise un gé-
nérateur congruentiel simple et permet de simuler la loi normale centrée réduite et grand une
gigantesque fonction permettant de choisir son générateur parmi une liste assez longue (avec
Mersenne Twister 19937 comme choix par défaut) et de simuler un nombre important de lois de
probabilités.
Pour la suite de ce chapitre, nous utiliserons uniquement la fonction prédéfinie rand et nous
travaillerons sous l’hypothèse que ses appels successifs fournissent des réalisations indépendantes
de variables uniformes sur [0 ; 1]. En général, on peut imposer un état initial au générateur de
nombres aléatoires en choisissant une graine (en anglais, seed). Voyons le comportement par
défaut de Scilab.
Exercice 1.1. rand, rand("seed", s), rand(m, n), histplot

1. Lancer Scilab et taper rand puis « Entrée » trois fois. Noter les résultats.
2. Quitter Scilab et recommencer l’expérience précédente.
3. Quitter Scilab (promis c’est la dernière fois), et entrer la commande rand("seed") dans la
console. Taper rand puis redemander la graine comme précédemment.
4. Entrer la commande rand("seed", 0), puis refaire un appel à rand
5. Faire le bilan. Quelle est l’influence de la graine ? Comment la modifier ? Comment retrouver
l’état initial par défaut du générateur rand ?
6. Que fait la commande rand(3, 5) ?
Université Paris 13 1/42 Préparation à l’agrégation 2017 - 2018

7. Que fait la commande histplot(linspace(0, 1, 10), rand(1, 10000)) ? Consulter l’aide

de histplot
8. Entrer plusieurs fois de suite la commande getdate(), et aller voir sa documentation.
9. Que fait le script suivant ?
x = getdate ( " s " ) ;

rand ( " seed " , x ) ;
2 Lois uniformes
Nous commençons par simuler la loi uniforme sur un intervalle réel ou entier, en donnant
une définition qui englobe ces deux cas et bien d’autres encore. Nous en profitons au passage
pour donner quelques rappels sur le calcul des moments de ces lois.
Définition 1. Soit (E, E, µ) un espace mesuré, et A ∈ E tel que 0 < µ(A) < ∞ . La loi uniforme
sur A (relativement à la mesure µ) est la loi de probabilité νA définie par :
µ (A ∩ B)
∀B ∈ E, νA (B) = .
µ (A)
Exemple 2. 1. Prenons E = R, E = B (R) et pour µ, la mesure de Lebesgue. La loi uniforme sur

1[a;b]
l’intervalle [a ; b] (ou ]a ; b], ou [a ; b[, ou ]a ; b[, cela ne change rien) est la loi de densité b−a
par rapport à la mesure de Lebesgue.
Qd
2. On peut, de même, caractériser la loi uniforme sur un pavé de Rd de la forme i=1 [ai ; bi ] par
sa densité par rapport à la mesure de Lebesgue :
d
Y 1[ai ;bi ] (xi )
(x1 , x2 , . . . , xd ) 7→ .
i=1
bi − ai
3. Prenons E = Z, muni de la tribu discrète 2Z et pour µ, la mesure de comptage (c’est-à-dire,

pour tout A ⊂ Z, µ(A) = Card(A)) . On peut définir ainsi la mesure uniforme sur toute partie
finie non vide de Z. Ceci est encore valable en remplaçant Z par tout ensemble dénombrable.
Exercice 2.1.
Soient a et b deux réels avec a < b. Soit U une variable aléatoire de loi uniforme sur [a ; b].
a+b (b−a)2
Montrer que l’espérance de U vaut 2 et que sa variance vaut 12 .
Exercice 2.2.
En utilisant la fonction rand, écrire une fonction runif(m,n,a,b) qui renvoie une matrice de
taille m × n dont les coefficients sont des réalisations de la loi uniforme sur [a ; b].
Exercice 2.3.
Soient a et b deux entiers, avec a < b. Notons Ja ; bK l’ensemble {a, a + 1, . . . , b} et n = b − a + 1.
Soit X une variable aléatoire de loi uniforme sur Ja ; bK.
1. Montrer que la fonction génératrice des moments (aussi appelée transformée de Laplace) de
X est égale à
sinh n2 t

1 a+b
h i
tX
MX (t) := E e = exp t .
n 2 sinh 1 t 2
2. En faisant un développement limité de MX à l’ordre 2 en 0, calculer les moments d’ordre 1

2
et 2 de X (c’est-à-dire E [X] et E X 2 ), puis montrer que la variance de X vaut n 12−1 .


3. En déduire une formule permettant de calculer la somme des carrés des n premiers entiers,
pour n ≥ 1. Donner une méthode pour calculer la somme des n premiers cubes.
Remarque 3. On pouvait aussi utiliser la fonction génératrice pour résoudre l’exercice 2.1. En
utilisant la fonction sinus cardinal hyperbolique, définie par
sinh (x) ex − e−x
sinch (x) = = ,
x 2x
on trouve, pour U uniforme sur [a ; b],
a+b b−a

MU (t) = exp t sinch t .
2 2
Exercice 2.4.
On considère deux entiers a et b avec a < b. En utilisant la fonction rand, écrire une fonction
runif_int(m, n, a, b) qui renvoie une matrice de taille m × n dont les coefficients sont des
réalisations indépendantes de la loi uniforme sur {a, a + 1, . . . , b}.
3 Lois discrètes
Pour continuer, nous allons simuler les principales lois de référence sur N et en rappeler les
définitions et les principales propriétés.
Définition 4 (Bernoulli, binomiale, géométrique). Soit p ∈ [0 ; 1], n ∈ {1, 2, . . . } et X une
variable aléatoire. On dit que X suit :
1. la loi de Bernoulli de paramètre p (notation : X ∼ B (1, p) ) si
P (X = 1) = 1 − P (X = 0) = p ;
2. la loi binomiale de paramètres n et p (notation : X ∼ B (n, p)) si X a la même loi que

Y1 + Y2 + · · · + Yn ,
où Y1 , Y2 , . . . , Yn , sont n variables de Bernoulli indépendantes de même paramètre p ;
3. la loi géométrique de paramètre p si X a la même loi que
inf {n ≥ 1 | Yn = 1} ,
où (Yn )n≥1 est une famille de variables de Bernoulli indépendantes de même paramètre p.
Exercice 3.1. Loi binomiale
1. Soient X et Y deux variables aléatoires indépendantes de carré intégrable. Montrer que
Var [X + Y ] = Var [X] + Var [Y ] et que pour tout réel a, Var [aX] = a2 Var [X].
2. Soit X une variable aléatoire de loi B (n, p), où n est un entier naturel non nul et p ∈ [0 ; 1].
Démontrer les formules suivantes :
!
n k
P (X = k) = p (1 − p)n−k , k = 0, 1, . . . , n, (1)
k
E [X] = np, (2)
Var [X] = np (1 − p) . (3)
Remarque 5. La réciproque de la première formule de la question 1 est, sauf cas très particulier
(voir les vecteurs gaussiens) complètement fausse. Pour un contre-exemple, considérer le cas
où ε est une variable uniforme sur {−1; 1} et X suit la loi normale centrée réduite. Alors, si
Y = ε × X, on peut montrer (exercice) que Y suit également la loi normale centrée réduite et
que Var [X + Y ] = Var [X] + Var [Y ] alors que X et Y ne sont pas indépendantes.

Exercice 3.2. Simulation de la loi binomiale

1. En utilisant la fonction rand, écrire une fonction rbernoulli(m,n, p) qui renvoie une ma-
trice de taille m × n dont les coefficients sont des réalisations indépendantes de la loi de
Bernoulli de paramètre p.
2. En utilisant la définition de la loi binomiale, écrire une fonction rbinomiale(m,n,r,p) qui
renvoie une matrice de taille m × n dont les coefficients sont des réalisations indépendantes
de la loi B (r, p) .
3. En utilisant la formule (1), écrire une fonction rbinomiale2(m,n,r,p) ayant le même objectif
que la fonction précédente mais n’utilisant que m × n appels à rand.
Exercice 3.3. Fonction génératrice des probabilités

Soit X une variable aléatoire à valeurs dans N et GX , sa fonction génératrice des probabilités,
définie pour s ∈ [0 ; 1] par
h i ∞
X
GX (s) = E sX = P (X = k) sk .
k=0
1. Calculer GX (1) et G0X (1). Pour le calcul de G0X (1), on pourra, pour intervertir les limites
utiliser le théorème de convergence monotone avec la mesure de comptage.
2. Soit n un entier supérieur ou égal à 1. On suppose que E [X n ] < ∞. Montrer que
(n)
E [X (X − 1) (X − 2) · · · (X − n + 1)] = GX (1) ,
(n)
où GX désigne la dérivée nème de GX .
3. Montrer que, si X est de carré intégrable, on a
2
Var [X] = G00X (1) + G0X (1) − G0X (1) .
Exercice 3.4. Loi géométrique

Soit p ∈ ]0 ; 1[ et T une variable aléatoire de loi géométrique de paramètre p. Soit q = 1 − p.
1. En utilisant la définition de la loi géométrique, démontrer que
P (T = k) = q k−1 p, k = 1, 2, . . . . (4)
2. Montrer que la fonction génératrice des probabilités de T vaut

ps
GT (s) = , s ∈ [0 ; 1]. (5)
1 − qs
3. Démontrer que l’espérance et la variance de T sont données par

1 q
E [T ] = et Var [T ] = . (6)
p p2
Exercice 3.5. Loi géométrique avec chaîne de Markov

Soit p ∈ ]0 ; 1[ et (Xk )k≥0 la chaîne de Markov dont l’ensemble d’états est {0; 1} telle que :
X0 = 0 p.s. et pour tout k ≥ 0,
P (Xk+1 = 1 | Xk = 0) = p
P (Xk+1 = 0 | Xk = 0) = 1 − p
P (Xk+1 = 1 | Xk = 1) = 1
Soit T le temps d’atteinte de l’état 1. On note q = 1 − p.

1. Montrer que T suit la loi géométrique de paramètre p.

2. En utilisant la propriété de Markov faible au temps 1, montrer que pour tout entier naturel
non nul k,
q k−1
!
h i X k h i
k
E T =1+ E Ti .
p i=0 i
3. Retrouver E [T ] et Var [T ].
Exercice 3.6. Simulation de la loi géométrique
Soit p ∈ ]0 ; 1[ et q = 1 − p .
1. En utilisant la définition, écrire une fonction rgeom(m,n,p) qui renvoie une matrice de taille
m × n dont les coefficients sont des réalisations indépendantes de la loi géométrique de para-
mètre p.
2. Soit U une variable aléatoire de loi uniforme sur ]0 ; 1[. Montrer que la variable aléatoire
ln (U )

+ 1,
ln (q)
où bxc désigne la partie entière du réel x, suit la loi géométrique de paramètre p.
3. En utilisant ce résultat, écrire une fonction rgeom2(m,n,p) ayant le même rôle que la fonction
précédente.
Définition 6. Soit λ > 0. On dit que la variable aléatoire N suit la loi de Poisson lorsque, pour
tout k ∈ N,
λk
P (N = k) = e−λ .
k!
Exercice 3.7. Fonction génératrice des probabilités de la loi de Poisson
1. Montrer que la fonction génératrice des probabilités de N est
GN (s) = eλ(s−1) .
2. En déduire l’égalité
E [N ] = Var [N ] = λ.

3. Soit Xn de loi B n, nλ , pour n assez grand. Montrer que Xn converge en loi vers N .
Exercice 3.8. Simulation de la loi de Poisson
La loi Γ (α, λ), où α et λ sont deux réels strictement positifs est la loi de densité
xα−1 α −λx
fα,λ (x) = 1[0;∞[ λ e .
Γ (α)
On rappelle que la fonction Γ est définie sur ]0 ; ∞[ par
Z ∞
Γ (x) = tx−1 e−t dt
0
et que, si n ∈ N∗ , Γ (n) = (n − 1) !. On fixe λ > 0.
1. Soient X1 , X2 , . . . des variables aléatoires indépendantes de même loi exponentielle de para-
mètre λ. Montrer que, pour tout n ∈ N∗ , la variable aléatoire X1 + X2 + · · · + Xn suit la loi
Γ (n, λ).
2. Soit (Uk )k≥1 une famille de variables aléatoires indépendantes de loi uniforme sur [0 ; 1]. Soit
n o
N = max n ≥ 1 U1 U2 · · · Un > e−λ ,
avec ici max ∅ = 0. Justifier que P (N ≥ n) = P (X1 + X2 + · · · + Xn ≤ 1).

3. En déduire que N suit la loi de Poisson de paramètre λ .
4. Écrire une fonction rpoisson(m,n, lambda) qui simule la loi de poisson de paramètre λ .

4 Inversion de la fonction de répartition

Nous passons maintenant à une méthode générale, simple et puissante pour simuler des lois
de probabilité sur R. Mais d’abord, quelques rappels sur la fonction de répartition sont de mise.
4.1 Fonction de répartition

Soit X une variable aléatoire à valeurs réelles. La fonction de répartition de X (en anglais,
cumulative distribution function) est définie sur R par
FX (x) = P (X ≤ x) .
En Scilab, la plupart des fonctions de répartitions des lois de références sont déjà programmées.
Elles sont en général préfixées par cdf, par exemple, cdfnor, cdfpoi, ... Voir l’exemple ci-dessous
(en particulier la bizarrerie consistant à entrer une chaîne de caractère « PQ » pour dire que l’on
veut la fonction de répartition).
// Exemple d ' utilisation de cdfnor

// pour la fonction de r é partition de la loi normale
mu = 0 ; sigma = 1;
// Calcul de F (1 ,96)
cdfnor ( " PQ " , 1.96 , mu , sigma )
// r é sultat ( bien connu ) : 0.9750021
Exercice 4.1. Rappels sur les fonctions croissantes et la fonction de répartition

1. Soit I un intervalle réel d’intérieur non vide et f : I → R une fonction croissante. Montrer
qu’en tout point x de l’intérieur de I, f admet une limite à gauche f (x−) et une limite à
droite f (x+), lesquelles vérifient
f (x−) ≤ f (x) ≤ f (x+) .
2. En déduire que l’ensemble des points de discontinuité de f est au plus dénombrable.

3. Soit X une variable aléatoire à valeurs réelles et F sa fonction de répartition. Montrer que F
est croissante, a pour limite 0 en −∞, 1 en +∞, admet des limites à gauche et est continue
à droite en tout réel.
4. Justifier que tout ouvert de R est une union dénombrable d’intervalles ouverts.
5. Soient a < b. Exprimer avec la fonction F la probabilité P (X ∈ ]a ; b[).
6. Déduire des deux points précédents que la fonction F caractérise la loi de X.
4.2 Inverse de Lévy

Notons maintenant, pour 0 < u ≤ x,
FX−1 (u) = inf {x ∈ R | FX (x) ≥ u} .
Cette fonction s’appelle l’inverse généralisée ou l’inverse de Lévy de FX ou encore fonction

quantile de X.
En Scilab, on utilise encore les fonctions cdf avec des paramètres différents. Voir l’exemple
ci-dessous (en particulier la bizarrerie consistant à entrer une chaîne de caractère « X » pour
dire que l’on veut l’inverse de la fonction de répartition et l’entrée des deux nombres p et 1-p).

// Inversement , je veux connaitre x

// tel que F ( x ) = 0.975
mu = 0 ; sigma = 1;
p = 0.975; q = 1 - p ;
cdfnor ( " X " , mu , sigma , p , q )
// r é sultat : 1.959964
Consulter l’aide et s’entraîner avant l’épreuve semble indispensable.

Théorème 7 (Inverse de la fonction de répartition). La fonction F −1 est croissante et continue
à gauche. De plus, si U est une variable aléatoire de loi uniforme sur [0 ; 1]. Alors la variable
aléatoire FX−1 (U ) a la même loi que X.
Exercice 4.2. Démonstration du théorème
Pour alléger les notations, on écrit F := FX et F −1 := FX−1 .
1. Montrer que F −1 est croissante.
2. Montrer que pour tout u ∈ ]0 ; 1[,

F F −1 (u) ≥ u. (7)
3. Montrer que F −1 est continue à gauche.

4. Montrer que pour tous x ∈ R et 0 < u ≤ 1,
u ≤ F (x) ⇐⇒ F −1 (u) ≤ x.
5. Calculer P F −1 (U ) ≤ x puis conclure.

4.3 Loi prenant un nombre dénombrable de valeurs

Ceci n’entre pas exactement dans le cadre précédent mais suit le même principe. Soit X une
variable aléatoire à valeurs dans un ensemble (fini ou) dénombrable
E = {x1 , x2 , . . . } .
Soit, pour tout k ≥ 1, pk = P (X = xk ). Soit la suite (sk )k≥0 définie par s0 = 0 et, pour tout
k ≥ 0,
sk+1 = sk + pk . 1
Alors, si U suit la loi uniforme sur [0 ; 1], la variable aléatoire Z définie par
Z = xk ⇐⇒ sk < U ≤ sk+1 ,
a la même loi que X. On remarquera que cette méthode fait intervenir un nombre aléatoire de
comparaisons. On a parfois intérêt à indexer les valeurs xk par ordre décroissant de probabilités
pour diminuer le nombre de tests.
En Scilab, on pensera dans ce genre de situation à utiliser les fonctions cumsum (pour la
somme cumulative) et dsearch (qui localise la classe d’une valeur, les classes étant données par
un vecteur trié).
Exercice 4.3. 1. En utilisant la méthode ci-dessus, proposer une nouvelle fonction pour la
simulation de la loi de Poisson.
2. Soit X une variable aléatoire à valeurs dans {0, 1, 2, . . . , N − 1}. Écrire une fonction en Scilab
rfinite(m,n,P) qui renvoie un échantillon de taille m × n suivant la loi de X, où la matrice
P contient les valeurs pk = P (X = k), pour k = 0, 1, . . . , N − 1.
1. Bref, les sk sont les sommes cumulatives des pk .

4.4 Application à la loi de Cauchy et la loi exponentielle

Exercice 4.4. 1. Calculer l’inverse de la fonction de répartition de la fonction exponentielle de
paramètre λ > 0, de densité
λ exp (−λx) .
2. En déduire une fonction rexp(m,n,lambda) pour simuler la loi exponentielle.
3. Faire de même pour la loi de Cauchy centrée réduite, de densité
1
f (x) = .
π (1 + x2 )
4. La loi de Cauchy admet-elle une espérance ?
5 Méthode du rejet
Ce chapitre s’inspire, plus encore que les autres, du livre [BC07] (partie 1.10).
5.1 Loi uniforme

Proposition 8 (Méthode du rejet). Soit µ une mesure sur un espace mesurable (E, A) et B ⊂ A
deux mesurables tels que
0 < µ (B) ≤ µ (A) < ∞.
Soit p = µ(B)
µ(A) . Soit (Uk )k≥1 une suite de variables aléatoires i.i.d. de loi (µ-)uniforme sur A
et Vk = 1{Uk ∈B} . Alors les variables aléatoires Vk sont i.i.d. de loi de Bernoulli de paramètre
p. Si de plus, T = inf {k ≥ 1 | Yk = 1}, alors T suit la loi géométrique de paramètre p et est
indépendante de XT , laquelle suit la loi uniforme sur B.
Démonstration. Par définition de la loi uniforme sur A, on a pour tout k ≥ 1,
µ (B)
P (Xk ∈ B) = ,
µ (A)
donc Yk suit la loi de Bernoulli de paramètre p et T suit la loi géométrique de paramètre p. Pour
vérifier la deuxième partie du théorème, on calcule, étant donné une partie mesurable C incluse
dans B, et un entier n, en utilisant l’indépendance des Xk :
P (T = n, XT ∈ C) = P (X1 ∈
/ C, . . . , Xn−1 ∈
/ C, Xn ∈ C)
= P (X1 ∈/ C) · · · P (Xn−1 ∈/ C) P (Xn ∈ C)
µ (C) µ (C) µ (B)
= (1 − p)n−1 = (1 − p)n−1
µ (A) µ (B) µ (A)
µ (C)
= (1 − p)n−1 p = P (T = n) P (Z ∈ C) ,
µ (B)
où Z suit la loi uniforme sur B.
L’exercice suivant, très connu (bien que peu efficace), permet de calculer π en utilisant la
méthode de Monte-Carlo.
Exercice 5.1. Calcul de π par la méthode de Monte-Carlo
1. En utilisant la loi uniforme sur le carré C = [0 ; 1]2 , simuler la loi uniforme sur le quart de
disque n o
Q = (x, y) ∈ C x2 + y 2 ≤ 1 .
Combien faut-il, en moyenne, simuler de variables uniformes C pour obtenir une uniforme
sur Q ?

2. Soient U1 , U2 , . . . i.i.d. de lois uniformes sur C. Justifier que

1 π
Card {1 ≤ k ≤ n | Uk ∈ Q} −−−→ .
n n→∞ 4
3. Écrire un programme utilisant cette méthode pour calculer π.
4. En utilisant l’approximation donnée par le théorème central limite, estimer le nombre minimal
de lancers à faire pour que la valeur donnée de π4 soit juste à 10−2 près, avec une probabilité
supérieure à 95%.
5.2 Application aux lois à densités

Pour une approche plus intuitive du théorème suivant, se référer à [BC07], partie 1.11.
Théorème 9 (Méthode du rejet comparatif). Soient µ et ν deux lois à densité sur R, de
fonctions de densité respectives f et g. On suppose qu’il existe une constante c > 0 telle f ≤ cg.
Soient X1 , X2 , . . . des variables i.i.d. de loi µ et soient U1 , U2 , . . . des variables i.i.d. de loi
uniforme sur [0 ; 1], indépendantes des Xk . Soit
f (Xk )

T := inf k ≥ 1 Uk < . (8)
cg (Xk )
1
Alors, T suit la loi géométrique de paramètre c et est indépendante de XT , laquelle a pour loi
ν.
Démonstration. On commence par remarquer que l’expression 8 a un sens avec probabilité 1,
car Z
P (g (Xk ) = 0) = 1{g(x)=0} g(x) dx = 0.
f (Xk )
Les événements Uk < cg(Xk ) sont clairement indépendants de même probabilité égale à
f (Xk ) f (Xk ) f (x) 1 1
Z Z
E P Uk < Xk =E = g(x) dx = f (x) dx = .
cg (Xk ) cg (Xk ) cg(x) c c
Cela montre que T suit la loi géométrique de paramètre 1c . Ensuite, soit A un borélien et un
entier n ≥ 1. On a (en considérant X et U indépendantes de loi respectives µ et U ([0 ; 1]) pour
alléger les notations) :
P (T = n, XT ∈ A) = P (T = n, Xn ∈ A) = P (T = n) P [Xn ∈ A | T = n]
f (X)

= P (T = n) P X ∈ A U ≤
cg (X)

f (X)
P X ∈ A, U ≤ cg(X)
= P (T = n)
f (X)

P U≤ cg(X)
Z f (x) Z
cg(x)
= P (T = n) × c g(x) dxdu.
u=0 x∈A
où, pour passer de la première à la deuxième ligne, on a décomposé l’événement T = n en
f (X1 ) f (X2 ) f (Xn−1 ) f (Xn )

{T = n} = U1 ≥ , U2 ≥ , . . . , Un−1 ≥ , Un < ,
g (X1 ) g (X2 ) g (Xn−1 ) g (Xn )
et l’on a remarqué que tous ces événements à l’exception du dernier étaient indépendants de
l’événement Xn ∈ A. En utilisant le théorème de Fubini et en intégrant d’abord suivant u la
dernière intégrale devient
f (x) 1 1
Z Z
g (x) dx = f (x) dx = ν (A) ,
x∈A cg (x) x∈A c c
ce qui achève cette démonstration.

L’exercice suivant montre une application (un peu trop facile et maladroite) de cette méthode.
Exercice 5.2. On considère la fonction f définie par :

(
3x2 si x ∈ [0 ; 1],
f (x) =
0 sinon.
1. Justifier que f est une fonction de densité de probabilité, calculer sa fonction de répartition
et l’inverse de sa fonction de répartition. En déduire un programme inv_rep permettant de
simuler la loi de densité f .
2. Expliquer comment, en utilisant la méthode du rejet comparatif, on peut simuler la loi de
densité f puis mettre en œuvre cette méthode dans une fonction rejet.
3. Pour un grand nombre de valeurs simulées, combien d’appels à la fonction rand votre pro-
gramme rejet fait-il en moyenne ?
6 Loi normale, vecteurs gaussiens

Sur ce sujet, on pourra consulter par exemple les premières pages de [LG12], le chapitre 6
de [BC07], le chapitre 7 de [Can13] et bien d’autres références encore.
6.1 Variables gaussiennes

Proposition 10. Soit f : R → R définie par
1 x2
f (x) = √ e− 2 .
2π
Alors f est une densité de probabilité.
R − x2
Exercice 6.1. Soit I = e 2 dx. Justifier que
x2 +y 2
ZZ
2
I = e− 2 dxdy.
Faire un changement de variable en coordonnées polaires, puis calculer I.
Définition 11. La loi normale centrée réduite, notée N (0, 1) sur R est la loi de densité f .
Si µ ∈ R et σ ≥ 0, on définit la loi normale N µ, σ 2 comme la loi de la variable aléatoire
X = µ + σN , où N suit la loi N (0, 1). Lorsqu’une variable aléatoire X suit une loi normale, on
dit que cette variable aléatoire est gaussienne.
Proposition 12. Soit N ∼ N (0, 1), m ∈ R, σ ≥ 0 et X ∼ N m, σ 2 .

1. L’espérance et la variance de N valent respectivement 0 et 1.

2. L’espérance et la variance de X valent respectivement m et σ 2 .
3. Si σ = 0, alors X a pour loi la mesure de Dirac en m, sinon X a pour densité la fonction
1 (x−m)2
x 7−→ √ e− 2σ2 .
2πσ
4. Les fonctions caractéristiques de N et de X valent respectivement

t2
h i
ϕN (t) := E eitN = e− 2 , (9)
2 2
− σ 2t
h i
ϕX (t) := E eitX = eimt e . (10)

5. Les fonctions génératrices des moments de N et de X sont finies sur R et valent respective-
ment
t2
h i
MN (t) := E etN = e 2 (11)
σ 2 t2
MX (t) = etm e 2 . (12)
6. La loi normale centrée réduite admet des moments de tous ordre. Ses moments d’ordre impair
sont tous nuls, tandis que ses moments d’ordre pair sont donnés par
h i (2n)!
E N 2n = .
2n n!
Exercice 6.2. Démonstration des résultats.
On indique pour chacun des points une stratégie de démonstration.
1. Pour l’espérance, utiliser un argument de symétrie. Pour la variance, faire une intégration
par parties.
2. Surtout, ne pas refaire le même calcul, utiliser la définition.
3. Considérer une fonction f mesurable positive, écrire
E [f (X)] = E [f (m + σN )] ,
et faire un changement de variable dans l’intégrale.

4. Commencer par remarquer que ϕN est à valeurs réelles pour des raisons de symétries et donc
que
x2
Z
ϕN (t) = cos (tx) e− 2 dx.
Dériver sous le signe intégrale puis intégrer par parties pour montrer que ϕN vérifie l’équation
différentielle y 0 = −ty puis conclure.
5. Faire un calcul direct (avec un changement de variables sous la forme d’une translation).
6. Utiliser, par exemple, un développement limité de la fonction caractéristique au voisinage de
0.
Corollaire 13. Toute combinaison linéaire de variables gaussiennes indépendantes est encore
gaussienne.
Exercice 6.3. Démontrer ce corollaire en utilisant la fonction caractéristique.
Le rôle central de la loi normale, tant en probabilités qu’en statistiques, vient de son ca-
ractère universel donné par le théorème central limite. La démonstration utilise les fonctions
caractéristiques. Elle n’est pas très difficile et est bien détaillée dans [Tou99].
Théorème 14 (théorème central limite). Soit µ une mesure de probabilité sur R, de variance
finie σ 2 et de moyenne m. Soient X1 , X2 , . . . des variables aléatoires i.i.d. de loi µ et pour tout
n ∈ N∗ Sn = X1 + · · · + Xn . Alors on a la convergence en loi
√ Sn

(d)

n − m −−−→ N 0, σ 2 .
n n→∞

6.2 Méthode de Box-Muller

Proposition 15 (Box, Muller). Soient Y et θ deux variables aléatoires indépendantes, avec Y
de loi exponentielle de paramètre 12 et θ uniforme sur [0 ; 2π]. Soit
√ √
(X, Y ) = Y cos θ, Y sin θ .
Alors, les variables aléatoires X et Y sont indépendantes, de même loi N (0, 1).
Exercice 6.4.
1. Démontrer la proposition précédente. On pourra considérer une fonction f mesurable positive,
et calculer h √ √ i
E f Y cos θ, Y sin θ
en faisant un premier changement de variable y = r2 , puis un passage en coordonnées polaires.

2. En déduire une fonction Scilab rnorm(m,n) pour simuler la loi normale centrée réduite.
Remarque 16. La méthode de Box-Muller est assez classique et utilisée en pratique, mais on
lui préfère souvent la méthode de la ziggurat, inventée par Marsaglia (un grand nom de la
simulation) et ses collaborateurs. Elle est basée sur une méthode de rejet très astucieuse. Voir
l’article original [MT+ 00].
6.3 Vecteurs aléatoires

On considère toujours un espace probabilisé (Ω, F, P). Un vecteur aléatoire X = (X1 , . . . , Xn )T
est une fonction mesurable de (Ω, F, P) vers un espace (Rn , B (Rn )). Ses composantes sont des
variables aléatoires et réciproquement, la donnée de n variables aléatoires définit un vecteur
aléatoire dans Rn .
On dit que le vecteur aléatoire X est carré-intégrable, lorsque chacune de ses composantes
l’est, c’est-à-dire lorsque h i
∀i ∈ J1 ; nK, E Xi2 < ∞.
Étant donné un vecteur aléatoire X carré-intégrable, on définit l’espérance de X comme le

vecteur (déterministe) de Rn ,
E [X] = (E [X1 ] , E [X2 ] , . . . , E [Xn ])T ,
et la matrice de covariance de X par
Cov (X) = (Cov (Xi , Xj ))1≤i,j≤n ,
où la covariance des deux variables aléatoires Xi et Xj est définie par
Cov (Xi , Xj ) = E [(Xi − E [Xi ]) (Xj − E [Xj ])] = E [Xi Xj ] − E [Xi Xj ] .
La covariance est ici bien définie car, d’après l’inégalité de Cauchy-Schwarz 2 ,

h i h i
|E [Xi Xj ]|2 ≤ E Xi2 E Xj2 < ∞.
Bien sûr, Cov (Xi , Xi ) = Var [Xi ]. Noter que la covariance est, comme la variance, invariante
par translation.
Voyons maintenant l’effet du produit par une matrice A.
2. Hermann Schwarz (1843 - 1921), mathématicien allemand, auteur également du théorème de Schwarz sur
les dérivées partielles, à ne pas confondre avec Laurent Schwartz (1915 - 2002), mathématicien français célèbre
pour ses travaux sur la théorie des distribution.

Proposition 17. Soit X un vecteur aléatoire carré-intégrable de Rn et A une matrice réelle de

taille m × n. Alors AX est un vecteur aléatoire de Rm d’espérance et de covariance
E [AX] = AE [X] et Cov (AX) = A Cov (X) AT .
Démonstration. La première égalité est très facile en utilisant la linéarité de l’espérance. Pour
la seconde, on peut, quitte à lui soustraire son espérance, supposer que E [X] est le vecteur nul,
sans changer ni Cov (X), ni Cov (AX) (par invariance par translation). Notons
A = (aij )1≤i≤m, 1≤j≤n et D = (dij )1≤i,j≤m = Cov (AX) .
On a alors, pour tous i, j ∈ J1 ; mK,

n n
" ! !#
X X X
dij = E air Xr ajs Xs = air E [Xr Xs ] ajs ,
r=1 s=1 1≤r,s≤n
qui est précisément le coefficient d’indice (i, j) dans la matrice
A Cov (X) AT .
Proposition 18. Une matrice de covariance est toujours symétrique et (semi-définie) positive.
Démonstration. La symétrie est évidente. Soit u un vecteur colonne de Rn . D’après la proposition

précédente, on a
uT Cov (X) u = Cov uT X = Var uT X ≥ 0.
Comme dans le cas des variables aléatoires (réelles), les vecteurs aléatoires (de Rn ) ont une
loi caractérisée par leurs fonctions caractéristiques (leurs transformées de Fourier), définies pour
tout t = (t1 , t2 , . . . , tn ) ∈ Rn par h i
ϕX (t) = E eiht,Xi ,
où h·, ·i désigne le produit scalaire canonique sur Rn .
Avec les mêmes notations que précédemment, on obtient immédiatement,

ϕAX (t) = ϕX AT t . (13)
6.4 Vecteurs gaussiens

Définition 19. Soit n un entier naturel non nul. Un vecteur gaussien de Rn est un vecteur
aléatoire X de Rn , tel que toute combinaison linéaire de ses composantes est une variable gaus-
sienne.
Commençons par donner quelques conséquences immédiates de cette définition.
Proposition 20. Soit n un entier naturel non nul et X un vecteur aléatoire de Rn . Les propo-
sitions suivantes sont équivalentes.
1. Le vecteur X est un vecteur gaussien.
2. Pour tout vecteur (déterministe) u de Rn , le produit scalaire hu, Xi est une variable gaus-
sienne.
3. Pour tout entier naturel non nul m, toute matrice réelle A de taille m × n et tout vecteur
déterministe b de Rm , le vecteur aléatoire de Rm AX + b est un vecteur gaussien.
Exemple 21. Si X1 , X2 , . . . , Xn sont n variables gaussiennes indépendantes, alors X = (X1 , . . . , Xn )

est un vecteur gaussien d’après le corollaire 13. D’après la proposition précédente, et en utilisant
les mêmes notations, il en est donc de même pour tous les vecteurs aléatoire de la forme AX + b.

Nous allons voir qu’en fait, tous les vecteurs gaussiens peuvent s’écrire de cette façon.
Remarque 22. Si X est un vecteur gaussien, alors ses composantes sont des variables gaussiennes.
La réciproque est fausse ! La remarque 5 fournit un contre-exemple.
Le calcul de la fonction caractéristique d’un vecteur gaussien est immédiat, mais d’une im-
portance cruciale.
Proposition 23 (fonction caractéristique d’un vecteur gaussien). Soit X un vecteur aléatoire
de Rn . Alors, X est un vecteur gaussien si et seulement si, il existe un vecteur déterministe µ
de Rn et une matrice carrée de taille n, C symétrique et positive tels que que pour tout vecteur
v de Rn ,
1
h i
ϕX (v) := E eihv,Xi = eihv,µi e− 2 hv,Cvi .
Dans ce cas, µ est l’espérance de X et C est sa matrice de covariance.
Démonstration. Pour l’implication directe, on pose µ = E [X] et C = Cov (X). La variable
aléatoire hv, Xi est gaussienne, d’espérance hv, µi et de variance hv, Cvi, d’après la proposition 17
(le réel hv, Cvi est bien positif car C est symétrique positive). D’après la formule 10, sa fonction
caractéristique est
1 2
ϕhv,Xi (t) = eithv,Xi e− 2 t hv,Cvi .
On remarque pour conclure que
ϕX (v) = ϕhv,Xi (1) .
Pour l’implication réciproque, on commence par considérer un vecteur v de Rn . Comme
1 2
ϕhv,Xi (t) = ϕX (tv) = eithv,µi e− 2 t hv,Cvi
,
la variable hv, Xi est gaussienne. Ceci étant valable pour tout vecteur v, X est bien un vecteur
gaussien.
Pour calculer son espérance, il suffit d’appliquer la formule précédents aux éléments e1 , . . . , en
de la base canonique de Rn . Enfin pour trouver sa matrice de covariance, on peut considérer les
vecteurs ei + ej et utiliser la symétrie de la matrice C.
Corollaire 24. 1. La loi d’un vecteur gaussien est entièrement caractérisée par son espérance
et sa matrice de covariance.
2. Si X = (X1 , . . . , Xn ) est un vecteur gaussien, alors les variables gaussiennes X1 , . . . , Xn sont
indépendantes si et seulement si la matrice de covariance de X est diagonale, c’est-à-dire si
et seulement si les variables (X1 , . . . , Xn ) sont deux-à-deux non corellées.
Proposition 25. Soit C une matrice symétrique positive de taille n et µ un vecteur déterministe
de Rn . Soit N un vecteur gaussien de Rn dont toutes les composantes sont indépendantes de loi
normale centrée réduite.
Alors il existe des matrices A telles que AAT = C et pour toute telle matrice, le vecteur
aléatoire AN + µ est un vecteur gaussien de matrice de covariance C et d’espérance µ.
Remarque 26. En vertu des deux résultats précédent, on peut parler, pour tout vecteur µ et toute
matrice symétrique positive C de « la loi normale de matrice de covariance C et d’espérance
µ ». On notera cette dernière N (µ, C).
Démonstration. La seule chose à montrer est l’existence de telles matrices A. Une première
possibilité est d’utiliser le fait que C, symétrique, est diagonalisable par le groupe orthogonale,
donc il existe Q orthogonale et ∆ diagonale telles que
C = Q−1 ∆Q = QT ∆Q.
Comme C est positive par hypothèse, tous √ les éléments de ∆ sont positifs, de sorte que la matrice
formée par leurs racines carrées, notée ici ∆ est bien définie. On peut alors poser
√
A = QT ∆.

Remarque 27. Nous verrons plus loin avec la décomposition de Cholesky qu’il n’est nullement
besoin de diagonaliser C pour trouver une telle matrice A.
Nous justifions ensuite la remarque qui suit l’exemple 21.
Proposition 28. Soit X un vecteur gaussien de Rn , de moyenne µ et de matrice de covariance

C. Soit r ≤ n le rang de la matrice C. Alors il existe r variables gaussiennes indépendantes
N1 , N2 , . . . , Nr indépendantes centrées et réduites et une matrice A de taille n × r telles que
X = A(N1 , N2 , . . . , Nr )T + µ.
Démonstration. Quitte à raisonner sur X − µ, on peut supposer que µ est le vecteur nul. Soient
λ1 ≥ λ2 ≥ · · · ≥ λr > 0 = λr+1 = · · · = λn
les valeurs propres de C comptées avec leurs multiplicités. Soient

 1 
λ1 − 2
 
λ1
 ..  
.. 

 . 


 . 

1
λr − 2
 
 λr   
∆= et D= .
 
0

1
   
 
..
 
..
 
.
 
.
   
 
0 1
Soit Q une matrice orthogonale telle que C = QT ∆Q. Le vecteur gaussien
N := DQX
a pour matrice de covariance

 
1
 .. 

 . 

T T T
 1 
Cov (N ) = DQCQ D = DQQ ∆QQ D = D∆D =  .
 
 0 
..
 
.
 
 
0
Par conséquent ses r premières composantes N1 , . . . , Nr sont bien gaussiennes centrées réduites
et indépendantes, d’après le corollaire 24, tandis que ses n−r dernières composantes sont presque
sûrement égales à 0.
Pour finir, on note que
√ 
λ1
 .. 
 .  
N1
 √ 
−1 −1

−1  λr 
  .. 
X=Q D N =Q   . .
 0 
 N
..

  r
 . 
0
Remarque 29. En Scilab, la fonction spec permet de diagonaliser une matrice. Ce n’est pas
totalement explicite dans la documentation, mais dans le cas d’une matrice symétrique, la base
de vecteurs propres renvoyée est bien orthogonale 3 .
3. La sous-routine fortran utilisée dans ce cas est DSYEV, dont la documentation est claire.

6.5 Décomposition de Cholesky

Comme promis plus haut, la décomposition de Cholesky permet d’obtenir une matrice C,
symétrique et positive comme un produit AAT avec la propriété très pratique que A est triangu-
laire. Étonnamment, il est plus simple de démontrer d’abord l’existence d’une telle décomposition
avant d’en donner l’algorithme.
Lemme 30 (Complément de Schur). Soit C une matrice symétrique définie positive d’ordre
n = p + q dont une écriture par blocs est
!
C11 C12
C=
C21 C22
Avec C11 carrée d’ordre p et C22 carrée d’ordre q. Alors C11 est inversible et le complément de
Schur
−1
S = C22 − C21 C11 C12
est une matrice carrée d’ordre p, symétrique définie positive.
Exercice 6.5. Démonstration du lemme 30
1. Soit (e1 , . . . , en ) la base canonique de Rn et E = Vect (e1 , . . . , ep ). Soit q la forme quadratique
associée à la matrice C. Quelle est la matrice associée à la restriction de q à E ×E ? En déduire
que C11 est inversible.
2. Soit x un vecteur non nul du sous-espace E. Soit y le vecteur de Rn dont une écriture par
blocs est !
−1
−C11 C12 x
y=
x
Calculer q (y). Conclure.
Théorème 31 (Décomposition de Cholesky). Soit C une matrice symétrique définie positive
d’ordre n. Alors il existe une unique matrice triangulaire inférieure R telle que
C = RRT .
Démonstration. On commence par démontrer l’existence, par récurrence sur la dimension n. Le
résultat est évident si n = 1. Supposons que n ≥ 2 et que le théorème est vrai pour les matrices
d’ordre n − 1. Notons !
c11 C12
C= ,
C21 C22
où c11 est un scalaire (strictement positif), C12 est une ligne de taille n−1, C21 = C12 T et C22 est
carrée d’ordre n − 1. Soit S = C22 − C21 c−111 C12 le complément de Schur associé à cette écriture
par blocs. D’après le lemme précédent, la matrice S est encore définie positive. Par hypothèse
de récurrence, il existe une matrice triangulaire inférieure RS telle que S = RS RS T .
Il suffit de maintenant de poser (toujours en écriture par blocs)
√ !
c11 0
R = √1
c11 C21 RS
pour avoir C = RRT .

Pour l’unicité, on écrit R = (rij )1≤i,j≤n et l’on déduit de RRT = C les égalités, pour
1≤i≤j≤n:
i
X
cij = rik rjk (14)
k=1
Xi
2
cii = rii . (15)
k=1

√ c1j
On calcule donc d’abord le terme r11 = c11 , puis ceux de la première colonne rj1 = r11 . Puis,
si les i − 1 premières colonnes ont été calculées, on calcule la colonne d’indice i par :
v
u
u i−1
X
rii = tcii − 2
rik
k=1
P
i−1
cij − k=1 rik rjk
rji = .
rii
Exercice 6.6.
1. Vérifier que C = RRT dans la démonstration précédente.
2. Mettre en œuvre la décomposition de Cholesky dans une fonction Scilab.
3. Tester votre fonction.
4. Ne plus jamais utiliser votre fonction, utiliser la fonction prédéfinie chol, bien plus rapide
(mais lire la documentation de cette fonction).
6.6 Simulation des vecteurs gaussiens

Tout a déjà été dit dans les sous-parties précédentes. Si l’on veut :
— simuler la loi N (µ, C), où C est symétrique définie positive d’ordre n : on utilise la dé-
composition de Cholesky pour obtenir une matrice A telle que C = AAT , puis on utilise
un tirage N de n variables gaussiennes centrées réduites. En Scilab, on peut donc utiliser
les commandes rand(n,1,"normal") ou bien grand(n,1,"nor",0,1) ;
— diagonaliser la matrice C dans le groupe orthogonal : on utilisera la fonction spec (avec
deux arguments de retour).
!
3 1
Exercice 6.7. Soit C = et µ = −2 −3 .
1 2
1. En utilisant la décomposition de Cholesky, écrire un script permettant d’obtenir un échan-
tillon X, de 1000 points de R2 de loi N (µ, C).
2. Compléter votre script de façon à tracer ces 1000 points dans un repère, en utilisant, par
exemple la commande plot(X(1,:), X(2,:), "r.").
3. Diagonaliser la matrice C, puis tracer (par exemple avec xarrows) les vecteurs ui = λi ei ,
pour i = 1, 2, où les λi sont les valeurs propres de la matrice Cet ei des vecteurs propres (de
norme 1) associés.
7 Quelques exemples de processus

7.1 Chaîne de Markov à temps discret et espace d’états fini
On considère un ensemble dénombrable d’états E, muni de la tribu discrète (c’est-à-dire
l’ensemble de ses parties 2E ) et une matrice stochastique P sur E, c’est-à-dire une application
P : E × E −→ [0 ; 1]
telle que X
∀x ∈ E, P (x, y) = 1.
y∈E
On dit que P est la matrice de transitions de la chaîne de Markov. Soit µ une mesure de
probabilité sur E appelée loi initiale de la chaîne.

On définit une variable aléatoire X à valeurs dans l’ensemble des suites E N muni de sa tribu
produit 4 par :
1. la loi de X0 est µ ;
2. pour tout n ≥ 0 et tous x, y ∈ E,
P (Xn+1 = y | Xn = x) = P (x, y) .
La mesure de probabilité ainsi définie sur l’espace des suites E N est notée Pµ .
En conditionnant par rapport aux différentes valeurs de X0 , on voit que la loi de X1 sous Pµ
est donnée par le produit µP :
X
∀y ∈ E, Pµ (X1 = y) = µP (y) = µ (x) P (x, y) .
x∈E
On définit les puissances de la matrice stochastique P de la façon usuelle :

X
∀x, y ∈ E, P k+1 (x, y) = P k (x, z) P (z, y) .
z∈E
Le même raisonnement que ci-dessus montre alors que la loi de Xn sous Pµ est donnée par le
produit µP n .
Supposons maintenant que l’espace d’états E est fini. À l’aide des techniques vues précé-
demment, on voit qu’à n fixé, il n’est pas difficile de simuler une trajectoire X0 , X1 , . . . , Xn de
la chaîne de Markov. On peut par exemple :
1. utiliser la méthode de l’inverse de la fonction de répartition pour simuler la loi µ pour obtenir
ainsi un certain x0 dans E ;
2. puis simuler de la même manière la loi de probabilité y 7→ P (x0 , y) pour obtenir un certain
x1 ;
3. et recommencer l’étape précédente, en remplaçant x0 par x1 , etc.
On peut bien sûr trouver de nombreuses améliorations à l’algorithme précédent (calculer toutes
les sommes cumulatives des lignes de P dès le départ, etc). Tout ceci est déjà programmé dans
Scilab, on utilise grand(n,"markov",P,x0) où X0 est presque sûrement constante
pour le cas
égale à x0 . Si l’on donne un vecteur X0 = x0 , . . . x0 comme dernier argument, alors grand←-
1 N
(n, "markov", P, X_0) renvoie N trajectoires simulées, d’états initiaux respectifs x10 , x20 , etc.
Signalons aussi la fonction genmarkov qui permet d’obtenir des matrices de transitions ayant
certaines propriétés.
Exercice 7.1. Le chat bagarreur
Pito le chat n’a que trois activités : dormir, manger et se battre avec les autres chats. À chaque
heure, il peut décider, ou non de changer d’activité. On note Xn l’activité du chat après n heures.
1. Après avoir dormi une heure, il continue à dormir avec probabilité 0, 7 ou bien va manger
avec probabilité 0, 2 et se bat avec probabilité 0, 1.
2. Après avoir mangé, il va se battre avec probabilité 0, 3 ou bien dort avec probabilité 0, 7.
3. Après s’être battu, il va manger avec probabilité 0, 5 ou dormir avec probabilité 0, 5.
On suppose qu’au début de l’observation, Pito dort.
1. Tracer le graphe, puis écrire la matrice de transition de cette chaîne de Markov.
2. Calculer numériquement avec Scilab la loi de l’activité de Pito après 5 heures.
3. Simuler avec Scilab 1000 observations d’une durée de 15 heures chacune. Donner les fréquences
empiriques associées à X15 .
4. Justifier que la chaîne de Markov a une unique probabilité invariante, la calculer et justifier
que la loi de Xn converge vers cette probabilité invariante.
4. c’est-à-dire la plus petite tribu qui rend les applications coordonnées ϕk : (xn )n≥0 7→ xk mesurables

7.2 Marche aléatoire

7.2.1 Marche aléatoire dans R
Soit η une mesure de probabilité sur R et Z1 , Z2 , . . . des variables aléatoires i.i.d. de loi η.
On définit la marche aléatoire de position initiale x0 ∈ R et de loi d’incréments η par X0 = x0
et, pour tout k ∈ N,
Xk+1 = Xk + Zk+1 .
Si l’on sait simuler la loi η, la simulation de cette marche aléatoire se fait de façon immédiate à
l’aide de la fonction cumsum.
Exercice 7.2. Soit η définie par

1
η ({−1}) = η ({1}) = ,
2
et (Xk )k≥0 la marche aléatoire de loi d’incréments η et de position initiale 0.
1. Simuler une trajectoire de X entre les temps 0 et 100.
2. On étend X à tous les temps t ∈ [0 ; +∞[ par interpolation linéaire. Tracer une trajectoire de
X.
3. Soit n ≥ 1 et t ∈ [0 ; +∞[. On définit le changement d’échelle Ytn de la marche par
Xnt
Ytn = √ ,
σ n
où σ est l’écart-type de la loi η. Tracer une trajectoire de Ytn , pour n = 1000 et t ∈ [0 ; 1].
4. À t fixé, est-ce que Ytn converge en loi lorsque n tend vers l’infini ? Si oui, vers quelle loi ?
7.2.2 Marche aléatoire dans R2

Exercice 7.3. Soit (Xk )k∈N une suite de points aléatoires du plan R2 telle que X0 = (0, 0) et
pour tout entier k,
Xk+1 − Xk ∈ {(1; 0) , (0; 1) , (−1; 0) , (0; −1)} ,
et ces 4 possibilités sont équiprobables.
Tracer avec Scilab une trajectoire de (Xk ) pour 0 ≤ k ≤ 10000.
7.3 Ruine du joueur

Soit p ∈ ]0 ; 1[ et q = 1 − p. Un joueur joue à un jeu de hasard avec une fortune initiale
a ∈ {1, 2, . . . }. À chaque tour, il gagne 1 euro avec probabilité p et perd un euro avec probabilité
q. Soit c un entier supérieur à a. Le joueur décide de continuer à jouer jusqu’à ce que sa fortune
soit égale à c ou bien qu’il soit ruiné.
On modélise cette situation par une chaîne de Markov Xn à espace d’états Z avec X0 = 0 et

p si y = x + 1,


P (Xn+1 = y | Xn = x) = q si y = x − 1,


0 sinon.
Pour x ∈ Z, on note
Tx = inf {k ≥ 0 | Xk = x} ,
avec inf ∅ = +∞. Soient a et b dans {1, 2, . . . } et
not.
τ = min {T−a , Tb } = T−a ∧ Tb .

On s’intéresse à la probabilité que le joueur s’arrête en ayant gagné b euros par rapport à sa
mise de départ :
α = P (τ = Tb ) .
On note pour tout n ≥ 0, Fn = σ (X0 , X1 , X2 , . . . , Xn ) la tribu engendrée par les variables
aléatoires X0 , X1 , . . . , Xn .
Exercice 7.4. Cas où p = q = 21 .
1. Faire en Scilab un nombre important de simulation dans le cas où a = b = 10 pour donner
une estimation de α.
2. Montrer que pour tout x ∈ Z, Tx est un temps d’arrêt par rapport à la filtration (Fn )n≥0 En
déduire que τ est également un temps d’arrêt.
3. Montrer que τ est presque sûrement fini (on pourra par exemple considérer la chaîne de
Markov a espace d’états {−a, −a + 1, . . . , b − 1, b} ayant les mêmes probabilités de transition
que la chaîne de départ sauf en b et en −a où elle est réfléchie, c’est-à-dire va respectivement
en b − 1 et en −a + 1 avec probabilité 1).
4. Montrer que (Xn )n≥0 est une martingale par rapport à la filtration Fn .
5. En déduire que (Xτ ∧n ) est une martingale. Montrer que cette dernière est uniformément
intégrable.
6. En déduire la valeur de α en fonction de a et b.
Exercice 7.5. Cas où p 6= q
On reprend les notations de l’exercice précédent.
1. Faire d’autres simulations, par exemple pour p = 0, 4 et a = b = 10.
2. Montrer que le processus, défini pour n ≥ 0 par
Xn
q
Mn =
p
est une martingale.
3. En déduire la valeur de α dans ce cas.
7.4 Processus de Galton-Watson

On commence par se donner une loi de probabilité η sur N, appelée loi de reproduction. Le
processus de Galton-Watson a été introduit au départ pour étudier la probabilité de survie des
noms de famille. Il modélise de façon assez simple l’évolution d’une population.
On considère qu’il n’y a, au temps 0, qu’un seul individu. Celui-ci se reproduit suivant une
certaine loi de probabilité appelée loi de reproduction. S’il n’a pas d’enfants, il y a extinction du
processus. S’il a des enfants, ceux-ci se reproduisent de façon indépendante, suivant la même loi
de reproduction.
On peut formaliser le modèle de la façon suivante. Soit η = (ηk )k≥0 une suite de réels positifs
de somme égale à 1. Avec un léger abus de notation, on assimile la suite η à la loi de probabilités
qu’elle définit sur N. Soit ξ une variable aléatoire de loi η, de sorte que, pour tout entier naturel
k,
P (ξ = k) = ηk .
On note, pour n dans N, Zn , le nombre d’individus à la génération n. Ainsi, Z0 = 1. On se donne
une famille (ξn,i )n≥0,i≥1 de variables i.i.d. de loi η, puis on définit par récurrence
Zn
X
Zn+1 = ξi,n n = 1, 2, . . . ,
i=1
où, bien sûre la somme sur un ensemble d’indices vide vaut 0.

Définition 32. Le processus (Zn )n≥0 ainsi défini est appelé processus de Galton-Watson de loi
de reproduction η.
Exercice 7.6. Simulation d’un processus de Galton-Watson

1. Nous commençons par simuler quelques lois de reproduction. Prenons le cas où η est à support
fini, c’est-à-dire qu’il existe un entier N tel que pour tout k ≥ 1, ηN +k = 0. Écrire une fonction
support_fini(m,n,eta) simulant une matrice aléatoire de taille m × n dont les coefficients
sont i.i.d. de loi η, et dont le paramètre eta est le vecteur (η0 , . . . , ηN ).
2. Soit maintenant η une loi de probabilité quelconque sur N. On suppose qu’on dispose d’une
fonction xi qui est telle que pour tout entier naturel n, xi(n) est un vecteur ligne dont les
coefficients sont des variables aléatoires i.i.d. de loi η.
Écrire une fonction galton_watson(n, xi) qui renvoie un échantillon du vecteur aléatoire
(Z0 , Z1 , Z2 , . . . , Zn ) de loi de reproduction η.
3. À l’aide des deux questions précédentes, faire une dizaine de simulations du processus (Z0 , . . . , Z15 )
et les tracer dans le cas
1 2 2
η0 = , η1 = , η2 = .
5 5 5
Estimer (avec un grand nombre de simulations) la probabilité P (Z15 = 0).
2
4. Soit f (s) = 1+2s+2s
5 la fonction génératrice des probabilités de la loi η. Calculer à l’aide de
Scilab le réel
f ◦ f ◦ · · · ◦ f (0) ,
| {z }
15 fois
et comparer avec le résultat obtenu à la question précedente. Pour une justification de ce

résultat, voir l’exercice suivant.
5. Reprendre (si l’envie vous en prend) tout ou partie des deux questions précédentes, dans les
cas suivants :
— η0 = 0, η1 = η2 = 1/2 ;
— η0 = 14 , η1 = 12 , η2 = 1
4 ;
— η est la loi géométrique de paramètre 1/2 ;
— η est la loi de poisson de paramètre 3.
Dégageons tout de même quelques résultats plus théoriques sur ce modèle. On garde les
mêmes notations que précédemment. Nous nous intéressons tout d’abord à la probabilité que le
processus s’éteigne, c’est-à-dire qu’il existe un entier n pour lequel Zn = 0.
Le personnage principal de cette histoire est la fonction génératrice des probabilités de la loi
η, que nous noterons f . Celle-ci est définie au moins sur l’intervalle [0 ; 1] par :
∞
X
f (s) = ηk sk .
k=0
Ses principales propriétés ont déjà été étudiées lors de l’exercice 3.3.
Lemme 33. Pour tout entier naturel n, la fonction
fn := f ◦ f ◦ · · · ◦ f
| {z }
n fois
est la fonction génératrice des probabilités de la variable aléatoire Zn .

Démonstration. Le résultat est vrai au rang n = 1, car la loi de Z1 est η. On raisonne ensuite
par récurrence. Soit n un entier naturel. Par définition, la fonction génératrice gn+1 de Zn+1
vaut h i
gn+1 (s) = E sZn+1 .
On utilise ensuite le fait que
Zn
X ∞
X k
X
Zn+1 = ξn,i = 1{Zn =k} ξn,i ,
i=1 k=1 i=1
pour obtenir en conditionnant par rapport aux valeurs possibles de Zn et en utilisant que les
variables aléatoires ξn,i , pour i entier naturel non nul, sont indépendantes et indépendantes de
Zn :
h i ∞ P
k

ξ
X
E sZn+1 = P (Zn = k) E s i=1 n,i Zn = k
k=0
X∞ k
Y h i
= P (Zn = k) E sξn,i
k=0 i=1
∞
P (Zn = k) (f (s))s
X
=
k=0
= fn (f (s)) ,
la dernière égalité découlant de l’hypothèse de récurrence car on a fait apparaître la fonction

génératrice de Zn .
A Rappels d’intégration et de probabilités : l’appendice dont

vous êtes le héros
Pour se rafraîchir la mémoire, on pourra consulter les ouvrages suivants :
[LG06] Polycopié sur l’intégration et les probabilités de Jean-François Le Gall, disponible en
ligne. Une référence.
[BP12] Ouvrage de référence pour la théorie de la mesure, l’intégration et la transformée de
Fourier aux niveaux L3 et M1. Très détaillé.
[Rud87] Pas toujours facile d’accès, mais un ouvrage d’analyse très complet et profond.
[Can13] Un excellent livre d’introduction aux probabilités, de niveau L3, très détaillé, avec des
exercices souvent abordables et corrigés.
[Dur10] Une très bonne référence en anglais qui aborde certains points parfois absents des
références précédentes.
[RS12] Le chapitre 10 liste, souvent sans démonstrations, les théorèmes classiques les plus
méconnus des candidats à l’agrégation.
[Tou99] Le premier chapitre est une bonne liste des résultats importants d’intégration et de
probabilité.
1. Dans la théorie de l’intégration de Lebesgue, il y a trois théorèmes sur l’interversion limite-
intégrale pour les suites de fonctions. Les énoncer. Donner un ordre possible de démonstration.
2. Donner un énoncé possible et une démonstration pour le théorème de continuité sous le signe
intégrale.
3. Faire de même pour le théorème de dérivabilité sous le signe intégrale.

4. Énoncer le théorème de Fubini-Tonelli et le théorème de Fubini-Lebesgue. Rappeler comment

on utilise en pratique le théorème de Fubini-Lebesgue.
5. Donner des conditions pratiques pour qu’une application C 1 entre deux ouverts de Rn soit
un C 1 -difféomorphisme.
6. Énoncer le théorème de changement de variable dans Rn pour l’intégrale de Lebesgue.
7. Énoncer et démontrer le changement de coordonnées en polaires.
8. Rappeler la définition de la convergence en loi. Donner un critère pour une suite de variables
aléatoires à valeurs discrètes. Donner un critère utilisant la fonction caractéristique. Donner
un critère utilisant la fonction de répartition.
9. Comment peut-on utiliser la fonction caractéristique pour démontrer que deux variables aléa-
toires sont indépendantes ?
10. Même question avec la fonction de répartition.
Sur la transformée de Laplace, une bonne référence est la partie 3.1 de [Tou99]. La transfor-
mée de Fourier est abondamment traitée dans [Can13] et [LG06]. Pour la fonction génératrice,
voir par exemple, [Can13] partie VI.3 et le petit paragraphe 8.2.4 dans [LG06].
B Solutions des exercices

Corrigé de l’exercice 1.1 Le générateur de nombre aléatoire par défaut de Scilab rand a un
état totalement décrit par un entier appelé la graine. Par défaut, Scilab initialise la graine à 0.
Ainsi, chaque fois que l’on redémarre le programme, on obtient exactement les mêmes résultats.
Contrairement à ce que l’on pourrait penser, c’est un très bon comportement par défaut car il
permet de déboguer les simulations en permettant la reproductabilité des expériences. À chaque
appel à la fonction rand, l’état du générateur change et donc la graine avec lui.
On peut choisir la graîne avec la syntaxe rand("seed", n), où n est une variable contenant
un entier.
On peut sauvegarder l’état actuel du générateur dans une variable avec la syntaxe n = ←-
rand("seed").
Le dernier script, utilisant la fonction getdate se sert de l’horloge interne de la machine
pour choisir « au hasard » une graine.
Corrigé de l’exercice 2.1 Il y a plusieurs façons de procéder. La plus simple est bien sûr
un calcul direct.
" #b
x2
Z b
1[a;b] (x) dx 1 a+b
Z
E [U ] = x dx = x = = .
R b−a a b−a b−a 2 a
2
On calcule de la même manière

1 b3 − a3 b2 + ab + a2
Z b
h i dx
E U2 = x2 = = .
a b−a 3 b−a b−a
Enfin,
h i (b − a)2
Var [U ] = E U 2 − (E [U ])2 = .
12
Une remarque au passage : la variance est invariante par translation, ce qui permet de détecter
d’éventuelles erreurs. Une autre possibilité est le calcul de la transformée de Laplace, voir la
remarque 3.
Corrigé de l’exercice 2.2 On utilise simplement que si V suit la loi uniforme sur [0 ; 1],
alors U := a + (b − a) V suit la loi uniforme sur [a ; b].

function U = runif (m , n , a , b )
y = a + ( b - a ) * rand (m , n ) ;
endfunction
Corrigé de l’exercice 2.3 Comme X est bornée, sa fonction génératrice des moments est dé-
finie sur R. Le calcul suivant est classique, il utilise la somme des termes d’une suite géométrique
et une méthode de symétrisation des quotients bien connue.
b a+b+1
h
tX
i 1X 1 eat − et(b+1) 1e 2 t sinh nt
E e = etk = = 2
.
n k=a n 1−e t n e 2t sinh 2t
Comme la transformée de Laplace est finie sur R, elle est analytique (et même entière) sur R et
admet pour développement en série entière
h i ∞
X h i tk
E etX = 1 + E Xk .
k=1
k!
h i
Ainsi, pour connaître les moments (c’est-à-dire les valeurs E X k ) de X, il suffit de faire un
développement limité de MX (t) au voisinage de 0.
On pose α = a+b
2 . Les formules classiques sur les développements limités donnent :
! ! −1
1 α2 2 n n3 3 1 1 3
MX (t) = 1 + αt + t + o t2 t+ t + o t3 t+ t + o t3
n 2 2 8×6 2 8×6
! !
α2 2 n2 1

= 1 + αt + t + o t2 1 + t2 + o t2 1 − t2 + o t2
2 24 24
(a + b)2 n2 − 1
!
t2
= 1 + αt + + + o t2 .
2 2 12
En identifiant les coefficients, on trouve bien
2
a+b a+b n2 − 1
h i
2
E [X] = et E X = + .
2 2 12
Donc la variance est bien celle donnée dans l’exercice.
On peut retrouver la formule donnant la somme des carrés des n premiers entiers, en prenant
a = 1, b = n et en remarquant que dans ce cas,
n 2 !
n2 − 1 n 2n2 + 3n + 1

n+1 n (n + 1) (2n + 1)
X h i
2 2
k =n×E X =n + = = .
k=1
2 12 6 6
On pourrait, bien sûr procéder de la même manière pour calculer la somme des cubes, il
faudrait alors faire un développement limité à l’ordre 3 de MX (t) (ce qui n’est pas beaucoup
plus long car il suffit de rajouter le terme d’ordre 3 dans l’exponentielle).
Remarque 34. Il y a bien sûr d’autres façons de faire. Rien qu’en restant avec ce point de vue
probabiliste, on aurait ici pu calculer la fonction génératrice des probabilités (car X est à valeurs
entières) et trouver son moment d’ordre 2 en utilisant sa dérivée seconde en 1.
Sans utiliser de probabilité, il existe beaucoup de méthodes élémentaires. Mais en général le
problème de la somme des puissances des entiers consécutifs est un problème qui n’est pas si
simple. Vous pouvez, par exemple, consulter la page Wikipedia sur les nombres de Bernoulli si
cela vous intéresse.
Corrigé de l’exercice 2.4 Il suffit de remarquer que, si U a la loi uniforme sur [0 ; 1], alors
V := a + b(b − a + 1) U c a la loi uniforme sur Ja ; bK. Une possibilité de mise en œuvre en Scilab
est la suivante.

function U = runif_int (m ,n ,a , b )
U = a + floor (( b - a +1) * rand (m , n ) ) ;
endfunction
Corrigé de l’exercice 3.1

1. Il suffit de remarquer que si X et Y sont indépendantes, de carré intégrable, alors XY est
intégrable (par l’inégalité de Cauchy-Schwarz) et E [XY ] = E [X] E [Y ], donc
h i
Var [X + Y ] = E (X + Y )2 − (E [X] + E [Y ])2
h i h i h i h i
= E X 2 + E Y 2 + 2E [X] E [Y ] − E X 2 − E Y 2 − 2E [X] E [Y ] .
= Var [X] + Var [Y ] .
Pn
2. Soient Y1 , . . . , Yn indépendantes de loi de Bernoulli de paramètre p et X = i=1 Yi . Alors X
a bien pour loi B (n, p). Par définition, on a l’égalité entre événements
G
{X = k} = {∀i ∈ A, Yi = 1 et ∀j ∈
/ A, Yj = 0} .
A⊂J1;nK, Card(A)=k
Chacun des événements de l’union disjointe a pour probabilité pk (1 − p)n−k . Il ne reste qu’à
compter le nombre de parties à k éléments de l’ensemble à n éléments J1 ; nK, ce qui, par
n
définition est k . Pour l’espérance de X, il suffit d’utiliser la linéarité :
n
X
E [X] = E [Yk ] = np.
k=1
Pour la variance, par indépendance,
n
X
Var [X] = Var [Yk ] = np (1 − p) .
k=1

1. Il suffit de remarquer que, si U suit la loi uniforme sur [0 ; 1], alors la variable aléatoire 1{U ≤p}
suit la loi de Bernoulli de paramètre p.
function Y = rbernoulli (m ,n , p )
Y = double ( rand (m , n ) < p ) ;
endfunction
Ici, on a utilisé la fonction double pour convertir les booléens en flottants à double précision.
Ce n’est pas obligatoire et même superflu dans la plupart des cas.
2. On utilise la définition. Si l’on ne cherche pas à vectoriser, voici une première solution avec
boucles, qui n’est pas très efficace en Scilab mais a le mérite d’être assez simple.
function Y = r b i n o m i a l e _ a v e c _ b o u c l e s (m ,n ,r , p )
Y = zeros (m , n ) ; // pr é allocation
// parcours de Y
for j = [1: n ]
for i = [1: m ]
Y (i , j ) = sum ( rand (1 , r ) < p ) ;
end
end
endfunction

Si l’on cherche à vectoriser, il faut en fait créer un tableau de variables de Bernoulli en 3 di-
mensions, sommer suivant la première dimension, et effacer la dimension inutile (la première)
en utilisant la fonction matrix qui sert à redimensionner une matrice.
function Y = rbinomiale (m ,n ,r , p )
bern = ( rand (r ,m , n ) < p ) ;
Y = matrix ( sum ( bern , " r " ) , m , n ) ;
endfunction
3. On commence par calculer la fonction de répartition F de la loi binomiale, avec le calcul de

sa loi d’abord 5 dans la fonction binomiale_probas puis un appel à cumsum. On utilise la
méthode de l’inverse de la fonction de répartition. Étant donné une variable aléatoire uniforme
sur [0 ; 1], on cherche (avec la fonction gsearch) le plus petit entier k tel que
U ≤ F (k) .
Cet entier suit la loi binomiale de paramètres r et p. En raisons de problèmes numériques,

ce programme ne fonctionne pas pour des valeurs trop grandes de r. Dans ce cas, utiliser la
méthode directe ci-dessus ou bien l’approximation par la loi normale.
function P = binomiale_probas (p , r )
// Calculer la loi B (r , p ) sous la forme d ' un
// tableau à r +1 é l é ments
P = zeros (1 , r +1) ;
P (1) = (1 - p ) ^ r ;
for i = [1: r ]
P ( i +1) = P ( i ) * (r - i +1) / i * p /(1 - p ) ;
end
endfunction
function Y = rbinomiale2 (m ,n ,r , p )
// simulation de la loi binomiale par l ' inverse
// de la fonction de r é partition
P = cumsum ( binomiale_probas (p , r ) ) ;
// /!\ Des valeurs trop elevees de r
// entrainent des valeurs " egales " dans P
// et donc une erreur dans dsearch
Y = dsearch ( rand (m , n ) , P ) ;
endfunction
Corrigé de l’exercice 3.3 On note, pour k ∈ N, pk = P (X = k).

1. Par définition,
n
X
GX (1) = pk = 1,
k=0
ce qui prouve que le rayon de convergence de la série entière GX est au moins égal à 1. Ceci
assure que pour tout s ∈ [0 ; 1[, GX est dérivable en s et
∞
G0X (s) =
X
kpk sk−1 .
k=1
5. On aurait aussi pu appeler directement la fonction binomial de Scilab.

Pour étudier la dérivée à gauche en 1, on écrit

∞ ∞
GX (s) − 1 X sk − 1 X
= pk = pk 1 + s + · · · + sk−1 .
s−1 k=1
s−1 k=1
Comme la fonction
s 7→ 1 + s + · · · + sk−1
est croissante (pour s > 0), on peut utiliser le théorème de convergence monotone avec la
mesure de comptage :
∞
GX (s) − 1 X
lim = lim pk 1 + s + · · · + sk−1
s↑1 s−1 s↑1
k=1
∞
X
= lim pk 1 + s + · · · + sk−1
s↑1
k=1
X∞
= kpk .
k=1
Ce raisonnement reste valable lorsque E [X] = +∞ .

2. On raisonne par récurrence. Soit n ≥ 1. Supposons le résultat vrai au rang n et E X n+1 < ∞.

Comme le rayon de convergence de GX est au moins égal à 1, pour tout s ∈ [0 ; s[,

+∞
(n) X
GX (s) = k (k − 1) · · · (k − n + 1) pk sn−k+1
k=n
et par hypothèse de récurrence, cette formule est encore vraie en s = 1. On raisonne ensuite
comme ci-dessus, par convergence monotone :
(n) (n) ∞
GX (s) − GX (1) X sk−n+1 − 1
lim = lim pk (k(k − 1) · · · (k − n + 1))
s↑1 s−1 s↑1
k=n+1
s−1
∞
X
= lim pk (k(k − 1) · · · (k − n + 1)) 1 + s + · · · + sk−n
s↑1
k=n+1
∞
X
= lim pk (k(k − 1) · · · (k − n + 1)) 1 + s + · · · + sk−n
s↑1
k=1
X∞
= pk (k(k − 1) · · · (k − n + 1)(k − n)) .
k=1
ceci est encore valable si tous les moments E [X] , E X 2 , . . . , E [X n ] sont finis

Au passage,
n+1
mais E X = +∞.
3. Le calcul de la variance se fait très simplement.
h i
Var [X] = E X 2 − (E [X])2 = E [X (X − 1)] + E [X] − E [X]2 .
Puis il suffit de remplacer E [X] par G0X (1) et E [X (X − 1)] par G00X (1).
Corrigé de l’exercice 3.4 Soit p ∈ ]0 ; 1[ et q = 1−p. Soient Y1 , Y2 , . . . des variables aléatoires

indépendantes de loi Bernoulli de paramètre p et
T = inf {k ≥ 1 | Yk = 1} .
Par définition, T suit la loi géométrique de paramètre p.

1. Soit k ≥ 1. On a l’égalité d’événements
{T = k} = {Y1 = 0, Y2 = 0, . . . , Yk−1 = 0, Yk = 1} ,
d’où le résultat, par indépendance. Au passage, il serait en théorie possible que T soit égale à
+∞, si toutes les variables de Bernoulli étaient nulles. Comme ∞ k−1 = 1, ceci n’arrive
P
k=1 pq
presque sûrement jamais.
2. On calcule la fonction génératrice avec la série géométrique. Pour tout s ∈ [0 ; 1],
∞
X ps
GT (s) = pq k−1 sk = .
k=1
1 − qs
3. On en déduit l’espérance et la variance par dérivations successives.

p
G0T (s) = ,
(1 − qs)2
p 1
donc E [T ] = (1−q)2
= p (résultat très intuitif !).
2pq
G00T (s) = ,
(1 − qs)3 )
donc
2q
E [T (T − 1)] = ,
p2
puis
2q 1 1 q
Var [T ] = 2
+ − 2 = 2.
p p p p
1. Soit k ≥ 1. En conditionnant successivement et en utilisant la propriété de Markov, on trouve :
P (T = k) = P (X1 = 0, X2 = 0, . . . , Xk−1 = 0, Xk = 1) = q × q × · · · × q × p = q k−1 p,
donc T a bien la loi géométrique de paramètre p.

2. On travaille avec la chaîne canonique et on note Θ l’opérateur de décalage. On note que sur
l’événement {X1 = 0}, on a T = 1 + T ◦ Θ, puis on écrit
h i h i
E T k X1 = E 1 × 1{X1 =1} + T k 1{X1 =0} X1
h i
= 1{X1 =1} + 1{X1 =0} E (1 + T ◦ Θ)k X1 .
D’après la propriété de Markov faible au temps 1,

h i h i
1{X1 =0} E (1 + T ◦ Θ)k X1 = 1{X1 =0} E (1 + T )k .
En prenant l’espérance, on trouve

h i h i
E T k = p + qE (1 + T )k ,
h i
d’où le résultat, en isolant E T k et en divisant par p.

3. Pour k = 1, on trouve directement

q 1
E [T ] = 1 + = .
p p
Pour k = 2, on trouve
h i q p2 + qp + 2q p2 + (1 − p) p + 2q p + 2q
E T2 = 1 + (1 + 2E [T ]) = 2
= 2
= .
p p p p2
La variance vaut donc h i p + 2q 1 q
E T 2 − (E [T ])2 = 2
− 2 = 2.
p p p
1. Une mise en œuvre possible du premier algorithme est la suivante.
function Y = rgeom (m ,n , p )
// Pre - allocation de Y
Y = zeros (m , n ) ;
// Parcours de Y
for j = [1: n ]
for i = [1: m ]
// simulation avec boucle tant que
k = 0; // nombre d ' essais
succes = %f ; // false
while (~ succes ) // ~ pour n é gation
k = k + 1;
if ( rand () < p ) // succes
Y (i , j ) = k ;
succes = %t ; // true
end
end
end
end
endfunction
2. On a l’égalité des événements :
ln (U ) ln (U )
n o
b c+1=k = k−1≤ < k = q k−1 ≥ U > q k
ln (q) ln (q)
Ainsi,
ln (U )

P b c + 1 = k = q k−1 − q k = q k−1 p,
ln (q)
donc b ln(U )
ln(q) c + 1 suit bien la loi géométrique de paramètre p.
3.
function Y = rgeom2 (m ,n , p )
Y = floor ( log ( rand (m , n ) ) / log (1 - p ) ) + 1;
endfunction


1. Soit s ∈ [0 ; 1].
∞
λk k
e−λ
X
GN (s) = s
k=0
k!
∞
(λs)k
= e−λ = e−λ eλs = eλ(s−1) .
X
k=0
k!
2. On dérive deux fois la fonction génératrice :

G0N (s) = λeλ(s−1) ,
G00N (s) = λ2 eλ(s−1) .
On en déduit que
E [N ] = G0N (1) = λ,
et
E [N (N − 1)] = G00N (1) = λ2 .
Finalement,
Var [N ] = λ2 + λ − λ2 = λ.
3. Soit k ∈ N fixé. Soit n assez grand pour avoir nλ < 1 et n ≥ k. Comme les lois considérées
sont à valeurs discrètes 6 , il suffit, pour étudier la convergence en loi de montrer que pour
tout k ∈ N,
P [Xn = k] −−−→ P [N = k] .
n→∞
Il suffit de calculer :
!
k n−k n−k
n λ λ n (n − 1) (n − 2) · · · (n − k + 1) λk λ

P [Xn = k] = 1− = 1− .
k n n k! nk n
On a les équivalents suivants :
n (n − 1) (n − 2) · · · (n − k + 1) ∼n→∞ nk ,
n−k
λ

1− −−−→ e−λ .
n n→∞
Finalement,
λk
P [Xn = k] −−−→ e−λ = P [N = k] .
n→∞ k!
1. Soit f une fonction mesurable positive. Par définition,
n n n
" !# Z Z !
Pn
xi λn e−λ xi
X X Y
E f Xi = ··· f i=1 1{xi ≥0} dxi .
i=1 i=1 i=1
On fait le changement de variable linéaire bijectif, de jacobien 1 :

x1 = x01
x2 = x02 − x01
..
.
xn = x0n − x0n−1 .
6. ATTENTION ! C’est totalement faux dans les autres cas !

Ainsi, l’intégrale devient

n
" !# Z Z
0
f x0n λn e−λxn 1{0≤x0 ≤x0 ≤···≤x0 } dx01 . . . dx0n .
X
E f Xi = ···
1 2 n
i=1
Une récurrence permet de calculer :
x0n n−1
Z Z
0 0 0
··· 1{0≤x0 ≤···≤x0 0
} dx1 dx2 · · · dxn−1 = ,
1 n−1 ≤xn (n − 1)!
finalement (en posant x = x0n ),

n
" !#
xn−1 n −λx
X Z
E f Xi = f (x) λ e 1{x≤0} dx
i=1
(n − 1)!
Pn
donc la loi de i=1 Xi est la loi Γ (n, λ).
2. Comme, pour tout i ≥ 1, on a Ui ≤ 1,
( n )
n o X − log (Ui )
−λ
{N ≥ n} = U1 · · · Un ≥ e = ≤1 .
i=1
λ
En posant Xi = − log(U i)
λ , on obtient des variables indépendantes de loi exponentielle de
paramètre λ. Finalement, P (N ≥ n) = P (Z ≤ 1), où Z suit la loi Γ (n, λ).
3. Si n ≥ 2, on obtient en intégrant par parties,
xn−1 n −λx
Z 1
P (N ≥ n) = λ e dx
0 (n − 1)!
λn−1 −λ xn−2 n−1 −λx
Z 1
=− e + λ e dx.
(n − 1)! 0 (n − 2)!
On en déduit que
λn−1 −λ
P (N ≥ n) − P (N ≥ n − 1) = − e ,
(n − 1)!
puis que
λn−1 −λ
P (N = n − 1) = e ,
(n − 1)!
donc N suit la loi de Poisson de paramètre λ.
4. Une mise en œuvre possible de cet algorithme est la suivante :
function Y = rpoisson (m , n , lambda )

barriere = exp ( - lambda ) ;
for i = [1: m ]
for j = [1: n ]
// simulation
nb_unifs = 0;
produit = 1.;
while %t
u = rand () ;
produit = produit * u ;
if ( produit <= barriere )

break ;
end
nb_unifs = nb_unifs + 1;
end
// nb_unifs suit la loi poisson ( lambda )
// fin de la simulation
Y (i , j ) = nb_unifs ;
end
end
endfunction
Corrigé de l’exercice 4.1 Éléments de correction.

1. On montre assez facilement, en utilisant la croissance de f que, pour tout x dans l’intérieur
de I,
lim f (y) = sup {f (y) | y < x} .
y→x
y<x
De même pour la limite à droite. L’inégalité est évidente.

Au passage, on peut raisonner de la même manière lorsque x est une borne (incluse ou non)
de l’intervalle I, mais les limites infinies peuvent survenir.
2. D’après ce qui précède, x est un point de discontinuité si et seulement si f (x−) < f (x+).
Dans ce cas, l’intervalle ]f (x−) ; f (x+)[ contient des rationnels. Choisissons-en un, disons
qx . Par croissance de f , le rationnel qx n’appartient à aucun autre intervalle du même type.
Ainsi les points de discontinuité de f s’injectent dans l’ensemble des rationnels et sont donc
en quantité dénombrable.
3. Monotonie Si x et x0 sont deux réels tels que x ≤ x0 , alors ]−∞ ; x] ⊂ ]−∞ ; x0 ] donc
PX (]−∞ ; x]) ≤ PX (]−∞ ; x0 ]).
Limites Comme ∞ N =1 ]−∞ ; −N ] = ∅, et que cette intersection est décroissante, on a
T
lim PX (]−∞ ; k]) = PX (∅) = 0.

k→−∞
On procède de même avec la réunion croissante

∞
[
]−∞ ; N ] = R.
N =1
Continuité à droite Comme F est croissante, il suffit de montrer que pour tout réel x,
1

F x+ −−−→ F (x) ,
n n→∞
ce qui découle (l’intersection étant décroissante) de

 
+∞
1 1
\
lim PX ]−∞ ; x + ] = PX  ]−∞ ; x + ] = PX (]−∞ ; x]) .
n→∞ n n≥1
n
Le fait que F admette des limites à gauche en tout point est une conséquence de sa
monotonie.
4. Soit U un ouvert. On écrit G
U= Ui ,
i∈I
où les Ui sont ses composantes connexes (donc disjointes). Étant ouvertes (et fermées) dans
l’ouvert U , elles sont ouvertes dans R et comme ce sont des connexes de R, ce sont des
intervalles ouverts. Chacune d’entre elles possède un rationnel. Comme elles sont disjointes
deux à deux, on peut injecter l’ensemble d’indices I dans l’ensemble dénombrable Q.

5. P (X ∈ ]a ; b[) = P (X ∈ ]−∞ ; b[) − P (X ∈ ]−∞ ; a])

∞
!
[ 1
= PX ]−∞ ; b − ] − F (a) = F (b−) − F (a) .
n=1
n
6. La fonction F caractérise la loi de X sur les intervalles ouverts d’après le point précédent,
donc sur les unions dénombrables disjointes de tels intervalles, c’est-à-dire tous les ouverts.
Comme l’ensemble des ouverts est une famille de parties stable par intersections finies qui
engendre la tribu borélienne, par le lemme de classe monotone, F caractérise la loi de X sur
la tribu borélienne.

1. Si 0 < u < u0 < 1, alors, par croissance de F ,
x ∈ R F (x) ≥ u0 ⊂ {x ∈ R | F (x) ≥ u} ,

donc
inf x ∈ R F (x) ≥ u0 ≥ inf {x ∈ R | F (x) ≥ u} ,

donc F −1 est croissante sur ]0 ; 1[.

2. Soit u ∈ ]0 ; 1[. Alors, pour tout entier n ≥ 1,
1
F −1 (u) + > F −1 (u) ,
n
donc par définition de F −1 ,
1

−1
F F (u) + ≥u
n
Comme F est continue à droite, on en déduit
1

−1 −1
F F (u) = lim F F (u) + ≥ u.
n→∞ n
3. Comme F −1 est croissante,
il suffit de montrer que F −1 (u−) = F (u) pour tout u ∈ ]0 ; 1[.
On sait que F −1 u − n1 tend en croissant vers F −1 (u−). Comme F est continue à droite,
on a alors
1

lim F F −1 u − = F F −1 (u−) .
n→∞ n
Comme pour tout n ≥ 1,
1

−1 −1
F (u−) ≥ F u− ,
n
On obtient, en appliquant la fonction croissante F , en utilisant l’inégalité de la question
précédente puis en faisant tendre n vers +∞,

F F −1 (u−) ≥ u.
Donc par définition, F −1 (u−) ≥ F −1 (u) et par conséquent, F −1 (u−) = F −1 (u).

4. L’implication directe découle directement de la définition de F −1 . Puis, si F −1 (u) ≤ x, en
appliquant la fonction croissante F et en utilisant l’inégalité de la question 3,

F (x) ≥ F F −1 (u) ≥ u.
5. Si U suit la loi uniforme sur [0 ; 1], alors, d’après la question précédente,

P F −1 (U ) ≤ x = P (U ≤ F (x)) = F (x) .
Ainsi, les variables aléatoires X et F −1 (U ) ayant la même fonction de répartition, ont la

même loi.


1. Pour la loi de Poisson (on évite de calculer les puissances et les factorielles à chaque étape
pour des raisons de performance), on peut procéder comme suit.
function Y = rpoisson2 (m ,n , lambda )

for i = [1: m ]
for j = [1: n ]
U = rand () ;
s = 0;
k = 0;
p = exp ( - lambda ) ;
while %t
if U < p
break ;
end
k = k + 1;
p = p * lambda / k ;
end
Y (i , j ) = k ;
end
end
endfunction
2. En utilisant les fonctions cumsum et dsearch, on peut vectoriser le code et obtenir un pro-
gramme très efficace.
function Y = rfinite (m , n , P )
C = cumsum ( P ) ;
U = rand (m , n ) ;
Y = dsearch (U , C ) ;
endfunction

1. Soit F la fonction de répartition de la loi exponentielle de paramètre λ > 0. Des calculs
élémentaires montrent que la fonction de répartition de la loi exponentielle vaut
F (x) = 1 − e−λx .
En résolvant l’équation d’inconnue u, F (x) = u, on trouve que la fonction F est bijective sur
]0 ; 1[, de « vraie » inverse
ln (1 − u)
F −1 (u) = − .
λ
2. On sait donc que si U est uniforme sur [0 ; 1], la variable aléatoire − ln(1−U
λ
)
suit la loi expo-
nentielle de paramètre λ. Bien sûr, comme 1 − U suit également la loi uniforme sur [0 ; 1], il
en est de même de − ln(U
λ
)
. Le programme fait une ligne.
function Y = rexp (m , n , lambda )

Y = - log ( rand (m , n ) ) / lambda ;
endfunction

3. Soit X suivant la loi de Cauchy centrée réduite et F sa fonction de répartition. En calculant

l’intégrale
1 x 1 1 π
Z
dt = arctan (x) − ,
π −∞ 1 + t2 π 2
on trouve que la fonction de répartition de X est F (x) = π1 arctan (x) + 12 . Elle est bijective
sur R et son inverse est
π

−1
F (u) = tan πu − .
2
4. D’après le critère de Riemann, l’intégrale
Z ∞
1
|t| dt
−∞ 1 + t2
diverge, donc la loi de Cauchy n’admet pas d’espérance.

1. Il suffit de tirer des variables uniformes sur C sous la forme d’un couple (U, V ) avec U et V
indépendantes de lois uniforme sur [0 ; 1] et d’accepter le tirage lorsque U 2 + V 2 ≤ 1. D’après
la partie précédente, la loi du temps d’attente d’un tirage accepté est la loi géométrique de
paramètre π/4, donc d’espérance 4/π.
2. Il suffit d’appliquer la loi des grands nombres avec les variables aléatoires indépendantes de
Bernoulli 1{Uk ∈Q} d’espérances π/4.
3. Voici un exemple de script possible.
// Calcul de pi par la m é thode de Monte - Carlo

num_tirages = 10000;
num_succes = 0;
for i = [1: num_tirages ]
U = rand (1 , 2) ;
if ( sum ( U .^ 2) <= 1 )
num_succes = num_succes + 1;
end
end
approx_pi = 4 * num_succes / num_tirages ;
printf ( " Approximation de pi pour %d tirages : %f " , num_tirages ←-
, approx_pi ) ;
En voici un autre, sans boucle (donc beaucoup plus rapide) :
// Calcul de pi par la m é thode de Monte - Carlo

// Version vectoris é e
num_tirages = 10000;
U = rand (1 , num_tirages ) ;
V = rand (1 , num_tirages ) ;
Z = U .^2 + V .^2;
num_succes = sum ( Z <= 1) ;
approx_pi = 4 * num_succes / num_tirages ;
printf ( " Approximation de pi pour %d tirages : %f " , num_tirages ←-
, approx_pi ) ;

4. Soit Y1 , Y2 , . . . , des variables i.i.d. de même loi de Bernoulli de paramètre π4 et Sn = nk=1 Yk .

P
Pour alléger les notations et gagner en généralité, on pose p = π4 , q = 1 − p, ε = 0, 01 et

α = 0, 05. On fait donc l’hypothèse que n est assez grand pour assimiler Snn à une loi normale
pq
N p, n . On cherche n ≥ 1 tel que
Sn

P − p ≤ ε ≥ 1 − α.
n
q
pq
On pose Snn = p + n N, où N suit la loi normale centrée réduite. L’inéquation ci-dessus est
équivalente à
n
r
P |N | ≥ ε ≥ 1 − α.
pq
Ce qui, pour des raisons de symétrie, équivaut encore à
n α
r
P N ≥ε ≥1− .
pq 2
On applique la fonction quantile FX−1 pour obtenir finalement

2
pq FN−1 1 − α
2
n≥ .
ε2
Ce qui est important à noter est la proportionnalité à ε−2 . L’application numérique donne
ici environ n ≈ 4600 avec le script suivant.
p = %pi / 4;
q = 1 - p;
eps = 0.01;
alpha = 0.05;
Fn_m1 = cdfnor ( " X " , 0 , 1 , 1 - alpha , alpha ) ;
n_min = p * q * Fn_m1 ^2 / eps ^2;
disp ( n_min ) ;

1. Il est évident que f est positive et d’intégrale 1. Sa primitive nulle en −∞, c’est-à-dire sa
fonction de répartition est F (X) = x3 et l’inverse de sa fonction de répartition est primitive
nulle en −∞, c’est-à-dire sa fonction de répartition est F (X) = x3 et l’inverse de sa fonction
√
de répartition est F −1 (X) = 3 x. Le programme suivant simule la loi de densité f .
function Y = inv_rep (m , n )
Y = rand (m , n ) .^ (1. / 3.) ;
endfunction
2. On doit trouver une fonction de densité g d’une loi connue et facilement simulable, et une
constante c la moins grande possible telle que pour tout x ∈ [0 ; 1], f (x) ≤ cg(x). On va
utiliser ici la loi uniforme sur [0 ; 1] (donc ici g est constante et égale à 1) et c = 3.
Pour simuler X, on commence par tirer V suivant la loi uniforme sur [0 ; 1] avec rand, puis
U une autre variable uniforme sur [0 ; 1] indépendante de V . On accepte V comme un tirage
de X si
f (V ) 3V 2
U≤ = = V 2.
cg (V ) 3×1

function Y = rejet (m , n )
for i = 1: m
for j = 1: n
while %t
essai = rand () ;
valeur_test = rand () ;
if valeur_test <= essai ^2
break ; // on accepte
end
end
Y (i , j ) = essai ;
end
end
endfunction
3. Le nombre d’itération dans la boucle while est aléatoire, de loi géométrique de paramètre 1c
donc son espérance est c = 3. Si l’on demande un nombre N important de valeurs, d’après la
loi des grands nombres, le nombre moyen d’itérations dans la boucle sera proche de 3. Comme
chaque boucle while fait appel deux fois à la fonction rand, cela fait en moyenne 6 appels à
rand par valeur renvoyée.
C = [3 1 ; 1 2];
A = chol ( C ) ;
mu = [ -2 ; -3];
n = 1000;
N = rand (2 , n , " normal " ) ;
X = A*N;
X (1 ,:) = X (1 ,:) + mu (1) ;
X (2 ,:) = X (2 ,:) + mu (2) ;
clf () ;
plot ( X (1 ,:) , X (2 ,:) , " r . " )
[Q , D ] = spec ( C ) ;
u = Q (: ,1) * sqrt ( D (1 ,1) ) ;
xarrows ([ mu (1) , u (1) ] , [ mu (2) , u (2) ]) ;
v = Q (: ,2) * sqrt ( D (2 ,2) ) ;
xarrows ([ mu (1) , v (1) ] , [ mu (2) , v (2) ]) ;
// marche al é atoire en dimension 2

t_max = 100000;
// d é placements : haut bas gauche droite equiprobables
D = [1 0 ; 0 1 ; -1 0 ; 0 -1];
// on choisit des nombres entre 1 et 4
// pour connaitre les d é placements
Alea = grand (1 , t_max , " uin " , 1 , 4) ;
// deplacements aleatoires

Z = zeros (2 , t_max ) ;
for i = [1:4]
Z (1 , Alea == i ) = D (i ,1) ;
Z (2 , Alea == i ) = D (i ,2) ;
end
Positions = cumsum (Z , " c " ) ;
plot ( Positions (1 ,:) , Positions (2 ,:) ) ;
Corrigé de l’exercice 7.6 Voici un exemple de ce que l’on pouvait faire. J’ai rangé les deux
fonctions principales dans un fichier de fonctions galtonwatson.sci que voici.
// galtonwatson . sci
// fichier de fonctions pour des experiences
// sur les processus de galton - watson
function Y = support_fini (m ,n , eta )

// simule la loi definie par le vecteur
// stochastique eta sur l ' ensemble fini
// 0 , 1 , ... , length ( eta ) - 1
Y = dsearch ( rand (m , n ) , cumsum ( eta ) ) ;
endfunction
function Z = galton_watson (n , xi )
// renvoie une simulation de
// ( Z_0 , Z_1 , ... , Z_n )
// xi ( N ) doit ê tre un vecteur al é atoire
// de taille N compose d ' entiers naturels
Z = zeros (1 , n +1) ; Z (1) = 1; // Z (1) = Z_0
for j = 1: n
Z ( j +1) = sum ( xi ( Z ( j ) ) ) ;
end
endfunction
Pour la suite, on peut séparer en plusieurs scripts. En voici pour le tracé d’une dizaine de courbes,
le fichier de script s’appelle gw_trace.sce. Au début du script, on appelle le fichier de fonction
pour en avoir toujours la dernière version (et être sûr qu’il soit chargé).
// gw_trace . sce
// Script pour experimenter sur les
// processus de Galton - Watson
exec ( ' galtonwatson . sci ') ;
// choix de la loi
eta = [.2 .4 .4];
function Y = xi ( n )
Y = support_fini (1 ,n , eta ) ;
endfunction
// un premier echantillon
n = 15
disp ( galton_watson (n , xi ) )

// Plusieurs experiences
num_exp = 20;
Z = zeros ( num_exp , n +1) ;
for i = 1: num_exp
Z (i ,:) = galton_watson (n , xi ) ;
end
// Trace des num_exp trajectoires
scf (1) ; clf () ;
plot (Z ') ; // trace chaque colonne de Z '
La figure 1 montre un exemple de ce que l’on pouvait obtenir. Pour la suite de l’exercice, on
Figure 1 – 20 trajectoires de (Z0 , Z1 , . . . , Z15 ) dans le cas où η0 = 1/5 et η1 = η2 = 2/5.
pouvait poursuivre avec le script suivant.
// gw_extinction . sce
// extimation de la probabilite d ' extinction
// au bout de n generations

n = 15
// choix de la loi
eta = [.2 .4 .4];
Y = support_fini (1 ,n , eta ) ;
endfunction
// estimation de P ( Z_ {15} = 0)
num_exp = 2000;
cpt_vivant = 0;
GW = zeros (1 , n +1) ;
for i = 1: num_exp
GW = galton_watson (n , xi ) ;
if ( GW (16) ~= 0)
cpt_vivant = cpt_vivant + 1;
end
end
disp ( " fr é quence empirique pour " ) ; disp ( num_exp ) ; disp ( " lancers←-
: ");
disp ( cpt_vivant / num_exp ) ;
// Comparaison avec la valeur theorique
function y = f ( x )
y = (1. + 2. * x + 2. * x .^2) / 5.;
endfunction
qn = 0;
for i = 1: n
qn = f ( qn ) ;
end
disp ( " valeur theorique : " ) ; disp (1 - qn ) ;
On devrait trouver des valeurs « assez proches » l’une de l’autre mais pour avoir une idée plus
théorique de l’ordre de grandeur de l’écart, il faudrait davantages de renseignements sur la loi
de Z15 (par exemple une majoration de son écart-type).
Il n’est pas très difficile de changer de loi de reproduction, il suffit simplement de changer la
fonction xi. Voici, par exemple comment on pouvait choisir η comme étant la loi de Poisson de
paramètre 3 en utilisant la fonction prédéfinie grand. On a dû réduire le nombre de générations
à 12 pour éviter un dépassement de la pile.
// gw_trace_poisson . sce
// loi de reproduction
Y = grand (1 , n , " poi " , 3)
endfunction
// Nombre de generations
n = 12;
num_exp = 20;
Z = zeros ( num_exp , n +1) ;
for i = 1: num_exp
Z (i ,:) = galton_watson (n , xi ) ;
end
// Trace des num_exp trajectoires

scf (1) ; clf () ;

plot (Z ') ; // trace chaque colonne de Z '
On peut observer dans la figure 2 que les ordres de grandeurs ne sont pas les mêmes que dans
l’exemple précédent.
Figure 2 – 20 trajectoires de (Z0 , Z1 , . . . , Z12 ) dans le cas où la loi de reproduction est la loi de
Poisson de paramètre 3.
Références
[BC07] Bernard Bercu and Djalil Chafaï. Modélisation stochastique et simulation. Dunod,
2007.
[BP12] Marc Briane and Gilles Pagès. Analyse, Théorie de l’intégration, Convolution et Trans-
formation de Fourier, 5ème éd. Vuibert, Paris, 365p, 2012.
[Can13] Bernard Candelpergher. Théorie des probabilités. Calvage & Mounet, 2013.
[Dur10] Rick Durrett. Probability : theory and examples. Cambridge university press, 2010.

[LG06] Jean-François Le Gall. Intégration, probabilités et processus aléatoires. http://www.

math.u-psud.fr/~jflegall/IPPA2.pdf, 2006. Polycopié de cours.
[LG12] Jean-François Le Gall. Mouvement brownien, martingales et calcul stochastique, vo-
lume 71. Springer Science & Business Media, 2012.
[MT+ 00] George Marsaglia, Wai Wan Tsang, et al. The ziggurat method for generating random
variables. Journal of statistical software, 5(8) :1–7, 2000.
[RS12] Vincent Rivoirard and Gilles Stoltz. Statistique mathématique en action. Vuibert,
2012.
[Rud87] Walter Rudin. Real and complex analysis. McGraw-Hill Higher Education, 1987.
[Tou99] Paul S. Toulouse. Thèmes de probabilités et statistiques. Dunod, Paris, 1999.

Simulation Agreg 1718

Transféré par

Droits d'auteur :

Formats disponibles

Simulation Agreg 1718

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Simulation Agreg 1718

Transféré par

Droits d'auteur :

Formats disponibles

Modélisation option A Simulation

Simulation de variables aléatoires

La dernière version de ce document est disponible sur le site :

Exercice 1.1. rand, rand("seed", s), rand(m, n), histplot

Université Paris 13 1/42 Préparation à l’agrégation 2017 - 2018

7. Que fait la commande histplot(linspace(0, 1, 10), rand(1, 10000)) ? Consulter l’aide

x = getdate ( " s " ) ;

Exemple 2. 1. Prenons E = R, E = B (R) et pour µ, la mesure de Lebesgue. La loi uniforme sur

3. Prenons E = Z, muni de la tribu discrète 2Z et pour µ, la mesure de comptage (c’est-à-dire,

2. En faisant un développement limité de MX à l’ordre 2 en 0, calculer les moments d’ordre 1

Université Paris 13 2/42 Préparation à l’agrégation 2017 - 2018

2. la loi binomiale de paramètres n et p (notation : X ∼ B (n, p)) si X a la même loi que

Université Paris 13 3/42 Préparation à l’agrégation 2017 - 2018

Exercice 3.2. Simulation de la loi binomiale

Exercice 3.3. Fonction génératrice des probabilités

Exercice 3.4. Loi géométrique

2. Montrer que la fonction génératrice des probabilités de T vaut

3. Démontrer que l’espérance et la variance de T sont données par

Exercice 3.5. Loi géométrique avec chaîne de Markov

Soit T le temps d’atteinte de l’état 1. On note q = 1 − p.

Université Paris 13 4/42 Préparation à l’agrégation 2017 - 2018

1. Montrer que T suit la loi géométrique de paramètre p.

avec ici max ∅ = 0. Justifier que P (N ≥ n) = P (X1 + X2 + · · · + Xn ≤ 1).

Université Paris 13 5/42 Préparation à l’agrégation 2017 - 2018

4 Inversion de la fonction de répartition

4.1 Fonction de répartition

// Exemple d ' utilisation de cdfnor

Exercice 4.1. Rappels sur les fonctions croissantes et la fonction de répartition

f (x−) ≤ f (x) ≤ f (x+) .

2. En déduire que l’ensemble des points de discontinuité de f est au plus dénombrable.

4.2 Inverse de Lévy

FX−1 (u) = inf {x ∈ R | FX (x) ≥ u} .

Cette fonction s’appelle l’inverse généralisée ou l’inverse de Lévy de FX ou encore fonction

Université Paris 13 6/42 Préparation à l’agrégation 2017 - 2018

// Inversement , je veux connaitre x

Consulter l’aide et s’entraîner avant l’épreuve semble indispensable.

3. Montrer que F −1 est continue à gauche.

5. Calculer P F −1 (U ) ≤ x puis conclure.

4.3 Loi prenant un nombre dénombrable de valeurs

Université Paris 13 7/42 Préparation à l’agrégation 2017 - 2018

4.4 Application à la loi de Cauchy et la loi exponentielle

5.1 Loi uniforme

Université Paris 13 8/42 Préparation à l’agrégation 2017 - 2018

2. Soient U1 , U2 , . . . i.i.d. de lois uniformes sur C. Justifier que

5.2 Application aux lois à densités

Université Paris 13 9/42 Préparation à l’agrégation 2017 - 2018

Exercice 5.2. On considère la fonction f définie par :

6 Loi normale, vecteurs gaussiens

6.1 Variables gaussiennes

Faire un changement de variable en coordonnées polaires, puis calculer I.

Proposition 12. Soit N ∼ N (0, 1), m ∈ R, σ ≥ 0 et X ∼ N m, σ 2 .

1. L’espérance et la variance de N valent respectivement 0 et 1.

4. Les fonctions caractéristiques de N et de X valent respectivement

Université Paris 13 10/42 Préparation à l’agrégation 2017 - 2018

et faire un changement de variable dans l’intégrale.

Exercice 6.3. Démontrer ce corollaire en utilisant la fonction caractéristique.

Université Paris 13 11/42 Préparation à l’agrégation 2017 - 2018

6.2 Méthode de Box-Muller

en faisant un premier changement de variable y = r2 , puis un passage en coordonnées polaires.