TP2: Simulation de la loi normale

Introduction

Dans le logiciel R, on peut simuler très simplement à peu près toutes les lois classiques. Dans ce TP on propose d’étudier la loi normale: simuler des échantillons, calculer leurs résumés numériques, ainsi que représenter graphiquement leur lois.

En général, pour une loi usuelle, on dispose des fonctions suivantes:

  • rloi : permet d’echantillonner/simuler selon la loi
  • dloi : permet de calculer la densité de la loi
  • ploi : permet de calculer la fonction de répartition de la loi
  • qloi : permet de calculer le quantile de la loi

Pour reproduire les résultats d’une expérience, il est recommandé de choisir une graine aléatoire:

set.seed(13)

Echantillonner:

n <- 30
ech <- rnorm(n,mean=0, sd=2) # simulation d'un échantillon ~ N(0, 2^2)
head(ech,n = 20)
##  [1]  1.1086539 -0.5605439  3.5503267  0.3746402  2.2850523  0.8310523
##  [7]  2.4590131  0.4733593 -0.7307655  2.2102885 -2.1871879  0.9237418
## [13] -2.7219691 -3.7120543 -0.8797108 -0.3878938  2.7928630  0.2013265
## [19] -0.2288776  1.4044505

Densité:

# Densité:
hist(ech, freq=F, xlim=c(-6,6),
     ylab="densité", xlab="",
     main="Histogramme / Densité")
curve(dnorm(x, mean = 0,sd = 2), 
      from=-6, to=6, ,add=T,
      col="red", lw=2)

Fonction de repartition:

plot(ecdf(ech), xlim=c(-6,6), main="Répartition")  # Fonction de repartition empirique
curve(pnorm(x,mean = 0,sd = 2),
      from=-6,to=6, lw=2,
      add=T,col="red")

Quantiles empiriques vs theoriques.

quantile.normale <- qnorm(p=c(0.05,0.1,0.25,0.5,0.75,0.9,0.95),
                          mean=0, sd=2)
quantile.echantillon <- quantile(x=ech,
                                 probs=c(0.05,0.1,0.25,0.5,0.75,0.9,0.95))
quantiles = rbind(quantile.echantillon,quantile.normale)
5% 10% 25% 50% 75% 90% 95%
quantile.echantillon -3.36 -2.74 -1.31 0.34 1.21 2.49 3.21
quantile.normale -3.29 -2.56 -1.35 0.00 1.35 2.56 3.29

Exercice 1: Distribution des valeurs

Pour simuler \(n\) valeurs aléatoirement suivant la loi normale de moyenne \(m\) et d’écart-type \(s\) on utilise sous R la commande

rnorm(20,0,1) # simulation d'un échantillon de taille 20 d'une loi N(0,1)
rnorm(20,1,1) # simulation d'un échantillon de taille 20 d'une loi N(1,1)
rnorm(20,2,1) # simulation d'un échantillon de taille 20 d'une loi N(2,1)
rnorm(20,0,2) # simulation d'un échantillon de taille 20 d'une loi N(0,1)
rnorm(20,0,4) # simulation d'un échantillon de taille 20 d'une loi N(0,4)

Pour chacun des échantillons ci-dessus remplir le tableau suivant :

\(n=20\) \(\mu=0\), \(\sigma=1\) \(\mu=1\), \(\sigma=1\) \(\mu=2\), \(\sigma=1\) \(\mu=0\), \(\sigma=2\) \(\mu=0\), \(\sigma=4\)
% de valeurs < -2
% de valeurs < 0
% de valeurs = 0
% de valeurs >2

Pour toute la suite du TP, voici un script R qui vous guidera pour répondre aux questions posées.

# Echantillons de taille n=2000 individus avec differents parametres
n <- 2000
nd1 <- rnorm(n,0,1)
nd2 <- rnorm(n,1,1)
nd3 <- rnorm(n,2,1)
nd4 <- rnorm(n,0,2)
nd5 <- rnorm(n,0,4)


nd_en_col <- cbind(nd1, nd2, nd3, nd4, nd5) # les echantillons sont mis en ensemble, colonne par colonne
# pour chaque colonne on compte le nombre de valeurs qui satisfont une condition + pourcentage
colSums(nd_en_col< -2) / n              
colSums(nd_en_col< 0) / n 
colSums(nd_en_col == 0) / n 
colSums(nd_en_col > 2) / n 

