Motivation : D’un problème pratique à la modélisation

Example : Packaging A ou packaging B

  • On demande à des consommateurs s’ils préfèrent, pour un produit de grande consommation qu’on veut relooker, le packaging A ou le packaging B.

  • On interroge n personnes dans un panel de consommateurs. Les résultats sont les suivants

\[ A, A, B, A, B, B, . . .\]

  • Problème métier : choisir entre deux packaging. Choisir le packaging qui se vend le mieux !

  • Idée : se baser sur des données pour prendre la décision

Quel est le travail du statisticien ?

  • Donner un modèle probabiliste
  • Proposer une méthodologie
  • Répondre au problème métier à partir des résultats obtenus en prenant en compte l’incertitude (on ne pas sure à 100 %, il faut prendre des risques et accepter de ne pas avoir toujours raison)

Codage : on code la préférence pour A par 0 et celle pour B par 1.

\[ 0, 0, 1, 0, 1, 1, . . .\]

La réponse d’un consommateur choisi au hasard est vue comme la valeur d’une variable aléatoire de loi Bernoulli.

On dit que la variable \(X\) suit une loi Bernoulli si\[P(X = 1) = p \quad {\rm{et}} \quad P(X = 0) = 1-p,\] où \(p\) est un nombre réel compris entre 0 et 1 appelé paramètre de la loi. On note cette loi \(\mathcal{B}(p)\).

  • Dans notre exemple, le paramètre \(p\) correspond à la proportion de consommateurs dans le grand panel en faveur de B.

  • Si le panel de consommateurs est vraiment très grand, les enquêteurs n’ont pas le temps d’interroger tout le monde et ne peuvent accéder à la vraie valeur de \(p\).

  • Il est important de remarquer que ici \(p\) est inconnu.

Exercice supposons qu’en 100 consommateurs sondés, on ait eu 60 votes en faveur de A et 40 en faveur de B. À partir de cet échantillon,

Quel est l’estimateur naturel de \(p\) ? Est-ce qu’on peut être confiant de notre réponse ? Est-ce qu’on peut dire que l’estimation s’améliore si on augmente le nombre des consomateurs ? Est-ce qu’on peut avoir une idée de la précision ?

Je prétends que la vraie valeur de p est 0.43. Est-ce que cela vous parait raisonnable? A nous de jouer !

1. Simulations

La simulation de données consiste typiquement à tirer un échantillon \((X_1,\dots,X_n)\) de variables indépendantes et identiquement distribuées (en abrégé i.i.d.) suivant la même loi qu’une variable aléatoire \(X\). Les fonctions utilisées en R pour simuler des données suivant une loi s’écrivent avec le préfixe r. Par exemple :

Echantillonage d’une loi Bernoulli

  1. Fixons p=0.43
p <- .43
X<- rbinom(1, 1, p)
X
## [1] 0

  1. Exécuter plusieurs fois les deux lignes précedentes. Que peut-on remarquer ? Regarder le résultat de votre voisin. Que peut-on remarquer ?

  2. Taper les trois lignes suivantes sur votre console et comparer avec votre voisin. Que peut-on remarquer ?

set.seed(42)
X<- rbinom(1, 1, p)
X
## [1] 1
  1. On interroge 30 personnes. Les résultats obtenus pour chaque personne sont donnés ci-dessous
X <- rbinom(30, 1, p)
X
##  [1] 1 0 1 1 0 1 0 1 1 0 1 1 0 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 1 1
  1. Representer la variable X et calculer la moyenne observée
tab_X <- table(X)
tab_X
## X
##  0  1 
## 14 16
barplot(tab_X)

mean(X)
## [1] 0.5333333

Regarder le résultat de votre voisin. Que peut-on remarquer ?

  1. On a répété 20 fois la procédure précedente et nous l’avons stocké dans un tableau nommé sims.rds. Nous avons affiché les résultats obtenus à l’aide d’un graphique.

Que peut on observer ?

  1. Pour chaque essai, on va calculer la moyenne empirique. Stocker les valeurs dans une variable. Faire une répresentation graphique de la variable
sims1 <- readRDS("sims1.RDS")

sims1_mean <- tapply(sims1$sim,INDEX=sims1$rep,FUN=mean)
sims1_mean
##         1         2         3         4         5         6         7 
## 0.4333333 0.3666667 0.5000000 0.2666667 0.5333333 0.3333333 0.5000000 
##         8         9        10        11        12        13        14 
## 0.5000000 0.4000000 0.4333333 0.5333333 0.4666667 0.5333333 0.4333333 
##        15        16        17        18        19        20 
## 0.2666667 0.6000000 0.4000000 0.4666667 0.4000000 0.5000000
table(round(sims1_mean,2))
## 
## 0.27 0.33 0.37  0.4 0.43 0.47  0.5 0.53  0.6 
##    2    1    1    3    3    2    4    3    1
barplot(table(round(sims1_mean,2)))

# 100 repetitions

sims2 <-readRDS("sims2.RDS")

sims2_mean <- tapply(sims2$sim,INDEX=sims2$rep,FUN=mean)
table(round(sims2_mean,2))
## 
##  0.2 0.27  0.3 0.33 0.37  0.4 0.43 0.47  0.5 0.53 0.57  0.6 0.63 
##    1    2    7    5   10   21   12   15   10    9    5    2    1
barplot(table(round(sims2_mean,2)))

# 1000 repetitions

sims3 <-readRDS("sims3.RDS")

