TP de Simulation Aléatoire I
Université Nationale des Sciences, Technologies, Ingénierie et Mathématiques
(UNSTIM)
École Nationale Supérieure de Génie Mathématique et Modélisation
(ENSGMM)
Filière : GMM 1 Année académique : 2024–2025
TCHINHOUN Thibaut - ANATO Kristland
I - Algorithme de Box-Muller
On utilise l’algorithme de Box-Muller pour générer 50 réalisations indépendantes suivant
une loi normale N (35, 5). Les formules sont :
q q
Z1 = −2 ln U1 cos(2πU2 ), Z2 = −2 ln U1 sin(2πU2 ),
où U1 , U2 ∼ U(0, 1), et X = 35 + 5Z.
Listing 1 – Code R pour Box-Muller
1 # Param tres
2 n <- 50 # Nombre de r a l i s a t i o n s
3 mu <- 35 # Moyenne
4 sigma <- 5 # cart - type
5
6 # Algorithme de Box - Muller
7 U1 <- runif ( n / 2)
8 U2 <- runif ( n / 2)
9 Z1 <- sqrt ( -2 * log ( U1 ) ) * cos (2 * pi * U2 )
10 Z2 <- sqrt ( -2 * log ( U1 ) ) * sin (2 * pi * U2 )
11 Z <- c ( Z1 , Z2 )
12 X <- mu + sigma * Z
13
14 # Histogramme
15 hist (X , prob = TRUE , main = " Histogramme ␣ de ␣ N (35 , ␣ 5) " , xlab = " Valeurs
" , col = " lightblue " )
16 curve ( dnorm (x , mean = mu , sd = sigma ) , add = TRUE , col = " red " , lwd =
2)
1
Figure 1 – Histogramme des réalisations de N (35, 5) avec densité superposée.
II - Loi de Weibull par inversion
La fonction de répartition de la loi de Weibull est :
F (x) = 1 − exp(−(λx)a ), x > 0.
Par inversion, pour U ∼ U(0, 1), on a :
(− ln U )1/a
X= .
λ
Listing 2 – Code R pour la loi de Weibull
1 # Param tres
2 n <- 1000
3 a <- 2
4 lambda <- 1
5
6 # Simulation par inversion
7 U <- runif ( n )
8 X <- (( - log ( U ) ) ^(1 / a ) ) / lambda
9
10 # Histogramme
2
11 hist (X , prob = TRUE , main = " Histogramme ␣ de ␣ la ␣ loi ␣ de ␣ Weibull " , xlab =
" Valeurs " , col = " lightblue " )
12 curve ( dweibull (x , shape = a , scale = 1 / lambda ) , add = TRUE , col = " red
" , lwd = 2)
Figure 2 – Histogramme des réalisations de la loi de Weibull.
III - Loi de Cauchy
1 - Fonction de répartition
La densité de la loi de Cauchy est :
1
g(x) = , x ∈ R.
π(1 + x2 )
La fonction de répartition est :
arctan(x) 1
G(x) = + .
π 2
3
2 - Simulation par inversion
Pour U ∼ U(0, 1), on a :
1
X = tan π U − .
2
On utilise u1 = 0.313, u2 = 0.03, u3 = 0.872.
Listing 3 – Code R pour la loi de Cauchy
1 # Valeurs de U
2 U <- c (0.313 , 0.03 , 0.872)
3
4 # Simulation
5 X <- tan ( pi * ( U - 0.5) )
6 print ( X )
Résultats : X1 ≈ −0.655, X2 ≈ −3.078, X3 ≈ 1.253.
IV - Intégrale par Monte-Carlo
On considère l’intégrale :
Z 4
1 y2
I= √ e− 2 dy = Φ(4) − Φ(0).
0 2π
1 - Valeur exacte
Listing 4 – Code R pour la valeur exacte
1 I _ exact <- pnorm (4) - pnorm (0)
2 print ( I _ exact )
Résultat : I ≈ 0.4999683.
2 et 3 - Estimation et IC
Listing 5 – Monte-Carlo et intervalle de confiance
1 n <- 10000
2 a <- 0; b <- 4
3 f <- function ( y ) (1 / sqrt (2 * pi ) ) * exp ( - y ^2 / 2)
4 Y <- runif (n , a , b )
5 f _ Y <- f ( Y )
6 I _ hat <- ( b - a ) * mean ( f _ Y )
7
8 s <- sd ( f _ Y )
9 se <- s / sqrt ( n )
10 IC <- c ( I _ hat - 1.96 * se , I _ hat + 1.96 * se )
11 print ( list ( Estimation = I _ hat , IC _ 95 = IC ) )
4
4 - Convergence
Listing 6 – Convergence de Monte-Carlo
1 n _ max <- 10000
2 I _ n <- numeric ( n _ max )
3 se _ n <- numeric ( n _ max )
4 for ( i in 1: n _ max ) {
5 Y <- runif (i , 0 , 4)
6 f _ Y <- f ( Y )
7 I _ n [ i ] <- 4 * mean ( f _ Y )
8 se _ n [ i ] <- sd ( f _ Y ) / sqrt ( i )
9 }
10 plot (1: n _ max , I _n , type = " l " , ylim = c (0.45 , 0.55) ,
11 xlab = " Taille ␣ de ␣l ’ chantillon " , ylab = " Estimation ␣ de ␣ I " , main
= " Convergence ␣ de ␣ Monte - Carlo " )
12 lines (1: n _ max , I _ n + 1.96 * se _n , col = " blue " , lty = 2)
13 lines (1: n _ max , I _ n - 1.96 * se _n , col = " blue " , lty = 2)
14 abline ( h = I _ exact , col = " red " , lwd = 2)
15 legend ( " topright " , legend = c ( " Estimation " , " IC ␣ 95% " , " Valeur ␣ exacte " )
,
16 col = c ( " black " , " blue " , " red " ) , lty = c (1 , 2 , 1) )
5
Figure 3 – Convergence de l’estimation de I avec intervalle de confiance à 95%.
V - Méthode de rejet
On simule une variable aléatoire avec la densité :
2√
f (x) = 1 − x2 , x ∈ [−1, 1],
π
1
en utilisant la densité instrumentale g(y) = 2
sur [−1, 1].
1 - Constante c
On cherche c tel que f (x) ≤ cg(x). On a :
√
f (x) 2
π
1 − x2 4√
= 1 = 1 − x2 .
g(x) 2
π
Le maximum est atteint en x = 0, soit π4 . Ainsi, c = π4 .
2 - Simulation par rejet
y+1 f (Y )
On génère Y ∼ U(−1, 1) avec G(y) = 2
, soit Y = 2U − 1. On accepte Y si U ′ ≤ cg(Y )
.
6
Listing 7 – Code R pour la méthode de rejet
1 # Param tres
2 n <- 10
3 c <- 4 / pi
4 f <- function ( x ) (2 / pi ) * sqrt (1 - x ^2)
5 g <- function ( x ) 1 / 2
6
7 # Simulation
8 X <- numeric (0)
9 while ( length ( X ) < n ) {
10 U1 <- runif (1) # Pour Y
11 U2 <- runif (1) # Pour le test de rejet
12 Y <- 2 * U1 - 1 # Simulation de Y selon g
13 if ( U2 <= f ( Y ) / ( c * g ( Y ) ) ) {
14 X <- c (X , Y )
15 }
16 }
17
18 # R sultats
19 print ( X )
20
21 # Histogramme pour v r i f i c a t i o n
22 hist (X , prob = TRUE , main = " Histogramme ␣ de ␣ la ␣ loi ␣ cible " , col = "
lightblue " )
23 curve ( f ( x ) , add = TRUE , col = " red " , lwd = 2)
Les réalisations simulées sont :
X ≈ {−0.449, −0.171, −0.796, −0.948, 0.794, −0.343, −0.784, −0.948, −0.893, 0.824}.
Ces valeurs sont tirées de l’exécution du code R avec ‘[Link](123)‘