Question 1 : avec les mêmes valeurs des moyennes et des écart-types, on simule maintenant des échantillons de taille \(n=2000\) individus.

Pour chacun de ces échantillons remplir le tableau suivant en utilisant les commandes R dans le script :

\(n=20\) \(\mu=0\), \(\sigma=1\) \(\mu=1\), \(\sigma=1\) \(\mu=2\), \(\sigma=1\) \(\mu=0\), \(\sigma=2\) \(\mu=0\), \(\sigma=4\)
% de valeurs < -2
% de valeurs < 0
% de valeurs = 0
% de valeurs >2

Question 2 : On considère uniquement le premier échantillon nd1. Que trace la fonction plot dans le script ? (axe des \(x=\) ? et axe des \(y=\) ?)

# Tracer les valeurs
plot(nd1, pch=19, cex=0.5, xlab="", ylab="")

Pour dessiner la répartition des valeurs on trace un histogramme.

Question 3 : Quelle est la différence entre les trois commandes ci-dessous?

hist(nd1)
hist(nd1,freq=TRUE,breaks=40)
hist(nd1,freq=TRUE,breaks=50)

Question 4 : on veut comparer les histogrammes des \(5\) échantillons. Lancer les commandes de Introduction/Densité qui permettent d’afficher les histogrammes et de superposer en rouge la densité théorique associée (pour mieux voir les graphiques cliquer sur zoom dans le menu Plots).

Comparer les graphiques des échantillons nd1, nd2 et nd3 : qu’est-ce qui change ?

Et puis comparer les graphiques des échantillons nd1, nd4 et nd5 : qu’est-ce qui change ?

Question 5 : on veut maintenant comparer les fonctions de répartition des \(5\) échantillons. Lancer les commandes de Introduction/Repartition pour les 5 échantillons.

Comparer les graphiques des échantillons nd1, nd2 et nd3 : qu’est-ce qui change ?

Et puis comparer les graphiques des échantillons nd1, nd4 et nd5 : qu’est-ce qui change ?

Question 6 : on veut montrer le lien entre la densité et la fonction de répartition.

Donner ci-dessous la densité de la loi normale de moyenne \(0\) et de variance \(1\) :

\(f(x) =\)

Quel est le lien entre la densité et la fonction de répartition ?

\(f(x)=\)

\(F(t)=\)

La valeur au point \(t\) de la fonction de répartition (théorique) \(F(t)\) d’une loi gaussienne de moyenne \(m\) et d’écart-type \(s\) est donnée par la commande pnorm(t,m,s) de R.

Lancer le script suivant et décrire ce qu’il permet d’observer.

# une fonction pour colorier l'aire sous la courbe d'une fonction <fun> entre <xmin> et <xmax>:
shade_under_curve <- function(fun, xmin, xmax, length=100){
  xvals <- seq(xmin, xmax, length=length)
  dvals <- match.fun(fun)(xvals)
  polygon(c(xvals,rev(xvals)),c(rep(0,length),rev(dvals)),col="gray")
}

curve(dnorm,from=-4,to=+4)
shade_under_curve(dnorm,-5,-1)

Question 7 : on a tracé ci dessous la courbe de densité d’une loi gaussienne centrée (de moyenne 0) et de variance \(4\) et l’on a hachuré différentes aires sous la courbe.

Compléter le tableau ci-dessous :

Aire de la figure Proba représentée par l’aire Ecriture avec \(F(t)\) Valeur numérique pnorm
1 \(P(X<-1)\) \(F(-1)\)
2
3
4

Question 8 : soit \(X \sim \mathcal N(15, 9)\) (une variable aléatoire gaussienne de moyenne \(15\) et de variance \(9\)).

Quelle est la probabilité que \(X\leq 17\) ? (utiliser la fonction pnorm)

Quelle est la probabilité que \(X>7\) ? Enfin quelle est la probabilité que \(7<X<17\) ?

Exercice 2: Lien entre lois normales

Soit \(X \sim \mathcal N(0, 1)\). On souhaite montrer que les variables \(X_1=X+1\), \(X_2=X-1\), \(X_3=2X\) et \(X_4=X/2\) suivent aussi une loi normale. Pour cela on simule un grand nombre de valeurs suivant la loi normale \(N(0,1)\) et l’on construit les nouvelles variables. Étudier la distribution de ces nouvelles variables (avec des histogrammes puis en fonction d’une distribution usuelle).