sims3_mean <- tapply(sims3$sim,INDEX=sims3$rep,FUN=mean)
table(round(sims3_mean,2))
## 
## 0.13  0.2 0.23 0.27  0.3 0.33 0.37  0.4 0.43 0.47  0.5 0.53 0.57  0.6 0.63 
##    1    8   17   30   43   91  116  148  153  131   94   72   48   27   11 
## 0.67  0.7 0.73 
##    5    4    1
barplot(table(round(sims3_mean,2)))

  1. On a répété 1000 fois la procédure précedente pour des échantillons de taille n=10, 30, 50, 100, 1000,5000. et nous l’avons stocké dans un tableau nommésims_mean_n.RDS. On a présenté les 6 figures dans une même fenêtre graphique. Comparer les résultats obtenus et interpréter.

L’étude de somme de variables aléatoires indépendantes et de même loi joue un rôle important en statistique.

La loi des grands nombres assure la convergence en probablilité de la moyenne empirique d’une suite de variables i.i.d. vers l’espérance de la loi commune.

Application TCL

On va d’illuster le théorème centrale limite : une somme de variables aléatoires indépendantes et de même loi renormalisé converge vers une variable aléatoire de loi normale.

La simulation proposée illustre le résultat fondamental du théorème de la limite centrale. Commenter avec votre enseignant le graphique ci-dessous.

Loi normale, \(\mathcal{N}(\mu,\sigma)\)

  1. Représenter la densité de la loi normale centrée réduite (on représentera ces fonctions sur l’intervalle \([-4;4]\)).
X <- seq(-4,4,0.1)
pdf <- dnorm(X)
plot(pdf, col='red', type = "l", lwd=2)
title(main='Densité de la loi normale')

  1. Que font ces lignes de commande ? Indiquer ce que R vous retorne.
qnorm(0.975)
## [1] 1.959964
dnorm(0)
## [1] 0.3989423
pnorm(1.96)
## [1] 0.9750021

Echantillonage d’une loi normale

  1. Que fait la commande suivante ?
rnorm(10,mean = 4,sd = 0.5)
##  [1] 4.440896 4.241102 4.482876 3.592715 4.141979 3.919151 4.967786
##  [8] 4.861615 4.179201 4.151215

Comparer le résultat obtenu avec votre voisin. Que peut-on remarquer ?

  1. Taper les deux commandes suivantes sur votre console et comparer avec votre voisin. Que peut-on remarquer ?
set.seed(42)
rnorm(10,mean = 4,sd = 0.5)
##  [1] 4.685479 3.717651 4.181564 4.316431 4.202134 3.946938 4.755761
##  [8] 3.952670 5.009212 3.968643

La fonction rnorm() dans R nous permet de générer un échantillon aléatoire de taille n d’une loi normale de moyenne \(\mu\) et d’écart type \(\sigma\).

  1. Ragarder que fait la commande help(rnorm)

Estimation

  1. Générer n variables aléatoires d’une variable X selon une loi normale de moyenne \(\mu = 35\) et d’écart-type \(\sigma=5\).
n <- 5
# génération
X <- rnorm(n,35,5)
  1. Donner une estimation de la moyenne et de l’écart-type.
# moyenne
mean(X)
## [1] 36.79054
# écart-type
sd(X)
## [1] 7.208785

Que peut-on remarquer ?

  1. Générer n = 5000 observations de la loi \(\mathcal{N}(\mu = 35; \sigma = 5)\). Puis calculer la moyenne et l’écart-type de votre echantillon.
n <- 5000
# génération
X <- rnorm(n,35,5)
mean(X)
## [1] 34.91985
sd(X)
## [1] 5.020453
  1. Produire un histogramme de densité (pour n=20, n=100, n=5000) et superimposer sur le même graphique la densité estimé et la fonction de densité de la normale
n <- 5000
# génération
X <- rnorm(n,35,5)
# histogramme de la densité
hist(X, probability=T, 
     col='lightblue') 
# estimation  de la loi par la méthode du noyau
lines(density(X), col="green", lwd=2)
# loi théorique
x=1:100
curve(dnorm(x,mean=35,sd=5),add=TRUE,col="red",lwd=2)

Echantillonage d’une loi Binomiale

  1. On jette dix fois de suite une pièce et on s’intéresse à la variable aléatoire \(X\) correspondant au nombre de Pile obtenus sur les 10 lancers .
X <- rbinom(1,10,0.43)
X
## [1] 2
  1. Quelle est la loi de \(X\) ?
  2. Représenter la loi de la variable \(X\).
  1. Que font les commandes suivantes
X <- rbinom(30,10,0.43)
X
##  [1] 5 4 7 2 3 5 3 3 6 5 3 6 6 5 2 6 3 4 6 4 4 5 3 4 2 4 3 3 5 6
  1. Donner une estimation de la moyenne et l’ecart-type
mean(X)
## [1] 4.233333
sd(X)
## [1] 1.406471
  1. Répéter 1000 fois la procédure précedente et donner une répresentation graphique
k <- 1000
n <- 10
p <- .43
X <- rbinom(k,n,p)
hist(X, 
     xlim=c(min(X),max(X)), probability=T, nclass=max(X)-min(X)+1, 
     col='lightblue',
     main='Loi binomiale, n=10, p=.5')

Approximation de la loi binomiale par la loi normale.

A nous de jouer ! Code à discuter avec votre enseignant.

p <- 0.43
n <- 100
K <- 10000
par(mfrow=c(1,2))
Y_n <- rbinom(K,n,p)
range(Y_n)
## [1] 24 64
hist(Y_n,probability=TRUE)
grille_x <- seq(range(Y_n)[1],range(Y_n)[2],by=0.01)
lines(grille_x,dnorm(grille_x,
                     mean=n*p,sd=sqrt(n*p*(1-p))),
      col="red",lwd=2)
Z_n <- (Y_n-n*p)/(sqrt(n*p*(1-p)))
grille_x <- seq(-4,4,by=0.01)
plot(density(Z_n))
lines(grille_x,dnorm(grille_x,0,1),type="l",col="red")