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)