x<-rnorm(10000,0,1)
x1<-x+1
x2<-x-1
x3<-2*x
x4<-x/2

Soit maintenant \(Y \sim \mathcal N(4, 3)\) quelle est la distribution de la variable \(Z=\frac{Y-2}{\sqrt{3}}\) ?

Exercice 3: Sommes de variables aléatoires gaussiennes

En vous aidant des résultats du script suivant, donner la loi de \(X+Y\) et de \(X-Y\) si \(X\sim N(m_1,\sigma_1^2)\) et \(Y\sim N(m_2,\sigma_2^2)\).

# Deux variables gaussiennes de meme loi
x <- rnorm(10000,0,1)
y <- rnorm(10000,0,1)
z <- x+y
hist(z, breaks=20, freq=FALSE, ylim=c(0,0.5), main="Z=X+Y | Même loi")
curve(dnorm, from=-4, to=4,
      add=TRUE, col="red",
      lw=2, xaxs="i") # la loi N(0,1)
x<-seq(-5,5,0.1)
lines(x,dnorm(x,0,sqrt(2)), col="blue", lw=2) #la loi proposée

# Deux variables gaussiennes de meme ecart type mais pas de meme esperance
x <- rnorm(10000,0,1)
y <- rnorm(10000,1,1)
z <- x + y
hist(z, breaks=20, freq=FALSE, ylim=c(0,0.5), main="Z=X+Y | Même variance")
curve(dnorm,from=-4, to=4, add=TRUE, col="red", lw=2)  # la loi N(0,1)
x <- seq(-5,5,0.1)
lines(x, dnorm(x,0,sqrt(2)), col="darkgreen", lw=2) #la 1ere loi proposée
# votre proposition ici:
# lines(x,dnorm(x, ??, ??) ,col="blue", lw=2) 

# Deux variables gaussiennes de meme esperance mais pas de meme ecart type
x<-rnorm(10000,0,1)
y<-rnorm(10000,0,2)
z<-x+y
hist(z, breaks=20, freq=FALSE, ylim=c(0,0.5), main="Z=X+Y | Même esperance")
curve(dnorm,from=-4, to=4, add=TRUE, col="red", lw=2)
x<-seq(-5,5,0.1)
## votre proposition ici:
## lines(x,dnorm(x, ??, ??),col="blue",lw=2) 

En général, pour \(X\sim N(m_1,\sigma_1^2)\) et \(Y\sim N(m_2,\sigma_2^2)\):

\(X+Y\sim\)

\(X-Y\sim\)

Déduire de ce qui précède la loi de la variable \(\bar{X}=\frac{1}{n}(X_1+X_2+X_3+...+X_n)\) où toutes les variables \(X_i\) sont gaussiennes centrées réduites. Écrire un script pour l’illustrer.

Votre script ici

Exercice 4: Somme de carrés de variables gaussiennes centrées réduites

Pour une somme de carrés de variables gaussiennes centrées réduites regarder la suite du script. Trouver le nom de la nouvelle loi des carrés et comment on choisit ses paramètres. (consulter l’aide de R)

N <- 10000
mu <- 0
sigma <- 1
v1 <- rnorm(N, mu, sigma)^2
v2 <- rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2
v3 <- rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2
v4 <- rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2 + rnorm(N, mu, sigma)^2

par(mfrow=c(2,2))
hist(v1, breaks=40, prob=TRUE, xlim=c(0,10), ylim=c(0,1.5))
curve(dchisq(x,1), from=0, to=10, add=TRUE, col="red")

hist(v2, breaks=40, prob=TRUE, xlim=c(0,10), ylim=c(0,0.5))
curve(dchisq(x,2),from=0, to= 10, add=TRUE, col="red")

hist(v3, breaks=40, prob=TRUE, xlim=c(0,10), ylim=c(0,0.3))
curve(dchisq(x,3), from=0, to=10, add=TRUE, col="red")

hist(v4, breaks=40, prob=TRUE, xlim=c(0,20), ylim=c(0,0.20))
curve(dchisq(x,4), from=0, to=10, add=TRUE, col="red")      

Question : Compléter la formule ci-dessous ?

\(X_1^2+X_2^2+X_3^2+...+X_n^2\sim\)