Statistique et Probabilités [M3201]
TP 5 : Intervalles de confiance
Introduction
Pour commencer, lisez la partie du cours relative aux intervalles de confiance.
Exercices
Exercice 1
Première Partie
On considère un échantillon de taille \(n\) tiré d’une loi normale, \(\mathcal N(\mu,\sigma^2)\). ( rnorm(n, mu, sigma)
).
On s’intéresse à l’intervalle de confiance de la moyenne \(\bar X=\frac{1}{n}\sum_{i=1}^nX_i\).
Pour cela, nous allons donner les trois intervalles suivants (définis dans le cours sur l’estimation):
En considérant que la variance est connue (ici \(\sigma=1\)). \[\rm{IC_1} = \left[ \bar X - u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}};\bar X + u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}} \right] \]
En considérant que la variance est inconnue et donc en se servant du quantile de la loi de Student avec une estimation de la variance \(S\).
\[ S^2=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)^2 \]
\[ \rm{IC_2} = \left[\bar X - t_{1-\alpha/2} \frac{S}{\sqrt{n}};\bar X + t_{1-\alpha/2} \frac{S}{\sqrt{n}} \right] \]
- En considérant que la variance est inconnue mais en remplaçant le quantile de la loi de Student par celui de la loi normale. \[ \rm{IC_3} = \left[\bar X - u_{1-\alpha/2} \frac{S}{\sqrt{n}};\bar X + u_{1-\alpha/2} \frac{S}{\sqrt{n}} \right] \]
Avec \(u_{1-\alpha/2}\) le quantile d’ordre \(1 − \alpha/2\) de la loi normale \(\mathcal N(0, 1)\) (qnrom
). et \(t_{1-\alpha/2}\) le quantile d’ordre \(1-\alpha/2\) de la loi de Student de paramètre \(n-1\) (qt
).
Remplir le tableau suivant:
ICs <- function(n, alpha, mu){
x <- rnorm(n, mean=mu, sd=1)
barx <- mean(x)
# IC1
err1 <- qnorm(1-alpha/2) / sqrt(n)
IC1 <- sprintf('[%.1e, %.1e, %.1e]', barx-err1, barx, barx+err1)
# IC2
s2 <- mean((x - barx)^2)
err2 = qt(1-alpha/2, df=n-1) /sqrt(n)* sqrt(s2)
IC2 <- sprintf('[%.1e, %.1e, %.1e]', barx-err2, barx, barx+err2)
# IC3
err3 = qnorm(1-alpha/2) /sqrt(n)* sqrt(s2)
IC3 <- sprintf('[%.1e, %.1e, %.1e]', barx-err3, barx, barx+err3)
return(rbind(IC1, IC2, IC3))
}
ns <- c(100, 1000)
alpha <- c(0.05, 0.01)
mu <- 0
result <- rbind(ICs(ns[1], alpha, mu), ICs(ns[2], alpha, mu))
colnames(result ) <- paste('alpha=', alpha)
icnames <- paste('IC', c(1,2,3), sep='')
rownames(result) <- c(paste(icnames, 'n=', ns[1]),
paste(icnames, 'n=', ns[2]))
library(knitr)
kable(result)
alpha= 0.05 | alpha= 0.01 | |
---|---|---|
IC1 n= 100 | [-2.5e-01, -5.9e-02, 1.4e-01] | [-3.2e-01, -5.9e-02, 2.0e-01] |
IC2 n= 100 | [-2.4e-01, -5.9e-02, 1.2e-01] | [-3.0e-01, -5.9e-02, 1.8e-01] |
IC3 n= 100 | [-2.4e-01, -5.9e-02, 1.2e-01] | [-2.9e-01, -5.9e-02, 1.8e-01] |
IC1 n= 1000 | [-2.9e-02, 3.3e-02, 9.5e-02] | [-4.8e-02, 3.3e-02, 1.1e-01] |
IC2 n= 1000 | [-2.9e-02, 3.3e-02, 9.5e-02] | [-4.8e-02, 3.3e-02, 1.2e-01] |
IC3 n= 1000 | [-2.8e-02, 3.3e-02, 9.5e-02] | [-4.8e-02, 3.3e-02, 1.1e-01] |
On visualise pour \(n=1000\) les intervalles de confiance \(\rm{IC_1}\) des niveaux 0.85, 0.90 et 0.99. Quelle est la couleur associée à chaque niveau de confiance?
La largeur de l’intervalle de confiance est \(2\, u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}}\) donc
A n fixe (le meme échantillon):
Plus \(\alpha\) (le risque d’erreur) est petit, plus large est l’intervalle de confiance.
De maniere equivalente, plus \(1-\alpha\) (le niveau de confiance) est large, plus large est l’intervalle de confiance.
\[ \begin{align} \forall (\alpha_1, \alpha_2)\in(0,1)^2, \quad \alpha_1 < \alpha_2 & \implies 1- \frac{\alpha_2}2 < 1-\frac{\alpha_1}2 \\ & \implies u_{1-\alpha_2/2}< u_{1-\alpha_1/2} \\ &\implies \rm{IC}(\alpha_2) \subset \rm{IC}(\alpha_1)\\ & \implies \rm{IC}(\alpha_1)^{\text{inf}} < \rm{IC}(\alpha_2)^{\text{inf}} < \rm{IC}(\alpha_2)^{\text{sup}} < \rm{IC}(\alpha_1)^{\text{sup}}. \end{align} \]
A \(\alpha\) fixe:
Plus grand est l’echantillon, plus petite est l’erreur \(u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}}\), ainsi l’intervalle est plus etroit quand \(n\) est large.
Par contre, on ne peut pas comparer les bornes des intervalles car la valeur de \(\bar X\) change avec \(n\).
Deuxième partie
On souhaite refaire la même chose mais cette fois-ci les échantillons sont pollués (comme c’est souvent le cas dans la vie réelle) par des données extrêmes. On suppose ici qu’il y a 1% de données polluées par échantillon.
La fonction \(\texttt{pollue_ech}\)(n, mu, sigma2, sigma2_pollue, p) retourne un échantillon gaussien, \(\mathcal N(\mu,\sigma^2)\) (par défaut \(\mu=0\), \(\sigma = 1\)) avec une proportion \(p\) (p=0.01) de données polluées issues d’un échantillon gaussien de même moyenne mais de variance beaucoup plus grande (sigma2_pollue = 16).
pollue_ech <- function(n, mu=0,
sigma2=1, sigma2_pollue=16,
p=0.01){
pollue <- rbinom(n,1,p) # 0 ou 1 ~ B(p)
return(pollue * rnorm(n,mu,sqrt(sigma2_pollue))
+ (1 - pollue) * rnorm(n, mu, sqrt(sigma2)))
}
Remplir le tableau pour des échantillons pollués à 1%:
ICs <- function(n, alpha, mu){
x <- pollue_ech(n, mu, p=0.01)
barx <- mean(x)
# IC1
err1 <- qnorm(1-alpha/2) / sqrt(n)
IC1 <- sprintf('[%.1e, %.1e, %.1e]', barx-err1, barx, barx+err1)
# IC2
s2 <- mean((x - barx)^2)
err2 = qt(1-alpha/2, df=n-1) /sqrt(n)* sqrt(s2)
IC2 <- sprintf('[%.1e, %.1e, %.1e]', barx-err2, barx, barx+err2)
# IC3
err3 = qnorm(1-alpha/2) /sqrt(n)* sqrt(s2)
IC3 <- sprintf('[%.1e, %.1e, %.1e]', barx-err3, barx, barx+err3)
return(rbind(IC1, IC2, IC3))
}
ns <- c(100, 1000)
alpha <- c(0.05, 0.01)
mu <- 0
result <- rbind(ICs(ns[1], alpha, mu), ICs(ns[2], alpha, mu))
colnames(result ) <- paste('alpha=', alpha)
icnames <- paste('IC', c(1,2,3), sep='')
rownames(result) <- c(paste(icnames, 'n=', ns[1]),
paste(icnames, 'n=', ns[2]))
library(knitr)
kable(result)
alpha= 0.05 | alpha= 0.01 | |
---|---|---|
IC1 n= 100 | [-1.8e-01, 1.3e-02, 2.1e-01] | [-2.4e-01, 1.3e-02, 2.7e-01] |
IC2 n= 100 | [-1.7e-01, 1.3e-02, 1.9e-01] | [-2.2e-01, 1.3e-02, 2.5e-01] |
IC3 n= 100 | [-1.6e-01, 1.3e-02, 1.9e-01] | [-2.2e-01, 1.3e-02, 2.4e-01] |
IC1 n= 1000 | [-8.0e-02, -1.8e-02, 4.4e-02] | [-1.0e-01, -1.8e-02, 6.3e-02] |
IC2 n= 1000 | [-8.9e-02, -1.8e-02, 5.2e-02] | [-1.1e-01, -1.8e-02, 7.4e-02] |
IC3 n= 1000 | [-8.9e-02, -1.8e-02, 5.2e-02] | [-1.1e-01, -1.8e-02, 7.4e-02] |
Exercice 2
Dans cet exercice, on suppose que \(X\) suit une loi gaussienne, avec une variance connue et une espérance inconnue. Le résultat visuel de ce qui est demandé est visible sur le serveur shiny de l’Université Grenoble Alpes.
1) Simuler un échantillon de taille \(n=50\) avec une espérance \(\mu=1\) et un écart-type \(\sigma=1\). Calculer la moyenne empirique \(\bar x\). Calculer l’intervalle de confiance de \(\mu\) au niveau de confiance \(0.95\).
2) Ecrire une fonction CI_mean
qui prend comme arguments, l’échantillon, l’écart-type \(\sigma\) (connu) et le niveau de confiance \(\alpha\), et retourne les deux bornes de l’intervalle de confiance. Tester cette fonction avec l’échantillon simulé avec des niveaux de confiance \(0.90\), \(0.95\) et \(0.99\). Commenter.
3) Un intervalle de confiance ne contient pas toujours la vraie valeur. Simuler 100 échantillons de taille \(n=50\) avec une espérance \(\mu=1\) et un écart-type \(\sigma=1\) et mettre cela dans une matrice notée \(X100\). Calculer les 100 intervalles de confiance au niveau \(0.95\) en utilisant votre fonction CI_mean
.
Faire une figure avec les 100 intervalles de confiance avec les vraies valeurs. On utilisera les fonctions segments
ou matplot
.
Tracer sur le même graphique la droite d’équation \(x=\mu\) avec abline
.
4) Répéter la même procédure avec \(N=1000\), \(N=10000\) et \(N=100000\) et compter le nombre d’intervalles de confiance ne contenant pas la vraie valeur \(\mu\). Commenter.
5) Simuler un échantillon \(X_1,X_2, \ldots, X_n\) de taille \(n=1000\). Pour chaque \(i=1,\ldots,n\), calculer la moyenne progressive \[ \bar X_i = \frac1i \sum_{j=1}^iX_j \] calculer les intervalles de confiance au niveau \(0.95\%\) pour chaque moyenne \(\bar X_i\). Représenter avec des points rouges les moyennes successives calculées et en bleus les bornes des intervalles de confiance. Commenter.
# 1-2)
n <- 50
mu <- 1
sigma <- 1
alpha <- 0.05
x <- rnorm(n , mu, sigma)
CI_mean <- function(x, sigma, alpha){
xbar <- mean (x)
err <- qnorm(1-alpha/2) * sigma / sqrt(length(x))
return(xbar + c(-1, 1) * err)
# return(c(xbar - err, xbar + err))
}
rbind(CI_mean(x, sigma, 0.1),
CI_mean(x, sigma, 0.05),
CI_mean(x, sigma, 0.01))
## [,1] [,2]
## [1,] 0.9122688 1.377504
## [2,] 0.8677055 1.422067
## [3,] 0.7806090 1.509164
# 3-4)
count_CIs <- function(N, alpha=0.05){
X <- replicate(N, rnorm(n , mu, sigma))
CI <- apply(X, 2, CI_mean, sigma=sigma, alpha=alpha)
# mu <= CI[2,] : 1 si mu est inferieure a la borne sup
# mu >= CI[1,] : 1 si mu est superieure a la borne inf
# le produit des deux conditions = 1 si et seulement si les deux sont vraies
# i.e. mu dans [CI[1,], CI[2,]]
# 1 - le produit = 1 si mu est à l'exterieur de [CI[1,], CI[2,]]
# on somme le vecteur des 1/0 pour compter les cas où mu est à l'extérieur
# ou on prend la moyenne pour estimer une proba.
mu_outside <- mean(1 - (mu <= CI[2,]) * (mu >= CI[1,]))* 100
# Segments trace les segmnets sur un plot existant,
# donc on initie le plot avec plot(c(xmin, xmax), c(ymin, ymax))
plot(c(min(CI), max(CI)), c(1, N),
cex=0, xlab = "", ylab = "",
main=sprintf('N=%d | %.2f %% des IC ne contiennent pas mu=%d',
N, mu_outside, mu))
segments(x0 = CI[1,],
y0 = seq(N),
x1 = CI[2,],
y1 = seq(N),
col='blue', lw=0.1)
# Alternativement à plot+segments on peut utiliser matplot (un peu plus lente)
# matplot(x=CI, y=rbind(seq(1,N), seq(1,N)), col='blue', lw=0.1, type='l',
# main=sprintf('N=%d | %.2f %% des IC ne contiennent pas mu=%d',
# N, mu_outside, mu))
# La droite x=mu:
abline(v=mu, lty=2, lw=2)
}
par(mfrow=c(3,1))
count_CIs(1e3)
count_CIs(1e4)
count_CIs(1e5)
L’objectif est de retrouver le resultat : \(P(\mu \not\in [IC^\text{inf},IC^\text{sup}]) = \alpha\) i.e 5% des intervalles ne contiennent pas \(\mu\). On a une meilleure estimation du risque d’erreur avec un large \(N\) (le nombre total des intervalles).
# 5)
# seq_along(x) est equivalent à seq(1, length(x))
n <- 1000
x <- rnorm(n, 1, 1)
# Les moyennes progressives:
cummean <- cumsum(x) / seq_along(x)
# Les intervalles de confiance progressives:
CIs <- cbind(cummean - qnorm(1 - 0.05/2) / sqrt(seq_along(x)),
cummean + qnorm(1 - 0.05/2) / sqrt(seq_along(x)))
par(mfrow=c(1,2))
# Les moyennes bar X_i:
plot(seq_along(x), cummean, 'p', col='red', cex=0.4, pch=19,
ylim=c(min(CIs), max(CIs)),
main='Les moyennes progressives\n
et leurs intervalles de confiance\n
en echelle lineaire',
xlab='i', ylab=bquote(bar(X[i])))
# Borne inf
points(seq_along(x), CIs[,1], col='blue', cex=0.4, pch=25)
# Borne sup:
points(seq_along(x), CIs[,2], col='blue', cex=0.4, pch=24)
# En echelle logarithmique pour mieux visualiser les effets de bord.
# Les moyennes bar X_i:
plot(log(seq_along(x)), cummean, 'p', col='red', cex=0.4, pch=19,
ylim=c(min(CIs), max(CIs)),
xlab='log(i)', ylab=bquote(bar(X[i])), main='En echelle logarithmique')
# Borne inf
points(log(seq_along(x)), CIs[,1], col='blue', cex=0.4, pch=25)
# Borne sup:
points(log(seq_along(x)), CIs[,2], col='blue', cex=0.4, pch=24)
Exercice 3
Afin d’évaluer l’impact d’une campagne média anti-tabac, on s’est intéressé à la proportion de fumeurs menant des actions pour essayer d’arrêter de fumer (diminution de la consommation, achat de patchs anti-tabac, consultations médicales, …), c’est-à-dire à la proportion de fumeurs “actifs” pour arrêter.
Un sondage “avant campagne” a été effectué auprès de \(n_1=1000\) fumeurs, et un sondage “après campagne” a été effectué auprès d’un autre échantillon de \(n_2=1000\) fumeurs ; les deux échantillons sont donc indépendants. Le premier sondage donne une proportion de \(p_1=0.15\) de fumeurs “actifs”, alors que le deuxième sondage donne une proportion de \(p_2=0.17\) de fumeurs “actifs”. On veut savoir si la campagne a été efficace ; autrement dit si la proportion de fumeurs “actifs” a augmenté après la campagne.
1) Déterminer un intervalle de confiance au niveau 95% de la proportion de fumeurs “actifs” avant la campagne.
Rappelons l’intervalle de confiance d’une proportion \(\hat p\): \[ \left[\hat p - u_{1-\alpha/2} \sqrt{\frac{\hat p (1-\hat p)}{n}};\hat p + u_{1-\alpha/2} \sqrt{\frac{\hat p (1-\hat p)}{n}} \right]. \] Avec \(\hat p\) la fréquence empirique et \(u_{1-\alpha/2}\) le quantile d’ordre \(1-\alpha/2\) de la loi normale \(\mathcal N(0,1)\).
2) De façon analogue, donner un intervalle de confiance au niveau 95% de la proportion de fumeurs “actifs” après la campagne.
3) Peut-on déduire de ces deux intervalles que la campagne a été efficace ?
4) On suppose maintenant que \(p_2=0,20\). Que pouvez-vous conclure ?
5) On suppose encore que \(p_2=0,20\) mais que cette proportion a été calculée sur un échantillon de \(n_2=500\) personnes. Que peut-on en conclure ?
alpha <- 0.05 # confiance 0.95
ICcampagne <- function(p, n){
p + c(-1, 1) * qnorm(1- alpha/2) * sqrt(p * (1-p)/n)
} # retourne en 1 la borne inf et en 2 la borne sup.
# 1) Avant la campagne
IC_avant <- ICcampagne(0.15, 1000)
# 2) Apres la campagne
IC_apres <- ICcampagne(0.17, 1000)
# >> Resultats peu concluants car IC_apres[1] < IC_avant[2].
# 4) Apres la campagne avec p2=0.2
IC_apres_p20 <- ICcampagne(0.2, 1000)
# >> la campagne etait efficace car IC_apres_20[1] > IC_avant[2].
# 5)
IC_apres_p20_n500 <- ICcampagne(0.2, 500)
# >> Resultats peu concluants car IC_apres_20_500[1] < IC_avant[2].
result <- rbind(IC_avant, IC_apres,
IC_avant, IC_apres_p20,
IC_avant, IC_apres_p20_n500)
colnames(result) <- c('Borne inf', 'Borne sup')
# Afficher le tableau:
kable(result, digits=3)
Borne inf | Borne sup | |
---|---|---|
IC_avant | 0.128 | 0.172 |
IC_apres | 0.147 | 0.193 |
IC_avant | 0.128 | 0.172 |
IC_apres_p20 | 0.175 | 0.225 |
IC_avant | 0.128 | 0.172 |
IC_apres_p20_n500 | 0.165 | 0.235 |
# Visualiser les intervalles de confiance:
matplot(t(result), rbind(seq(1, 6), seq(1,6)),
type='o', lty=1, pch=3, lwd=2,
col=c('black', 'red', 'black', 'green', 'black', 'blue'),
xlab='CI', ylab='Experience', ylim=c(0, 6.3))
text(.5*(result[1,1]+result[1,2]),
.2+c(1,3,5),
rep('Avant, p=0.15, n=1000', 3))
text(.5*(result[c(2,4,6),1]+result[c(2,4,6),2]),
.2 + c(2,4,6),
c('Apres, p=0.17, n=1000',
'Apres, p=0.2, n=1000',
'Apres, p=0.2, n=500'))