# Pour puovoir reproduire les resultats :
set.seed(1234)
library(tidyverse)
library(plotly)

Objectifs du TP

Données simulées

\[ y = \beta_{1} + \beta_{2}x + \epsilon \]

Avec \(x\), \(\epsilon\) indépendantes.

On fixe \(\beta_{1} = -1\), \(\beta_{2} = 3\) et \(\sigma^2 = 2\).

Définition des paramètres et de la fonction.

beta1 <- -1.0
beta2 <- 3.0
sigma2 <- 2.0

modele <- function(x, eps) beta1 + beta2*x + eps
modele <- Vectorize(modele)

Génération de données

  1. Simuler \(n = 100\) réalisations indépendantes \((x_{i} , y_{i} )\) du couple \((x, y)\) défini par (1).
epsilon_chap <- rnorm(100, mean = 0, sd = sigma2^0.5)
x_i <- runif(100)
y_i <- modele(x_i, epsilon_chap)

Estimation ponctuelle des paramètres

  1. À l’aide du seul échantillon \((x_{i} , y_{i} )\), estimer par moindres carrés les paramètres \(\beta_{1}\) et \(\beta_{2}\) du modèle (1), ainsi que la variance résiduelle \(\sigma^2\). Comparer les valeurs estimées aux valeurs exactes de ces paramètres.
beta2_chap <- cov(x_i, y_i) / var(x_i)
beta1_chap <- mean(y_i) - beta2_chap*mean(x_i)
#glue::glue("beta1 ~ {beta1_chap}, beta2 ~ {beta2_chap}")
lm1 <- lm(y_i ~ x_i)
summary(lm1)
## 
## Call:
## lm(formula = y_i ~ x_i)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8382 -1.1232 -0.1934  0.9606  3.9752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.7379     0.3014  -2.448 0.016147 *  
## x_i           2.0346     0.5322   3.823 0.000231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.404 on 98 degrees of freedom
## Multiple R-squared:  0.1298, Adjusted R-squared:  0.1209 
## F-statistic: 14.62 on 1 and 98 DF,  p-value: 0.0002315
var(lm1$residuals)
## [1] 1.952109
Paramètre Réel \(\theta\) Estimation \(\hat{\theta}\)
pente \(\beta_{2}\) 3 2.0346201
ordonnée à l’origine \(\beta_{1}\) -1 -0.7378618

Représentation visuelle des données et du modèle estimé

donnees <- tibble::tibble(y=y_i, x=x_i)
plot1 <- donnees %>%
  ggplot(mapping = aes(x=x, y=y)) +
    geom_point() +
      geom_abline(
        slope = lm1$coefficients[["x_i"]],
        intercept = lm1$coefficients[["(Intercept)"]],
        colour = "blue"
      )
ggplotly(plot1)

Simulation

Estimation sur plusieurs expériences indépendantes

  1. Répéter les questions précédentes en changeant la valeur de la graine du générateur de nombres aléatoires (utiliser la fonction set.seed). Que constate-t-on ?
  • Cette question sera répondue avec la suivante.
  1. À l’aide d’un script R, stockez les valeurs estimées \(\hat{\beta}_{1}\), \(\hat{\beta}_{1}\), \(\hat{\sigma}^2\) obtenues pour \(N = 200\) échantillons indépendants.
  • Estimez la moyenne, la médiane, la variance de \(\hat{\beta}_{1}\). Commentez.
  • Représentez la loi de \(\hat{\beta}_{1}\) sous forme d’histogramme. Commentez.
  • Reprenez les deux questions ci-dessus en remplaçant \(\hat{\beta}_{1}\) par \(\hat{\beta}_{2}\), puis \(\hat{\sigma}^2\).
# source(here::here("utils", "shorthand.R"))
n.echantillons <- 200
experiences <- data.frame(beta1=double(), beta2=double(), sigma2=double())
for(i in 1:n.echantillons){
  set.seed(i)
  epsilon_chap <- rnorm(100, mean = 0, sd = sigma2^0.5)
  x_i <- runif(100)
  y_i <- modele(x_i, epsilon_chap)
  lm1 <- lm(y_i ~ x_i)
  experiences <- rbind(
    experiences, 
    data.frame(
      beta1 = lm1$coefficients[["(Intercept)"]],
      beta2 = lm1$coefficients[["x_i"]],
      sigma2= var(lm1$residuals) #/ (length(x_i) - 2)
    )
  )
} 
experiences <- tibble::tibble(experiences)
Moyenne Médiane Variance
\(\hat{\beta}_{1}\) -0.9716082 -0.9676745 0.0863733
\(\hat{\beta}_{2}\) 2.9688219 2.9652754 0.2634283
\(\hat{\sigma}^2\) 1.9582711 1.9538695 0.0793315

