Simulation de loi de probabilité sur R
Gelbart Simon
2023-12-13
Exercice 2
Soit \(U\) une variable aléatoire suivant la loi uniforme sur \([0, 1]\). Nous allons montrer que la variable aléatoire \(X = -\frac{1}{\lambda} \ln |1 - U|\) suit une loi exponentielle de paramètre \(\lambda\).
Preuve
La fonction de répartition de \(X\) est donnée par :
\[ F_X(x) = P(X \leq x) = P\left(-\frac{1}{\lambda} \ln |1 - U| \leq x\right) \]
Simplifions cette expression en utilisant le fait que \(\ln |1 - U| \geq 0\) :
\[ F_X(x) = P\left(\ln |1 - U| \geq -\lambda x\right) = P\left(|1 - U| \geq e^{-\lambda x}\right) \]
\[ F_X(x) = P\left(U \leq 1 - e^{-\lambda x}\right) = 1 - e^{-\lambda x} \]
Ceci est la fonction de répartition d’une variable exponentielle de paramètre \(\lambda\), confirmant que \(X\) suit une loi exponentielle. On peut la représenter sur R tel que :
#Question 7
expo <- function(lambda) {
u <- runif(1)
x <- -1/lambda * log(1 - u)
return(x)
} On peut ensuite obtenir sa distribution avec son histogramme à l’aide de la fonction hist pour s nombres:
#Question 8
histoexpo <- function(s, lambda) {
resultat <- replicate(s, expo(lambda))
hist(resultat, main = "Distribution de la loi exponentielle", xlab = "Valeurs", ylab = "Fréquence", col = "lightblue", border = "blue", prob = TRUE,breaks=25) #la fonction hist permet de créer un histogramme
curve(dexp(x, lambda),from=0.01, add = TRUE, col = "black", lwd = 3) #cela permet de représenter sa densité et de la superposer avec notre histogramme
}Lorsqu’on essaye cette fonction comme exemple avec un lambda qui vaut 0,75 et s qui vaut 10 000 et on obtient cette histogramme :
histoexpo(100000, 0.75)
La densité de probabilité de la loi de Cauchy est donnée par :
\[ f(x; x_0, \alpha) = \frac{1}{\pi \alpha \left[1 + \left(\frac{x - x_0}{\alpha}\right)^2\right]} \]
qu’on peut “simplifier” par :
\[ f(x; x_0, \alpha) = \frac{1}{\pi} * \frac{a}{(x-x_0)^2+a^2}\]
elle a la particularité de ressembler à la loi normale,avec une forme de cloche, mais elle est bien plus aplatie. On pourrait utiliser la commande rcauchy pour la représenter mais ici on va utiliser runif et passer par une transformation : \[ X = \tan(\pi \cdot (U - 0.5)) \] Sa fonction sur R peut être représenter telle que :
cauchy <- function(s,n) {
u <- runif(n)
x <- tan(pi * (u - 0.5))
return(x)
}De même on peut illustrer sa distribution avec son histogramme pour s nombres tel que :
histocauchy <- function(t, s) {
set.seed(101) #on peut fixer une seed pour obtenir des résultats similaire à chaque execution
X <- cauchy(scale, s)
# Histogramme pr cauchy
hist(X, prob = TRUE, col = "lightgreen", main = "Histogramme de la loi de Cauchy ", breaks = 10)
# Densité de Cauchy + superposer la courbe à l'histogramme
curve(dcauchy(x, location = 0, scale = 100), add = TRUE, col = "darkgreen", lwd = 2)
legend("topright", legend = c("Simulation", "Densité"), fill = c("lightgreen", "darkgreen"))
}scale_cauchy <- 20 # impacte la largeur
histocauchy(scale_cauchy, 1000)
Exercice 3
Question 11 : \(A_0\) est notre point de départ il est donc toujours égale à 0, en revanche on ne peut pas vraiment parler de loi de probabilité pour \(A_0\) puisque ce n’est pas une variable aléatoire mais une valeur fixe qui vaut 0. On pourrait donc dire que on a une probabilité 1 d’avoir 0.
Question 12 : La loi de probabilité de \(A_N\) dépend de la séquence de pas effectués. Chaque pas peut être 1 ou -1 avec une probabilité égale, la somme des N pas suit une distribution binomiale et donc \(A_N\) suit une distribution binomiale avec paramètres N (nombre de pas) et p=0.5 (probabilité de chaque pas). Donc les lois de probabilité des variables aléatoires A1, A2, A3, A4 et A5 sont des distributions binomiales respectives avec les paramètres N=1,2,3,4 et 5.
Prenons comme exemple \(A_3\) : Il s’agit de la somme de 3 pas donc on a une loi binomiale de paramètres B(3,0.5). On peut représenter \(A_3\) tel que :
n <- 3 # SOMME DE 3 pas puisque c'est A3
p <- 0.5 # Probabilité de succès (un pas positif)
x<-0:n
proba <- dbinom(x, size = n, prob = p)
#histogramme de A3
barplot(proba, names.arg = x, col = "coral4", main = "Distribution binomiale B(3, 0.5)", xlab = "Nombre de pas positifs", ylab = "Probabilité")
On observe que \(A_3\) peut prendre 4 valeurs :
Si aucun pas positif se réalise -3
Si un seul pas positif se réalise -1
Si deux pas positifs se réalisent 1
Si les trois pas sont positifs 3
On peut observer la même chose pour \(A_1, A_2,A_4 et A_5\)
Par exemple, \(A_5\) prend les valeurs {-5, -3, -1, 1, 3, 5} selon le nombre de pas positifs
Question 13 : On a pu voir que \(A_N\) suivait une loi binomiale de paramètres n=N et p=0,5. On peut facilement exprimer la position moyenne du point noir qui est donc l’espérance de \(A_N\).
\[ E[A_N] = N \cdot E[B(1, 0.5)] = N \cdot 0.5 \]
Cependant on a ici une marche aléatoire avec la même probabilité de faire un pas en avant (positif) ou en arrière (négatif) on peut donc écrire:
\[ E[A_N] = 0.5 \times 1 + 0.5 \times (-1) = 0.5 \times 0 = 0 \]
La position moyenne est donc le point 0 :
#N Nombre de pas
N <- 100
n_simulations <- 1000
posf <- rep(0, n_simulations)
for (i in 1:n_simulations) {
posf[i] <- sum(sample(c(-1, 1), N, replace = TRUE))
}
posm <- mean(posf)
# Affichage du graphique
hist(posf, main = "Distribution des Positions Finales",
xlab = "Position Finale", col = "lightblue", xlim = c(-N, N))
abline(v = posm, col = "red", lwd = 2)
legend("topright", legend = c("Position Moyenne"), col = c("red"), lty = 1, cex = 0.8)
Question 14 :
Comme on a pu voir précédemment, pour une marche aléatoire simple, la position après N pas suit une distribution binomiale avec les paramètres N (nombre de pas) et p (probabilité de chaque pas). Avec dans ce cas, p=0.5 car chaque pas a une probabilité égale d’être vers l’avant ou vers l’arrière.
Question 15 :
Nous savons que : \(E[A_N] = E[p_1 + p_2 + \ldots + p_N]\), et donc \(E[A_N^2] = E[(p_1 + p_2 + \ldots + p_N)^2]\).
On peut donc écrire que : \(E[A_N^2] = E[p_1^2 + p_2^2 + \ldots + p_N^2 + 2(p_1p_2 + p_1p_3 + p_2p_3 + \ldots + p_{N-1}p_N)]\)
Or Tous les termes dans \(2(p_1p_2 + p_1p_3 + p_2p_3 + \ldots + p_{N-1}p_N)]\) sont nuls car, comme prouvé à la question 13, tous ont une espérance mathématique égale à 0, on peut donc les éliminer.
Donc, \[ E[A_N^2] = E[p_1^2 + p_2^2 + \ldots + p_N^2] \]
Comme \(p_i\) prend soit 1 soit -1 comme valeur, \(p_i^2 = 1\).
Donc, \[ E[A_N^2] = E[1 + 1 + \ldots + 1] = N \]
En prenant les racines, on obtient l’égalité demandée : \[\sqrt{E[A_N^2]} = \sqrt{N}\].
#On prend un nombre de pas de 25
N <- 25
n_simulations <- 1000
resultats <- c(n_simulations)
for (i in 1:n_simulations) {
pos <- sum(sample(c(-1, 1), N, replace = TRUE))
resultats[i] <- sum(pos)^2
}
# Calcul de l'espérance de E[A_N^2]
esperance_A_N2 <- mean(resultats)
print(paste("Espérance de E[A_N^2] :", esperance_A_N2))## [1] "Espérance de E[A_N^2] : 24.04"print(sqrt(esperance_A_N2))## [1] 4.90306print(sqrt(N))## [1] 5On observe que pour N=25 on retrouve bien \(E[A_N^2] \approx \sqrt{N}= 25\) ce qui confirme notre démonstration. Plus on augmente le nombre de simulations plus on devrait se rapprocher de 5.
Généralisation en dimension 2
Question 16 :
\(B_0\) est la position initiale c’est donc l’origine (0,0), de la même manière que la question 11 c’est une valeur fixée donc la probabilité que \(B_0\) soit (0,0) est de 1.
Question 17 :
Pour \(B_1\), on a 4 solutions possibles : \((1, 0); (-1, 0); (0, 1);(0, -1)\) pour chacune de ces positions on a une probabilité de \(1/4\) d’y arriver.\(B_1\) représente une loi uniforme discrète, on pourrait donc facilement la simuler avec la fonction sample qui renvoit une des directions choisis.
Pour \(B_2\), on exprime la position du point après 2 pas ou chaque pas est pris dans les 4 solutions possibles:\((1, 0); (-1, 0); (0, 1);(0, -1)\) avec une probabilité égale. On peut donc ainsi obtenir : \((2, 0); (-2, 0); (1, 1);(-1,-1);(1, -1);(-1;1);(0,0)\). On peut simuler B2 tel que :
set.seed(97)#on peut changer la seed pour avoir d'autres resultats
simulB <- function() {
dir <- sample(c("avant", "arrière", "haut", "bas"), 1)
if (dir == "avant") return(c(1, 0))
else if (dir == "arrière") return(c(-1, 0))
else if (dir == "haut") return(c(0, 1))
else return(c(0, -1))
}
B1<- simulB()
B2<-B1+simulB()
B1## [1] 0 -1B2## [1] 1 -1Question 18 : A l’instar de la partie à 1 dimension, on devrait aussi ici avoir une espérance égale à 0 puisqu’on a encore une fois une symétrie. On va donc chercher à le montrer avec du code en R
# fonction utilisée question avant
simulB <- function() {
dir <- sample(c("avant", "arrière", "haut", "bas"), 1)
if (dir == "avant") return(c(1, 0))
else if (dir == "arrière") return(c(-1, 0))
else if (dir == "haut") return(c(0, 1))
else return(c(0, -1))
}
# simulation de simulB pr N pas
simulBN <- function(N, n_rep) {
# matrice de taille nrep;2 avec ttes les valeurs initialisées à 0 pr stocker les positions finales
posf <- matrix(0, n_rep, 2)
# on repete avec le nombre de repetitions qu'on a choisit
for (i in 1:n_rep) {
# on place le point en (0,0)
position <- c(0, 0)
# on génere des positions selon le nombre de pas choisis
for (j in 1:N) {
# en gros c'est un pas
deplacement <- simulB()
# On actualise la position
position <- position + deplacement
}
# Stockage de la position finale dans la matrice ce qui va nous servir à faire la moyenne
posf[i, ] <- position
}
# Retourne la matrice avec ttes les positions finales
return(posf)
}
# Nombre de pas et de repetitions choisis
N <- 100
n_repetitions <- 1000
simulation <- simulBN(N, n_repetitions)
position_moyenne <- colMeans(simulation) #colMeans ca fait la moyenne des colonnes d'une matrice
#Affichage de la moyenne littéral et graphiquement (pch=16 c'est pour faire des points remplis car sans ça ça fait des points vides)
cat("Position moyenne après", N, "pas :", position_moyenne, "\n")## Position moyenne après 100 pas : -0.204 0.144plot(simulation, col = "cadetblue1", pch = 16, main = "Position Moyenne", xlab = "X", ylab = "Y")
points(position_moyenne[1], position_moyenne[2], col = "brown1", pch = 16)
Comme on peut voir sur le graphique, la position moyenne est en (0,0) ce qui ce qu’on avait déduit théoriquement.
Question 19 :
On peut maintenant simplement réutiliser nos fonctions que l’on a écrite pour simuler \(B_N\) avec \(N=10\) et \(N=25\):
simulB <- function() {
dir <- sample(c("avant", "arrière", "haut", "bas"), 1)
if (dir == "avant") return(c(1, 0))
else if (dir == "arrière") return(c(-1, 0))
else if (dir == "haut") return(c(0, 1))
else return(c(0, -1))
}
simulBN <- function(N, n_rep) {
posf <- matrix(0, n_rep, 2)
for (i in 1:n_rep) {
position <- c(0, 0)
for (j in 1:N) {
deplacement <- simulB()
position <- position + deplacement
}
posf[i, ] <- position
}
return(posf)
}
N <- 10
n_repetitions <- 100
simulation <- simulBN(N, n_repetitions)
plot(simulation, col = "darkorchid", pch = 16, main = "", xlab = "X", ylab = "Y")
position_moyenne <- colMeans(simulation)
points(position_moyenne[1], position_moyenne[2], col = "brown1", pch = 16)
N <- 25
n_repetitions <- 100
simulation <- simulBN(N, n_repetitions)
plot(simulation, col = "darkorchid", pch = 16, main = "", xlab = "X", ylab = "Y")
position_moyenne <- colMeans(simulation)
points(position_moyenne[1], position_moyenne[2], col = "brown1", pch = 16)
On remarque que pour N=25 on peut obtenir des points un peu plus éparpillés que pour N=10 mais la moyenne est encore une fois proche ou égale à 0.
Question 20 :
On sait que : \(Var(B_N) = E((B_N)^2) - (E(B_N))^2\)
Mais \((E(X))=0\) comme montré dans les questions précédentes donc \(Var(B_N) = E((B_N)^2)\)
On a donc \(\sqrt{(Var(B_N)} = \sqrt{E((B_N)^2)}\).
Or si on suppose que chaque pas est indépendants et identiquement distribués la variance de chaque composante est la même et vaut 1 : \(Var(y_i)=Var(x_i)=1\) puisqu’on réalise un déplacement unitaire (qui vaut 1).
Enfin, \(Var(B_N) = N⋅Var(x_i) = N⋅Var(y_i)=N⋅1\) Donc : \[\sqrt{(Var(B_N)} = \sqrt{E((B_N)^2)}=\sqrt{N}\]
On devrait retrouver comme pour la question 15 avec la distance moyenne quadratique et donc\(\sqrt{E[||B_N||_2^2]} = \sqrt(N)\).Pour 25 pas on devrait alors avoir \(\sqrt{E[||B_N||_2^2]}=5\)
# fonctions utilisées question avant
simulB <- function() {
dir <- sample(c("avant", "arrière", "haut", "bas"), 1)
if (dir == "avant") return(c(1, 0))
else if (dir == "arrière") return(c(-1, 0))
else if (dir == "haut") return(c(0, 1))
else return(c(0, -1))
}
simulBN <- function(N, n_rep) {
posf <- matrix(0, n_rep, 2)
for (i in 1:n_rep) {
position <- c(0, 0)
for (j in 1:N) {
deplacement <- simulB()
position <- position + deplacement
}
posf[i, ] <- position
}
return(posf)
}
n_rep <- 1000
# Simulation pour N=25
positions_N25 <- simulBN(25, n_rep)
# Calcul de ||B_N||^2, espérance et racine de l'espérance pour N=25
distance_N25 <- rowSums(positions_N25^2)
esperance_distance_N25 <- mean(distance_N25)
racine_esperance_N25 <- sqrt(esperance_distance_N25)
cat("Espérance de ||B_N||^2 pour N=25:", esperance_distance_N25, "\n")## Espérance de ||B_N||^2 pour N=25: 25.38cat("√Espérance de ||B_N||^2 pour N=25:", racine_esperance_N25, "\n")## √Espérance de ||B_N||^2 pour N=25: 5.037857On retrouve bien \(\sqrt{E[||B_N||_2^2]}\approx5\), plus on augmente le nombre de répétitions plus on devrait tendre vers 5.