Remerciez-le!

Remerciez @Admin pour avoir partagé cet document gratuitement, de la manière la plus simple, en partageant sur les réseaux sociaux.

TP 3 – SY02 Estimation et théorème central limite Corrigé Les questions/section

TP 3 – SY02 Estimation et théorème central limite Corrigé Les questions/sections marquées par un sont des questions plus ouvertes qui peuvent être abordées dans un deuxième temps. Les questions/sections marquées par un sont des questions qui sont prévues pour être traitées en autonomie en dehors de la séance de TP. 1 Estimation par la méthode des moments On modélise un phénomène par une variable aléatoire X suivant une loi uniforme sur l’intervalle [0, a] avec a entier. Aﬁn d’estimer le paramètre a, on cherche un estimateur avec la méthode des moments. On utilisera la fonction suivante pour générer un échantillon de taille n : runifa <- function(n) { if(!exists("param")) param <<- sample(10:20, 1) runif(n, min = 0, max = param) } 1 Sachant que E(X) = a 2, créer une fonction estim qui prend en argument un échantillon de taille quelconque issu de X et renvoie l’estimation, par la méthode des moments, du paramètre a. Estimateur des moments de a. On inverse le système en exprimant les paramètres en fonction des moments. On trouve a = 2E(X). On propose donc l’estimateur b a = 2X. Fonction estim : runifa <- function(n) { if (!exists("param")) param <<- sample(10:20, 1) runif(n, min = 0, max = param) } estim <- function(x) { 2 * mean(x) } (1/14) TP 3 – SY02 Printemps 2022 2 Lancer plusieurs fois la fonction précédente avec un échantillon de taille 100. Quel semble être le paramètre a ? n <- 100 estim(runifa(n)) [1] 17.89559 estim(runifa(n)) [1] 15.31492 estim(runifa(n)) [1] 18.1151 estim(runifa(n)) [1] 17.68204 Le paramètre entier inconnu semble être ici 18 (il est diﬀérent à chaque nouvelle session). Au lieu d’exécuter manuellement plusieurs fois l’instruction estim(runifa(n)), on va utiliser une fonction R qui le fait à notre place et accumule les diﬀérents résultats : il s’agit de la fonction replicate qu’on utilise comme suit : a <- replicate(1000, estim(runifa(n))) Le vecteur a stocke les 1000 résultats de l’instruction estim(runifa(n)). 3 Faire un diagramme en boîte de 1000 estimations successives de a. Quel est le paramètre inconnu ? Vériﬁer qu’il est en accord avec le vrai paramètre choisi au hasard, utilisé pour générer les observations et stocké dans param. a <- replicate(1000, estim(runifa(n))) boxplot(a) 15 16 17 18 19 20 21 Le paramètre recherché est donc 18. param [1] 18 4 On admet que les moments de X sont : (2/14) TP 3 – SY02 Printemps 2022 mk = E(Xk) = ak k + 1 pour tout k ≥1. Trouver pour chaque moment l’estimateur ˆ ak correspondant. Refaire l’étude précédente. Quel est l’estimateur le plus précis ? Les estimateurs sont : ˆ ak = ((k + 1) ˆ mk)(1/k) estimk <- function(x, k) { ((k + 1) * mean(x^k))^(1/k) } a1 <- replicate(1000, estimk(runifa(n), 1)) a2 <- replicate(1000, estimk(runifa(n), 2)) a5 <- replicate(1000, estimk(runifa(n), 5)) boxplot(a1, a2, a5, names = c("k=1", "k=2", "k=5")) k=1 k=2 k=5 15 16 17 18 19 20 21 La variance des estimateurs baisse. Lorsque k est grand, on trouve l’estimateur max i Xi. 2 Théorème central limite On souhaite vériﬁer expérimentalement le théorème central limite. Pour cela, on va choisir une variable aléatoire X de loi quelconque, d’espérance µ et d’écart-type σ. On utilisera la fonction suivante pour générer un échantillon de loi inconnue de taille n : runknown <- function(n) { bn <- rbinom(n, 1, 0.2) bn * rnorm(n, mean=-4, sd=1) + (1 - bn) * rnorm(n, mean=10, sd=1) } 5 Vériﬁer expérimentalement que l’espérance de X vaut µ = 7.2 et que l’écart-type vaut σ = √ 32.36. (3/14) TP 3 – SY02 Printemps 2022 runknown <- function(n) { bn <- rbinom(n, 1, 0.2) bn * rnorm(n, mean = -4, sd = 1) + (1 - bn) * rnorm(n, mean = 10, sd = 1) } n <- 1000 x <- runknown(n) mean(x) [1] 7.094625 sd(x) [1] 5.775383 sqrt(32.36) [1] 5.688585 6 Tracer un histogramme d’un échantillon de taille 1000. n <- 1000 x <- runknown(n) hist(x) Histogram of x x Frequency −5 0 5 10 0 100 200 300 400 7 Tracer la fonction de répartition empirique de la loi inconnue avec un échantillon de taille 1000 à l’aide des fonctions ecdf et plot. Commenter la cohérence du graphe de la fonction de répartition avec celui de l’histogramme. plot(ecdf(x)) (4/14) TP 3 – SY02 Printemps 2022 −5 0 5 10 15 0.0 0.2 0.4 0.6 0.8 1.0 ecdf(x) x Fn(x) La fonction ecdf est un peu spéciale : elle renvoie elle-même une fonction. Pour calculer la fonction de répartition empirique en x = 2, on écrit donc ecdf(x)(2) [1] 0.216 La fonction de répartition empirique présente un palier sur un intervalle inclus dans [−3, 9] : c’est cohérent avec l’histogramme qui indique que cet intervalle ne contient pas de valeurs simulées selon la loi runknown. Le théorème central limite dit que si on a n variables aléatoires iid X1, . . . , Xn de v.a. parente X, alors la variable aléatoire suivante, T = X −µ σ/√n , converge en loi lorsque n augmente vers une loi normale centrée réduite. On va donc vériﬁer que des réalisations issues de la loi de T sont distribuées comme une gaussienne centrée réduite. 8 À partir de l’échantillon de taille n de loi celle de X, calculer une seule réalisation de la variable aléatoire T. mu <- 7.2 sigma <- sqrt(32.36) (mean(x) - mu)/(sigma/sqrt(n)) [1] -1.241367 9 Réitérer l’opération précédente plusieurs fois en regénérant à chaque fois un nouvel échantillon de taille n de X. À quoi semble ressembler la distribution obtenue ? x <- runknown(1000) (mean(x) - mu)/(sigma/sqrt(n)) [1] 0.1960422 (5/14) TP 3 – SY02 Printemps 2022 x <- runknown(1000) (mean(x) - mu)/(sigma/sqrt(n)) [1] 1.228306 x <- runknown(1000) (mean(x) - mu)/(sigma/sqrt(n)) [1] -1.070602 x <- runknown(1000) (mean(x) - mu)/(sigma/sqrt(n)) [1] 1.371143 La distribution semble se concentrer autour de 0 avec un faible écart type. Pour générer un nombre quelconque de réalisations de T, on va à nouveau utiliser la fonction replicate. Il faut d’abord créer une fonction qui regroupe le calcul de la réalisation et retourne le résulat. random.T <- function(n) { # Générer le vecteur x de taille n # Calculer une réalisation de la loi T return(valeur) } Pour avoir une réalisation de la loi T, il suﬃt maintenant d’appeler la fonction random.T comme ceci random.T(n) Ensuite, pour appeler random.T un nombre quelconque de fois et stocker les diﬀérents résulats dans un vecteur, on exécute t.10 <- replicate(10, random.T(n)) qui stocke dans le vecteur t.10 10 appels successifs de la fonction random.T avec l’argument n. 10 Déﬁnir la fonction random.T et générer un vecteur t.1000 contenant 1000 réalisations de la variable aléatoire T. Vériﬁer que la moyenne empirique et la variance empirique sont en accord avec une loi normale centrée réduite. random.T <- function(n) { x <- runknown(n) (mean(x) - mu)/(sigma/sqrt(n)) } t.1000 <- replicate(1000, random.T(n)) mean(t.1000) [1] -0.01044766 var(t.1000) [1] 1.045376 La moyenne empirique semble nulle et la variance empirique égale à 1. On observe bien des estimations compatibles avec une loi normale centrée réduite. 11 Tracer la fonction de répartition empirique de l’échantillon t.1000 ; comparer la avec celle d’un échantillon x vu à la question 7. (6/14) TP 3 – SY02 Printemps 2022 t.1000 <- replicate(1000, random.T(n)) plot(ecdf(t.1000)) −4 −2 0 2 0.0 0.2 0.4 0.6 0.8 1.0 ecdf(t.1000) x Fn(x) La fonction de répartition de T pour n = 1000 est lissée et ressemble à celle d’une gaussienne centrée réduite. 12 Avec l’instruction curve(pnorm, add=TRUE), superposer au graphe précédent, la fonc- tion de répartition théorique d’une gaussienne centrée réduite. Qu’observez-vous ? t.1000 <- replicate(1000, random.T(n)) plot(ecdf(t.1000)) curve(pnorm, add = TRUE) −3 −2 −1 0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 1.0 ecdf(t.1000) x Fn(x) La variable aléatoire T est distribuée comme une variable normale centrée réduite. 13 En reprenant la question précédente, illustrer graphiquement le théorème central limite. (7/14) TP 3 – SY02 Printemps 2022 Ts <- replicate(1000, random.T(1)) plot(ecdf(Ts), main = "n=1") curve(pnorm, add = TRUE) Ts <- replicate(1000, random.T(2)) plot(ecdf(Ts), main = "n=2") curve(pnorm, add = TRUE) Ts <- replicate(1000, random.T(5)) plot(ecdf(Ts), main = "n=5") curve(pnorm, add = TRUE) Ts <- replicate(1000, random.T(10)) plot(ecdf(Ts), main = "n=10") curve(pnorm, add = TRUE) Ts <- replicate(1000, random.T(20)) plot(ecdf(Ts), main = "n=20") curve(pnorm, add = TRUE) Ts <- replicate(1000, random.T(50)) plot(ecdf(Ts), main = "n=50") curve(pnorm, add = TRUE) −2 −1 0 1 0.0 0.2 0.4 0.6 0.8 1.0 n=1 x Fn(x) −3 −2 −1 0 1 0.0 0.2 0.4 0.6 0.8 1.0 n=2 x Fn(x) −4 −3 −2 −1 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 n=5 x Fn(x) −5 −4 −3 −2 −1 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 n=10 x Fn(x) −4 −3 −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 n=20 x Fn(x) −4 −3 −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 n=50 x Fn(x) On observe une convergence simple de la fonction de répartition empirique vers la fonction de répartition d’une loi uploads/Religion/ tp3-estimation-corrige.pdf