L'objectif de cette fiche est de se familiariser avec la loi normale.
Commençons par tracer la densité de la la loi normale. Voici une première version avec les outils de base:
plot(dnorm, xlim = c(-5, 5))
Et une deuxième version avec ggplot2:
library(ggplot2)
ggplot(data.frame(x = c(-5, 5)), aes(x)) + stat_function(fun = dnorm, args = list(mean = 0,
sd = 1)) + theme_bw()
J'espère qu'il est bien clair que l'histogramme d'un échantillon est une approximation de la fonction de densité.
Traçon la même chose en ajoutant avec la fonction de répartition (i.e., la cumulative):
ggplot(data.frame(x = c(-5, 5)), aes(x)) + stat_function(fun = dnorm) + stat_function(fun = pnorm,
color = "red") + theme_bw()
De même l'équivalent de la fonction de répartition pour un échantillon est la fonction ecdf.
Et enfin, l'inverse de la fonction de répartition (notez bien le domaine de définition):
ggplot(data.frame(x = c(0, 1)), aes(x)) + stat_function(fun = qnorm) + theme_bw()
Évidemment, pour une loi non bornée, on a deux asymptotes en 0 et en 1…
Maintenant, même si ces variables aléatoires ne sont pas bornées, de grandes déviations de la moyenne restent très très rares:
pnorm(1) - pnorm(-1)
## [1] 0.6827
pnorm(2) - pnorm(-2)
## [1] 0.9545
pnorm(3) - pnorm(-3)
## [1] 0.9973
Wow. Moins de 3 pour mille au delà de 3…
x = rnorm(1000)
y = rnorm(1000)
“Vérifions” la linéarité.
ggplot(data.frame(x = c(4 + 3 * x)), aes(x)) + geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm, args = list(mean = 4, sd = 3)) + theme_bw()
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
“Vérifions” la convolution:
ggplot(data.frame(x = c(4 + 3 * x + 2 * y)), aes(x)) + geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm, args = list(mean = 4, sd = sqrt(9 + 4))) + theme_bw()
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Et tentons une comparaison différente de celle de l'histogramme, la comparaison des quantiles:
ggplot(data.frame(x = c(4 + 3 * x + 2 * y)), aes(sample = x)) + stat_qq() +
theme_bw() + coord_fixed() + geom_vline(xintercept = 0) + geom_hline(yintercept = 0)
n = 1000
x = seq(0, 0, length.out = n)
for (i in (1:12)) {
x = x + runif(n)
}
x = x - 6
Regardons l'histogramme:
ggplot(data.frame(x = x), aes(x)) + geom_histogram(aes(y = ..density..)) + stat_function(fun = dnorm) +
theme_bw()
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Et comparons les quantiles:
ggplot(data.frame(x = x), aes(sample = x)) + stat_qq() + theme_bw() + coord_fixed() +
geom_vline(xintercept = 0) + geom_hline(yintercept = 0)
Surprenant, hein ? C'est la magie du théorème central limite que nous verrons la prochaine fois…
Tirons donc des points au hasard dans le plan:
x = runif(10000, min = -1, max = 1)
y = runif(10000, min = -1, max = 1)
ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) + stat_bin2d() + coord_fixed() +
geom_vline(xintercept = 0) + geom_hline(yintercept = 0)
Pas de surprise, c'est uniforme dans le carré \( [0,1]\times[0,1] \). Que se passe-t-il si on fait pareil avec la loi uniforme ?
x = rnorm(10000)
y = rnorm(10000)
df = data.frame(x = x, y = y)
ggplot(df, aes(x = x, y = y)) + stat_bin2d() + coord_fixed() + geom_vline(xintercept = 0) +
geom_hline(yintercept = 0)
Étonnant, c'est circulaire… C'est donc invariant par rotation! Et donc, l'angle en polaire est uniforme dans \( [0,2\pi] \). Et le carré du rayon, à quoi ressemble-t-il ?
df$r = sqrt(df$x^2 + df$y^2)
ggplot(df, aes(x = r^2)) + geom_histogram(aes(y = ..density..)) + stat_function(fun = dexp,
args = list(rate = 1/2)) + theme_bw()
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Donc soient \( U_1 \) et \( U_2 \) deux variables aléatoires indépendantes uniformément distribuées dans [0,1]. Définissons \( Z_0 = R \cos(\Theta) =\sqrt{-2 \ln U_1} \cos(2 \pi U_2)\ \) et \( Z_1 = R \sin(\Theta) = \sqrt{-2 \ln U_1} \sin(2 \pi U_2).\, \). Alors \( Z_0 \) et \( Z_1 \) sont des variables aléatoires indépendantes suivant une loi normale de variance 1.
u1 = runif(10000)
u2 = runif(10000)
df = data.frame(u1 = u1, u2 = u2)
df$r = sqrt(-2 * log(df$u1))
df$theta = 2 * pi * df$u2
df$x = df$r * cos(df$theta)
df$y = df$r * sin(df$theta)
Regardons l'histogramme:
ggplot(data.frame(x = c(df$x, df$y)), aes(x)) + geom_histogram(aes(y = ..density..)) +
stat_function(fun = dnorm) + theme_bw()
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Pas mal. Cependant, on utilise peu cette méthode en pratique car elle utilise des fonctions trigonométriques, coûteuses en temps de calcul. On préfère des méthodes de rejet…