Les distributions des estimations des paramètres parmi les simulations semblent de ne pas être biaisées, du fait que leurs moyennes et leurs médianes sont assez proches. Avec \(N = 200\) simulations, les estimations des paramètres s’approchent des vraies valeurs fixées avant de construire le modèle. En comparant les simulations à la première estimation ponctuelle, la possibilité d’une première estimation biaisée est mise en évidence.

Pour complémenter l’analyse, une estimation d’un intervalle de confiance peut être faite. Afin d’avoir la certitude que l’hypothèse de normalité n’est pas contredite par les données, la fonction shapiro.test() est utilisée. Puis, si les p-valeurs ne se trouvent dans la zone de rejet, la fonction t.test() est utilisée pour calculer un \(IC95\%\).

## 
##  Shapiro-Wilk normality test
## 
## data:  experiences$beta1
## W = 0.99459, p-value = 0.6888
## 
##  Shapiro-Wilk normality test
## 
## data:  experiences$beta2
## W = 0.99532, p-value = 0.7961
## 
##  Shapiro-Wilk normality test
## 
## data:  experiences$sigma2
## W = 0.99065, p-value = 0.2216

Les p-valeurs ne permettent pas le rejet de \(H_{0}\) au seuil de \(\alpha = 0.05\). En supposant que les paramètres suivent une loi normale, les IC sont calculés :

## 
##  One Sample t-test
## 
## data:  experiences$beta1
## t = -46.754, df = 199, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0125883 -0.9306282
## sample estimates:
##  mean of x 
## -0.9716082
## 
##  One Sample t-test
## 
## data:  experiences$beta2
## t = 81.803, df = 199, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  2.897255 3.040389
## sample estimates:
## mean of x 
##  2.968822
## 
##  One Sample t-test
## 
## data:  experiences$sigma2
## t = 98.325, df = 199, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.918997 1.997545
## sample estimates:
## mean of x 
##  1.958271

Rapport moyenne-variance

\(abs(s^2(\theta) / \bar{x}(\theta))\)
\(\hat{\beta}_{1}\) 0.0888973
\(\hat{\beta}_{2}\) 0.0887316
\(\hat{\sigma}^2\) 0.040511

Les intervalles de confiance, avec le rapport \(abs(s^2(\theta) / \bar{x}(\theta))\) indiquent que avec plusieurs échantillons, il est possible d’obtenir une bonne estimation des valeurs des paramètres.

Les histogrammes correspondent bien à l’analyse fait précédemment. Seul le caractère apparemment bimodale de la distribution de \(\hat{\sigma}^2\) saut à l’œil comme un point important à discuter.

# Le code est le même pour les autres histogrammes 
source(here::here("utils", "plotting-funcs.R"))

plot.beta1 <- experiences %>%
  ggplot(mapping = aes(x = beta1)) +
    geom_histogram(
        aes(y = stat(density)), colour="black", fill="white", 
        binwidth = optim.binwidth(experiences$beta1)
    ) + 
        geom_density(alpha=.2, fill="blue") +
    ggtitle(glue::glue("B1 : {round(mean(experiences$beta1), 2)} vs {beta1}"))
ggplotly(plot.beta1)

\(x_{i}\) fixe, \(y_{i}\) recalculé

  1. Même question en générant une fois pour toutes les \(x_{i}\) (ici, seuls les \(y_{i}\) changent d’un échantillon à l’autre).
n.echantillons <- 200
experiences <- data.frame(beta1=double(), beta2=double(), sigma2=double())
# le vecteur x_i est crée hors de la boucle
x_i <- runif(100)
for(i in 1:n.echantillons){
  set.seed(i)
  epsilon_chap <- rnorm(100, mean = 0, sd = sigma2^0.5)
  y_i <- modele(x_i, epsilon_chap)
  lm1 <- lm(y_i ~ x_i)
  experiences <- rbind(
    experiences, 
    data.frame(
      beta1 = lm1$coefficients[["(Intercept)"]],
      beta2 = lm1$coefficients[["x_i"]],
      sigma2= var(lm1$residuals) #/ (length(x_i) - 2)
    )
  )
} 
experiences <- tibble::tibble(experiences)
Moyenne Médiane Variance
\(\hat{\beta}_{1}\) -0.9835802 -0.9851364 0.0920093
\(\hat{\beta}_{2}\) 2.9902003 2.9992879 0.2717644
\(\hat{\sigma}^2\) 1.961041 1.9564633 0.0790992