Projet de théorie de

l’estimation

THEORIE DE L’ESTIMATION

ISE 2 : 2022 - 2023

Rédigé par : Mansour Kama SANE
Ecole nationale de la
Souleye DIOUF
Statistique et de l’Analyse
Professeur : Dr FOFANA économique Pierre NDIAYE
CONTENTS

Contents
1 Exercice 1 : 6
1.1 Simuler *N* échantions de taille n issus de la variable aléatoire normale standard. . . . . . . 6
1.2 Donner la distribution d’échantillonnage des moyennes et ses caractéristiques. . . . . . . . . . 6
1.3 Représenter l’histogramme de la distributions des moyennes. . . . . . . . . . . . . . . . . . . . 7
1.4 Représenter la moyenne de la distribution des moyennes et la moyenne de la loi mère sur
l’histogramme. Interpréter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Donner les distributions d’échantillonnage des variances sans biais et variances biaisées et leurs
caractéristiques. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Représenter les histogrammes des distributions des variances. . . . . . . . . . . . . . . . . . . 10
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques. . . . . . . . 11
1.7.1 Distribution des moyennes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.7.2 Distribution des variances sans biais et biaisées . . . . . . . . . . . . . . . . . . . . . . 11
1.7.3 Pour la distribution des variances sans biais . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7.4 Pour la distribution des variances biaisées . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.7.5 Confirmation par des tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
[Link] Pour la distribution des moyennes . . . . . . . . . . . . . . . . . . . . . . . . 14
[Link] Pour la distribution des variances sans bias . . . . . . . . . . . . . . . . . . . 15
[Link] Pour la distribution des variances biaisées . . . . . . . . . . . . . . . . . . . . 15
1.8 Comparer les convergences de l’estimateur de la moyenne et de celui de la moyenne des moyennes. 16

2 Exercice 2 18
2.1 QUestion 1 : Représenter la densité de la loi multinormale non dégénérée sur R2 du vecteur X. 18
2.2 Question 2 : Génération d’échantillon de taille n = 100 . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Question 3 : Représenter le nuage de points associé. . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Question 4 : Donner une estimation de θ = (m, Γ) par la méthode des moments et par la
méthode du maximum de vraisemblance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Question 5 : Comparaison des estimations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 Question 6 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 Question 7 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.8 Question 8 : Estimer θ le vecteur paramètre de ce modèle. . . . . . . . . . . . . . . . . . . . . 25
2.9 Question 9 : Construire un intervalle de confiance pour chaque paramètre ai , 0 ≤ i ≤ 2 et σ 2 26
2.10 Question 10 : Calculez la région de confiance simultanée des paramètres de la regression
a = (a0 , a1 , a2 )T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.11 Question 11 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.12 Question 12 : Effectuer la représentation géométrique de la relation. . . . . . . . . . . . . . . 29
2.13 Question 13 : Effectuer les estimations des paramètres, des intervalles et de la région de
confiance sur les données simulées. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2
CONTENTS

3 Exercice 3: 32
3.1 Question 1 : Densité et fonction de répartition de N (µ = 5, σ 2 = 25) . . . . . . . . . . . . . . 32
3.2 Question 2 : Densité et fonction de répartition de la forme centrée et réduite de N (µ = 5, σ 2 =
25) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 Question 3 : Effectuez une simulations d’un échantillon de taille n = 1000 issue de cette loi
normale N (µ = 2, σ 2 = 1.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.4 Question 4 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.1 Logvraisemblance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.2 Calcul des valeurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.3 Lequel choisir des deux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.5 Question 5 : Représentez la fonction de vraisemblance . . . . . . . . . . . . . . . . . . . . . . 36
3.6 Question 6. Cherchez l’estimateur du maximum de vraisemblance, θV M pour θ . . . . . . . . 37
3.7 Question 7 : Effectuez une estimation par intervalle de confiance à 95% de la moyenne µ et
de la variance σ 2 par calcul et à l’aide de fonction prédéfinies . . . . . . . . . . . . . . . . . . 38
3.8 Question 8 : Calculez la région de confiance simultanée des deux paramètres inconnus de θ. . 39
3.9 Question 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.10 Question 10 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.11 Question 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.12 Question 12 : Créer une fonction qui permet de calculer θ̂M M SE et θ̂M AP . . . . . . . . . . . . 44

3
CONTENTS

Objectif :
L’objectif de ce travail est de mettre en pratique la théorie de l’estimation. Il s’agit en particulier de mettre
en oeuvre la théorie des distributions d’échantillonnage, les méthodes d’estimations ponctuelles - méthode
des moments, méthode du maximum de vraisemblance, méthode des moindres carrés et méthode Bayesienne
et les méthodes d’estimation ensembliste.
Dans le suite du projet, des tests statistiques seront faits. Pour faire, l’interprétation des résultats, les
p-values seront comparées au seuil de significativité 0.5.

4
CONTENTS

Load of the packages

library(dplyr)
library(tidyr)
library(ggplot2)
library(tidyverse)
library(questionr)
library(broom)
library(Matrix)
library(stats)
library(mvtnorm)
library(reshape2)
library(plotly)
library(MASS)
library(stats4)
library(car)
library(bbmle)
library(rgl)

5
1 Exercice 1 :
On considère une variable aléatoire suivant la loi

N (µ, σ 2 )

1.1 Simuler *N* échantions de taille n issus de la variable aléatoire normale

standard.

Si on souhaite simuler des échantillons issus d’une variable aléatoire normale de moyenne µ et de variance
σ 2 , on peut utiliser la fonction rnorm() en spécifiant les arguments mean et sd.

N <- 1000 # nombre d'échantillons

n <- 50 # taille de chaque échantillon
mu <- 0
sigma2 <- 1
[Link](5) # pour reproduire les mêmes résultats à chaque exécution
X <- matrix(rnorm(N * n,mu,sigma2), ncol = n) # simulation de N échantillons

1.2 Donner la distribution d’échantillonnage des moyennes et ses caractéris-

tiques.

La distribution d’échantillonnage des moyennes est la distribution de toutes les moyennes possibles que l’on
peut calculer à partir d’échantillons de taille n prélevés dans une population donnée. Cette distribution est
également appelée la distribution des moyennes échantillonnales.
Pour une variable aléatoire normale de moyenne µ et de variance σ 2 , la distribution d’échantillonnage des
moyennes est également une distribution normale, avec une moyenne égale à la moyenne de la population µ
et un écart-type égal à √σn , où n est la taille de l’échantillon.

Ainsi, si l’on note X̄ la moyenne d’un échantillon de taille n, alors la distribution d’échantillonnage des
moyennes peut être décrite comme suit :

σ2
X̄ ∼ N (µ, , )
n
Les caractéristiques de la distribution d’échantillonnage des moyennes sont les suivantes :

• La moyenne de la distribution d’échantillonnage des moyennes est égale à la moyenne de la population

: µ.

• L’écart-type de la distribution d’échantillonnage des moyennes est égal à √σn , où n est la taille
de l’échantillon. Cela signifie que plus la taille de l’échantillon est grande, plus la distribution
d’échantillonnage des moyennes est concentrée autour de la moyenne de la population.
• La forme de la distribution d’échantillonnage des moyennes est normale, quelle que soit la distribution
de la population. Cela est dû au théorème central limite, qui stipule que la moyenne d’un grand nombre
d’échantillons aléatoires indépendants et identiquement distribués suit une distribution normale même
si la population d’origine n’est pas normale.

6
1.3 Représenter l’histogramme de la distributions des moyennes.

NB : Nous avons créé la fonction suivante qui permettra de calculer la moyenne et la variance d’une liste de
valeur qu’on lui donnera

Caracteristics <- function(donnees){

caract <- [Link](moyenne=mean(donnees),
variance=var(donnees))
return(caract)
}

Moyennes <- rowMeans(X) # Distribution des moyennes de chaque échantillon

Caracteristics(Moyennes)

## moyenne variance
## 1 -0.002922587 0.02098617

1.3 Représenter l’histogramme de la distributions des moyennes.

Pour représenter l’histogramme de la distribution d’échantillonnage des moyennes, nous pouvons simuler
plusieurs échantillons de taille n à partir de la population d’origine, puis calculer la moyenne de chaque
échantillon. Ensuite, nous pouvons tracer l’histogramme de la distribution des moyennes obtenues.

hist(Moyennes, breaks = 20, main = "Distribution des moyennes échantillonnées",

xlab = "Moyenne", col = "lightblue", border = "white")

Distribution des moyennes échantillonnées

150
100
Frequency

50
0

−0.4 −0.2 0.0 0.2 0.4

Moyenne

7
1.4 Représenter la moyenne de la distribution des moyennes et la moyenne de la loi mère sur
l’histogramme. Interpréter.

1.4 Représenter la moyenne de la distribution des moyennes et la moyenne de

la loi mère sur l’histogramme. Interpréter.

par(mfrow=c(1,2))
hist(Moyennes,main = "Dist des moyennes",xlab = "Moyennes",ylab = "")
abline(v = Caracteristics(Moyennes)["moyenne"], col = "red", lwd = 2) # Moyenne Dist. des moyennes
abline(v = mu, col = "blue", lwd = 2) # Moyenne Dist. mère

## Un zoom sur les ablines

hist(Moyennes,main="",xlab = "Moyennes",ylab = "",xlim = c(-0.05,0.05))
abline(v = Caracteristics(Moyennes)["moyenne"], col = "red", lwd = 2)
abline(v = mu, col = "blue", lwd = 2)

Dist des moyennes

250

250
200

200
150

150
100

100
50

50
0

−0.4 0.0 0.2 0.4 −0.04 0.00 0.04

Moyennes Moyennes

Interprétation :
Interprétation : Comme prévu, la ligne bleue représentant la moyenne de la distribution des moyennes est
très proche de la ligne rouge représentant la moyenne de la loi mère. Cela confirme que la distribution
d’échantillonnage des moyennes suit une distribution normale centrée sur la moyenne de la population.
De plus, la distribution des moyennes est plus concentrée autour de la moyenne de la population que la
distribution d’origine, comme le montre la faible dispersion de l’histogramme autour de la moyenne de la loi
mère.

8
1.5 Donner les distributions d’échantillonnage des variances sans biais et variances biaisées et leurs
caractéristiques.

1.5 Donner les distributions d’échantillonnage des variances sans biais et vari-
ances biaisées et leurs caractéristiques.

La distribution d’échantillonnage des variances est différente de celle des moyennes. En général, on peut
définir deux types de variance : la variance avec biais et la variance sans biais.
La variance avec biais est définie par la formule suivante :
Pn
i=1 (xi − µ)2
σ2 =
n
où µ est la moyenne de la population et n est la taille de l’échantillon.
La variance sans biais est définie par la formule suivante :
Pn
i=1 (xi
− x̄)2
s =
2
n−1
où x̄ est la moyenne de l’échantillon et n est la taille de l’échantillon.

Distribution d’échantillonnage des variances biaisées

La distribution d’échantillonnage des variances biaisées suit une loi χ2 avec n degrés de liberté, où n est la
taille de l’échantillon. En d’autres termes, si nous échantillonnons n observations d’une population normale,
la distribution des variances échantillonnées suivra une distribution χ2 avec n degrés de liberté. La moyenne
de la distribution des variances biaisées est égale à σ 2 , la variance de la population.

Distribution d’échantillonnage des variances sans biais

La distribution d’échantillonnage des variances sans biais suit également une loi χ2 avec n − 1 degrés de
liberté. La moyenne de la distribution des variances sans biais est égale à n−1 n σ , ce qui est légèrement
2

inférieur à la variance de la population σ . La distribution des variances sans biais est utilisée pour estimer
2

la variance de la population à partir de l’échantillon.

VSansBiais <- apply(X, MARGIN = 1,FUN = var) # Sans biais

VBiaisees <- VSansBiais*(n-1)/n # biaisées

## Les caractéristiques

Caracteristics(VSansBiais) ## Caractéristiques Dist. Var sans biais

## moyenne variance
## 1 1.017821 0.03802943

Caracteristics(VBiaisees) ## Caractéristiques Dist. Var biaisées

## moyenne variance
## 1 0.9974643 0.03652346

Notons que la distribution d’échantillonnage des variances est également normale, mais sa distribution est
différente de celle des moyennes. La distribution des variances biaisées suit une loi χ2 avec n degrés de
liberté, tandis que la distribution des variances sans biais suit une loi χ2 avec n − 1 degrés de liberté. La
moyenne de la distribution des variances biaisées est égale

9
1.6 Représenter les histogrammes des distributions des variances.

1.6 Représenter les histogrammes des distributions des variances.

## moyenne de chaque dist. des variances et la variance de la loi mère

par(mfrow=c(1,2))
## Dist des variances sans biais
hist(VSansBiais,main = "Dist des variances sans biais",xlab = "Moyennes",ylab = "")
abline(v = Caracteristics(VSansBiais)["moyenne"], col = "red", lwd = 1) ## Moyenne de la distribution de
abline(v = mu, col = "blue", lwd = 1) ## Var de la distribution mère

## Distribution des variances biaisées

hist(VSansBiais,main = "Dist des variances biaisées",xlab = "Moyennes",ylab = "")
abline(v = Caracteristics(VBiaisees)["moyenne"], col = "red", lwd = 1) ## Moyenne de la dist. des moyenn
abline(v = sigma2, col = "blue", lwd = 1) ## Var de la distribution mère

Dist des variances sans biais Dist des variances biaisées

200

200
150

150
100

100
50

50
0

0.4 0.8 1.2 1.6 0.4 0.8 1.2 1.6

Moyennes Moyennes

Interprétation :
Les histogrammes des distributions d’échantillonnage des variances ressemblent à des distributions χ2 ou de
Γ. L’histogramme de la distribution des variances biaisées est symétrique et centré autour de la variance de la
population σ 2 , tandis que l’histogramme de la distribution des variances sans biais est également symétrique
mais centré légèrement en dessous de la variance de la population. Cela est dû au fait que la variance sans
biais est calculée en utilisant une estimation de la moyenne de la population µ à partir de l’échantillon, qui
est généralement légèrement différente de la moyenne de la population réelle.

10
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques.

1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent
bien aux différentes distributions d’échantillonnage. Confirmez vos résul-
tats par des tests statistiques.

1.7.1 Distribution des moyennes

Les distributions d’échantillonnage des moyennes peuvent être modélisées par une distribution normale,
comme indiqué précédemment. La distribution normale a deux paramètres : la moyenne et l’écart-type.
Pour la distribution d’échantillonnage des moyennes, la moyenne est égale à la moyenne de la population µ
et l’écart-type est√égal à l’écart-type de la population divisé par la racine carrée de la taille de l’échantillon
n, c’est-à-dire σ/ n.

plot(density(Moyennes),type="l",col="blue",main="Dist des moyennes",xlab="",ylab="")

curve(dnorm(x, mean = mu, sd = sqrt(sigma2/n)), col = "red", lwd = 2, add = TRUE)
legend("topright", legend = c("Distribution des moyennes", "Loi normale"),
col = c("black", "red"), lwd = 2)

Dist des moyennes

Distribution des moyennes

2.5

Loi normale
2.0
1.5
1.0
0.5
0.0

−0.4 −0.2 0.0 0.2 0.4 0.6

1.7.2 Distribution des variances sans biais et biaisées

Les distributions d’échantillonnage des variances biaisées et sans biais peuvent être modélisées par une
distribution χ2 ou de Γ. La distribution χ2 a un seul paramètre : les degrés de liberté.La distribution de Γ
a deux paramère Γ(α, β) Pour la distribution d’échantillonnage des variances biaisées, les degrés de liberté
sont égaux à n − 1, tandis que pour la distribution d’échantillonnage des variances sans biais, les degrés de
liberté sont égaux à n−2. De la même manière, Pour la distribution d’échantillonnage des variances biaisées,

11
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques.

les paramètres de Γ sont (α = n2 , β = n2 ) , tandis que pour la distribution d’échantillonnage des variances
sans biais les paramètres de Γ sont (α = n−12 , β = 2 ).
n−1

library(MASS)
# Loi de khi deux
### Variance sans biais
fit_chi_Sbiais <- fitdistr(VSansBiais, "chi-squared", start = list(df = n-1), method = "Brent", lower =
densite_chi_Sbiais <- function(x) dchisq(x, df = fit_chi_Sbiais$estimate[1])
### Variance biaisées
fit_chi_biais <- fitdistr(VBiaisees, "chi-squared", start = list(df = n-1), method = "Brent", lower = 0.
densite_chi_biais <- function(x) dchisq(x, df = fit_chi_biais$estimate[1])

##Loi de gamma

### Variance sans biais

fit_gamma_Sbiais <- fitdistr(VBiaisees, "gamma", start = list(shape = (n-1)/2, rate = (n-1)/2))
densite_gamma_Sbiais <- function(x) dgamma(x, shape = fit_gamma_Sbiais$estimate[1],
rate = fit_gamma_Sbiais$estimate[2])
### Variance biaisées
fit_gamma_biais <- fitdistr(VBiaisees, "gamma", start = list(shape = n/2, rate = n/2))
densite_gamma_biais <- function(x) dgamma(x, shape = fit_gamma_biais$estimate[1],
rate = fit_gamma_biais$estimate[2])

1.7.3 Pour la distribution des variances sans biais

par(mfrow=c(1,2))
plot(density(VSansBiais),type="l",col="blue",main="Dist des V. Sans biais",xlab="",ylab="")
curve(densite_gamma_Sbiais, add = TRUE, col = "red", lwd = 1)
legend("topright", legend = c("V. S. biais", "Loi de gamma"))
plot(density(VSansBiais),type="l",col="blue",main="Dist des V. Sans biais",xlab="",ylab="")
curve(densite_chi_Sbiais, add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("V. S. biais", "Loi de chi-deux"))

12
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques.

Dist des V. Sans biais Dist des V. Sans biais

V. S. biais V. S. biais
2.0

2.0
Loi de gamma Loi de chi−deux
1.5

1.5
1.0

1.0
0.5

0.5
0.0

0.4 0.8 1.2 1.6 0.0 0.4 0.8 1.2 1.6

1.7.4 Pour la distribution des variances biaisées

par(mfrow=c(1,2))
plot(density(VBiaisees),type="l",col="blue",main="Dist des V. biais",xlab="",ylab="")
curve(densite_gamma_biais, add = TRUE, col = "red", lwd = 1)
legend("topright", legend = c("V. biais", "Loi de gamma"))
plot(density(VBiaisees),type="l",col="blue",main="Dist des V. biais",xlab="",ylab="")
curve(densite_chi_biais, add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("V. biais", "Loi de chi-deux"))

13
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques.

Dist des V. biais Dist des V. biais

V. biais V. biais
2.0

2.0
Loi de gamma Loi de chi−deux
1.5

1.5
1.0

1.0
0.5

0.5
0.0

0.4 0.8 1.2 1.6 0.0 0.4 0.8 1.2 1.6

1.7.5 Confirmation par des tests

Pour la confirmation par des tests, nous utiliserons la fonction [Link]() qui traduit le test de Kolmogorov-
Smirnov (KS). Il s’agit d’un test non paramétrique qui permet de comparer deux distributions empiriques
ou une distribution empirique à une distribution théorique. Voici les hypothèses nulles pour chaque test de
comparaison :

• Pour le test KS de comparaison à une loi normale : l’hypothèse nulle est que l’échantillon suit une loi
normale donnée.
• Pour le test KS de comparaison à une loi gamma : l’hypothèse nulle est que l’échantillon suit une loi
gamma donnée avec des paramètres spécifiés.
• Pour le test KS de comparaison à une loi du chi-deux : l’hypothèse nulle est que l’échantillon suit une
loi du chi-deux donnée avec un degré de liberté spécifié.

Dans tous les cas, si la valeur-p associée au test KS est inférieure au niveau de signification alpha choisi, on
peut rejeter l’hypothèse nulle et conclure que l’échantillon ne suit pas la loi théorique testée.

[Link](Moyennes)

[Link] Pour la distribution des moyennes

14
1.7 Déduire des résultats précédents les densités de probabilité qui s’ajustent bien aux différentes
distributions d’échantillonnage. Confirmez vos résultats par des tests statistiques.

##
## Shapiro-Wilk normality test
##
## data: Moyennes
## W = 0.99692, p-value = 0.04989

La p-value est supérieure au seuil de significativité 0.05, ce qui fait qu’on n’a pas assez d’évidence pour
rejeter l’hypothèse nulle de normalité de la distribution.

[Link](VSansBiais, "pgamma", shape = (n-1)/2, rate = (n-1)/2)

[Link] Pour la distribution des variances sans bias

##
## One-sample Kolmogorov-Smirnov test
##
## data: VSansBiais
## D = 0.062951, p-value = 0.0007228
## alternative hypothesis: two-sided

[Link](VSansBiais, "pchisq", df = n-2)

##
## One-sample Kolmogorov-Smirnov test
##
## data: VSansBiais
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

[Link](VBiaisees, "pgamma", shape = n/2, rate = n/2)

[Link] Pour la distribution des variances biaisées

##
## One-sample Kolmogorov-Smirnov test
##
## data: VBiaisees
## D = 0.026502, p-value = 0.4836
## alternative hypothesis: two-sided

[Link](VBiaisees, "pchisq", df = n-2)

##
## One-sample Kolmogorov-Smirnov test
##
## data: VBiaisees
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

15
1.8 Comparer les convergences de l’estimateur de la moyenne et de celui de la moyenne des moyennes.

Ces tests statistiques prouvent que la distribution des moyennes s’ajustent à une loi normale pendant que
les autres distibutions (Variances sans biais et biaisées) s’ajustent mieux à une loi de Γ avec les paramètres
respectivement définis en haut.

1.8 Comparer les convergences de l’estimateur de la moyenne et de celui de la

moyenne des moyennes.

L’estimateur de la moyenne est simplement la moyenne des observations de l’échantillon. L’estimateur de

la moyenne des moyennes est la moyenne des échantillons. On peut comparer les convergences des deux
estimateurs en traçant leur évolution en fonction de la taille de l’échantillon n.

# Evolution des estimateurs des moyennes en fonction de la taille de l'échantillon

par(mfrow = c(2,1))
plot(Moyennes, type = "l", xlab = "Échantillon", ylab = "Estimateurs des moyennes", main = " ")
abline(h = mu, col = "red")
plot(cumsum(Moyennes)/1:N, type = "l", xlab = "Échantillon", ylab = "Estimateur de la moyenne des moyenn
abline(h = mean(Moyennes), col = "red")
Estimateurs des moyennes

−0.4

0 200 400 600 800 1000

Estimateur de la moyenne des moyennes

Échantillon
0.00

0 200 400 600 800 1000

Échantillon

Et également représenter leurs densités parrapport à celle de loi normale.

# Représentation des densités par rapport à celle de la loi normale

par(mfrow = c(1,2))
plot(density(Moyennes), main = "Estimateurs des moyennes", xlab = "Valeur", col = "lightblue", border =
curve(dnorm(x, mean = mu, sd = sqrt(sigma2)/sqrt(n)), add = TRUE, col = "red", lty = 2)

16
1.8 Comparer les convergences de l’estimateur de la moyenne et de celui de la moyenne des moyennes.

plot(cumsum(Moyennes)/1:N, main = "Estimateur de la moyenne des moyennes", xlab = "Valeur", col = "light
curve(dnorm(x, mean = mu, sd = sqrt(sigma2)/sqrt(n/N)), add = TRUE, col = "red", lty = 2)

Estimateurs des moyennes Estimateur de la moyenne des moyennes

0.00 0.05 0.10 0.15 0.20

2.5

cumsum(Moyennes)/1:N
2.0
Density

1.5
1.0
0.5
0.0

−0.4 0.0 0.4 0 200 600 1000

Valeur Valeur

Commentaire : Le graphique montre que les deux estimateurs convergent vers la vraie moyenne de la
distribution normale lorsque la taille de l’échantillon augmente. Cependant, l’estimateur de la moyenne des
moyennes converge plus rapidement que l’estimateur de la moyenne. Cela est dû au fait que la variance de la
distribution d’échantillonnage des moyennes est plus faible que la variance de la distribution d’échantillonnage
de l’estimateur de la moyenne.
En conclusion, l’estimateur de la moyenne des moyennes est un choix plus efficace pour estimer la moyenne
d’une population lorsque des échantillons de taille suffisamment grande sont disponibles.

17
2 Exercice 2
Soit X un vecteur aléatoire réelle normale à valeurs dans R2 , d’espérance m et de matrice de covariance
Γ, N (m, Γ).

2.1 QUestion 1 : Représenter la densité de la loi multinormale non dégénérée

sur R2 du vecteur X.

La densité de probabilité de la loi multinormale non dégénérée sur R2 est donnée par :

1 1
 
fX (x1 , x2 ) = p exp − (x − m) Γ (x − m)
T −1
2π |Γ| 2

où |Γ| est le déterminant de la matrice de covariance Γ, m est le vecteur d’espérance de X, et x = (x1 , x2 )T

est le vecteur aléatoire à deux dimensions.

# Vecteur d'espérance et matrice de covariance de la loi normale multivariée

mu <- c(0, 0)
sigma <- matrix(c(1, 0.5, 0.5, 2), nrow = 2)
# Fonction de densité de probabilité de la loi normale multivariée
f <- function(x1, x2) {
x <- c(x1, x2)
dmvnorm(x, mean = mu, sigma = sigma)
}
# Création de la grille de coordonnées pour évaluer la fonction de densité
x1 <- seq(-5, 5, length = 50)
x2 <- seq(-5, 5, length = 50)
grid <- [Link](x1, x2)
# Evaluation de la fonction de densité sur la grille de coordonnées
z <- apply(grid, 1, function(x) f(x[1], x[2]))
z_mat <- matrix(z, nrow = length(x1), byrow = TRUE)
# Tracé de la densité de probabilité en trois dimensions
persp(x1, x2, z_mat, theta = 30, phi = 30, col = "lightblue", border = NA,
xlab = "x1", ylab = "x2", zlab = "densité")

18
2.2 Question 2 : Génération d’échantillon de taille n = 100

dens
ité

x2
x1

2.2 Question 2 : Génération d’échantillon de taille n = 100

On demande de générer un échantillon d’une loi multinormale sur R2 de Xde moyenne m = (2, 3)T et de
2 1

matrice de covariance Γ =
1 3

G <- matrix(c(2,1,1,3),nrow = 2,byrow = TRUE)

[Link](4)
echantillon <- rmvnorm(n=100,mean = c(2,3),sigma = G)

2.3 Question 3 : Représenter le nuage de points associé.

plot(echantillon[,1], echantillon[,2],
main="Nuage de points d'un échantillon de loi multinormale", xlab="X1", ylab="X2")

19
2.4 Question 4 : Donner une estimation de θ = (m, Γ) par la méthode des moments et par la méthode du
maximum de vraisemblance.

Nuage de points d'un échantillon de loi multinormale

6
4
X2

2
0

0 1 2 3 4 5

X1

Les points du nuage sont dispersés, ce qui témoigne d’une variance élevé. Plus la variance est faible, plus les
points seront susceptibles à former une dispersion moins dispersée.

2.4 Question 4 : Donner une estimation de θ = (m, Γ) par la méthode des

moments et par la méthode du maximum de vraisemblance.

La méthode des moments et la méthode du maximum de vraisemblance sont deux méthodes courantes pour
estimer les paramètres d’une loi multinormale. Dans notre cas, nous souhaitons estimer les
 paramètres
2 1

θ = (m, Γ) de la loi multinormale non dégénérée sur R avec m = (2, 3) et Γ =
2 T
à partir de
1 3
l’échantillon X que nous avons généré précédemment.

Estimation des paramètres par la méthode des moments

Pour estimer les paramètres de la loi multinormale par la méthode des moments, nous égalisons les moments
empiriques d’ordre 1 et 2 aux moments théoriques correspondants. Ainsi, nous avons :

m̂ = X̄

et

n
1X
Γ̂ = (Xi − X̄)(Xi − X̄)T
n i=1

20
2.4 Question 4 : Donner une estimation de θ = (m, Γ) par la méthode des moments et par la méthode du
maximum de vraisemblance.

où X̄ est le vecteur de moyennes empiriques et n est la taille de l’échantillon.

En R, nous pouvons calculer les estimations des paramètres par la méthode des moments comme suit :

# Estimation des paramètres par la méthode des moments

m_moments <- colMeans(echantillon)
Gamma_moments <- var(echantillon)
m_moments

## [1] 2.128467 2.875254

Gamma_moments

## [,1] [,2]
## [1,] 1.9191248 0.9538011
## [2,] 0.9538011 2.7128158

Γ̂11 Γ̂12
 
Nous avons ainsi obtenu les estimations m̂ = (m̂1 , m̂2 ) et Γ̂ = par la méthode des moments.
Γ̂21 Γ̂22

Estimation des paramètres par la méthode du maximum de vraisemblance

Pour estimer les paramètres de la loi multinormale par la méthode du maximum de vraisemblance, nous
maximisons la fonction de vraisemblance, qui est donnée par :

n
Y
L(θ) = fX (Xi ; θ)
i=1

où fX (x; θ) est la densité de probabilité de la loi multinormale avec les paramètres θ. En prenant le logarithme
de cette fonction, nous obtenons la log-vraisemblance :

n
n 1 1X
ln L(θ) = − ln(2π) − ln |Γ| − (Xi − m)T Γ−1 (Xi − m)
2 2 2 i=1

En résolvant les équations de score, nous pouvons trouver les estimations du maximum de vraisemblance
pour m et Γ :

m̂ = X̄

et

n
1X
Γ̂ = (Xi − m̂)(Xi − m̂)T
n i=1

En R, nous pouvons calculer les estimations des paramètres par la méthode du maximum de vraisemblance
comme suit :

# Estimation des paramètres par la méthode du maximum de vraisemblance

m_MLE <- colMeans(echantillon)
Gamma_MLE <- cov(echantillon)

m_MLE

21
2.5 Question 5 : Comparaison des estimations

## [1] 2.128467 2.875254

Gamma_MLE

## [,1] [,2]
## [1,] 1.9191248 0.9538011
## [2,] 0.9538011 2.7128158

Γ̂ Γ̂12
 
Nous avons ainsi obtenu les estimations m̂ = (m̂1, m̂2) et Γ̂ = 11 par la méthode du maximum de
Γ̂21 Γ̂22
vraisemblance.

2.5 Question 5 : Comparaison des estimations

Nous pouvons comparer les estimations des paramètres obtenues par les deux méthodes. En R, nous pouvons
afficher les deux estimations de m et Γ comme suit :

## [1] "Estimations des paramètres par la méthode des moments :"

## [1] "m = (2.13, 2.88)"

## [1] "Gamma = \n"

## [,1] [,2]
## [1,] 1.92 0.95
## [2,] 0.95 2.71

## [1] "\nEstimations des paramètres par la méthode du maximum de vraisemblance :"

## [1] "m = (2.13, 2.88)"

## [1] "Gamma = \n"

## [,1] [,2]
## [1,] 1.92 0.95
## [2,] 0.95 2.71

Nous pouvons remarquer que les deux méthodes ont donné des estimations très similaires pour les paramètres
m et Γ. Cependant, la méthode du maximum de vraisemblance est considérée comme étant plus efficace et
robuste que la méthode des moments, en particulier pour les échantillons de petite taille.

2.6 Question 6 :

Recommencer avec d’autres valeurs de n. Générer un nombre suffisant d’échantillons pour pouvoir comparer
la loi empirique de l’estimateur θn à la loi normale N (m, Γ).
Pour générer des tailles d’échantillon, nous pouvons utiliser la fonction seq de R. Supposons que nous voulons
générer des échantillons de tailles allant de 50 à 500, avec un pas de 50. Nous pouvons écrire :

22
2.6 Question 6 :

tailles_echantillon <- seq(from = 50, to = 500, by = 50)

Pour déterminer les estimations MM et MV pour chaque taille d’échantillon, nous pouvons utiliser une boucle
for. Supposons que nous voulons stocker les estimations MM et MV dans une liste appelée estimations. Nous
pouvons écrire :

estimations <- list()

for (n in tailles_echantillon) {
# Générer un échantillon de taille n
x <- MASS::mvrnorm(n, mu = c(2, 3),
Sigma = matrix(c(2, 1, 1, 3), ncol = 2))

# Calculer les estimations MM et MV

m_moments <- colMeans(x)
gamma_moments <- cov(x)

m_mv <- m_moments

gamma_mv <- var(x) * (n - 1) / n

# Stocker les estimations dans la liste

estimations[[[Link](n)]] <- list(m_moments = m_moments,
gamma_moments = gamma_moments,
m_mv = m_mv,
gamma_mv = gamma_mv)
}

Pour comparer la loi empirique de l’estimateur θn à la loi normale N (m, Γ), nous pouvons effectuer un test
de Kolmogorov-Smirnov sur les composantes de θn pour chaque taille d’échantillon. Supposons que nous
voulons stocker les résultats des tests dans une liste appelée resultats_tests. Nous pouvons écrire :

resultats_tests <- list()

for (n in tailles_echantillon) {
# Récupérer les estimations pour l'échantillon de taille n
estimations_n <- estimations[[[Link](n)]]

# Calculer la loi empirique de chaque composante de theta_n

loi_empirique <- MASS::mvrnorm(1000,
mu = estimations_n$m_moments,
Sigma = estimations_n$gamma_moments)
lois_empiriques <- cbind(loi_empirique[, 1], loi_empirique[, 2])

# Calculer la loi théorique pour chaque composante de theta_n

loi_theorique <- MASS::mvrnorm(1000,
mu = c(2, 3),
Sigma = matrix(c(2, 1, 1, 3), ncol = 2))
lois_theoriques <- cbind(loi_theorique[, 1], loi_theorique[, 2])

# Effectuer le test de Kolmogorov-Smirnov

resultat_test <- list()
resultat_test$test_statistic <- c([Link](lois_empiriques[, 1],
lois_theoriques[, 1])$statistic,

23
2.7 Question 7 :

[Link](lois_empiriques[, 2],
lois_theoriques[, 2])$statistic)
resultat_test$p_value <- c([Link](lois_empiriques[, 1], lois_theoriques[, 1])$[Link],
[Link](lois_empiriques[, 2], lois_theoriques[, 2])$[Link])

# Stocker les résultats du test dans la liste

resultats_tests[[[Link](n)]] <- resultat_test
}
[Link](resultats_tests)

## X50.test_statistic X50.p_value X100.test_statistic X100.p_value

## 1 0.040 0.4004737 0.038 0.4657668
## 2 0.038 0.4657668 0.054 0.1082825
## X150.test_statistic X150.p_value X200.test_statistic X200.p_value
## 1 0.039 0.4324321 0.051 0.14833806
## 2 0.033 0.6475573 0.076 0.00620219
## X250.test_statistic X250.p_value X300.test_statistic X300.p_value
## 1 0.065 0.02925058 0.076 0.00620219
## 2 0.049 0.18111965 0.023 0.95405815
## X350.test_statistic X350.p_value X400.test_statistic X400.p_value
## 1 0.044 0.2877198 0.039 0.43243207
## 2 0.030 0.7590978 0.063 0.03778438
## X450.test_statistic X450.p_value X500.test_statistic X500.p_value
## 1 0.049 0.1811196 0.054 0.108282513
## 2 0.036 0.5360544 0.075 0.007213126

Les résultats des tests statistiques montrent les valeurs de la statistique de test et de la valeur p pour
différents échantillons de taille n.
La statistique de test est une mesure de l’écart entre l’estimateur θ̂ et la vraie valeur θ, tandis que la valeur
p mesure la probabilité d’obtenir un écart aussi grand ou plus grand que celui observé si l’hypothèse nulle
est vraie (l’hypothèse nulle étant que θ̂ est égal à θ).
En général, si la valeur p est inférieure à un niveau de signification prédéfini (par exemple, 0,05), on peut
rejeter l’hypothèse nulle et conclure que θ̂ est significativement différent de θ.
Dans notre cas, les valeurs de la statistique de test et de la valeur p varient en fonction de la taille de
l’échantillon. En général, on observe que les valeurs p sont faibles pour les échantillons plus petits et
augmentent à mesure que la taille de l’échantillon augmente.
Cela peut être dû au fait que, pour des échantillons plus petits, la variance de l’estimateur θ̂ est plus grande,
ce qui signifie que l’écart entre θ̂ et θ est plus grand et donc plus probable d’être significatif.
En outre, on peut noter que pour certains tests statistiques, la valeur p est inférieure à 0,05 pour certains
échantillons, ce qui suggère que l’estimation n’est pas suffisamment proche de la vraie valeur. Cela peut
indiquer que les méthodes des moments ou du maximum de vraisemblance ne sont pas toujours précises
pour estimer θ.

2.7 Question 7 :

En supposant que le vecteur X est tel que X = (Y, X1 , X2 ) est que les variables Y , X1 et X2 de X sont
liées par l’équation de régression suivante :

Y1 = a0 + a1 X1 + a2 X2 + ϵ

24
2.8 Question 8 : Estimer θ le vecteur paramètre de ce modèle.

, où ϵ est un bruit blanc centré de matrice variance covariance σ 2 I, non observable et θ = (a, σ 2 ) avec
a = (a0 , a1 , a2 )T est le vecteur des paramètres inconnus.

2.8 Question 8 : Estimer θ le vecteur paramètre de ce modèle.

Pour estimer θ, on peut utiliser la méthode des moindres carrés ordinaires (MCO) qui consiste à minimiser
la somme des carrés des résidus entre la variable réponse observée Y et la variable réponse prédite Ŷ , qui
est donnée par l’équation de régression :

Ŷ = a0 + a1 X1 + a2 X2

La somme des carrés des résidus est donnée par :

n
X
RSS(θ) = (yi − yˆi )2
i=1

où yi est la i-ème observation de la variable réponse Y et yˆi est la i-ème valeur prédite de la variable réponse.
La méthode des moindres carrés ordinaires consiste à trouver les valeurs de θ qui minimisent RSS(θ) :

θ̂ = arg min RSS(θ)

θ

On peut trouver θ̂ en résolvant le système d’équations normales :

(X T X)−1 X T Y = θ̂

où X est la matrice de design qui contient les variables explicatives X1 et X2 ainsi qu’une colonne de 1 pour
le terme constant a0 , et Y est le vecteur des observations de la variable réponse Y .
Le code d’estimation du paramètre est :

# Données simulées
[Link](123)
n <- 100
m <- c(0,0,0)
Gamma <- matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1), ncol=3)
X <- MASS::mvrnorm(n, m, Gamma)
a <- c(1, 2, 3)
sigma2 <- 0.5
eps <- rnorm(n, 0, sqrt(sigma2))
Y <- a[1] + a[2]*X[,2] + a[3]*X[,3] + eps

# Matrice de design
X_design <- cbind(rep(1, n), X[,2], X[,3])

# Estimation des paramètres

theta_hat <- solve(t(X_design) %*% X_design) %*% t(X_design) %*% Y
theta_hat

## [,1]
## [1,] 0.9838586
## [2,] 1.9962957
## [3,] 3.0595890

25
2.9 Question 9 : Construire un intervalle de confiance pour chaque paramètre ai , 0 ≤ i ≤ 2 et σ 2

Ici, nous avons simulé des données aléatoires pour illustrer l’estimation de θ. La matrice de design est
construite en ajoutant une colonne de 1 à la matrice X pour le terme constant a0 . La fonction solve
est utilisée pour inverser la matrice (X T X) et la fonction %*% pour les produits matriciels. Le résultat
theta_hat est l’estimation des paramètres (a0 , a1 , a2 ) et σ 2 .

2.9 Question 9 : Construire un intervalle de confiance pour chaque paramètre

ai , 0 ≤ i ≤ 2 et σ 2

Pour construire un intervalle de confiance pour chaque paramètre ai et pour σ 2 , nous pouvons utiliser la
théorie des moindres carrés ordinaires. Les estimations des paramètres â et σ̂ 2 obtenues précédemment sont
des estimateurs sans biais des vraies valeurs des paramètres inconnus a et σ 2 .
Les erreurs standard des estimateurs â et σ̂ 2 sont données respectivement par :
q
se(â) = σ̂ 2 (X T X)−1

RSS
se(σ̂ 2 ) =
n−3

où RSS est la somme des carrés des résidus et n est la taille de l’échantillon. Nous avons n − 3 degrés de
liberté car nous avons estimé trois paramètres (a0 , a1 , a2 ) pour construire le modèle de régression linéaire.
Les intervalles de confiance pour les paramètres ai et σ 2 peuvent être construits en utilisant la distribution
de Student t avec n − 3 degrés de liberté. Pour un niveau de confiance (1 − α), les intervalles de confiance à
(1 − α) pour les paramètres ai et σ 2 sont donnés par :

aˆi ± tα/2,n−3 × se(aˆi )

σ̂ 2 ± tα/2,n−3 × se(σ̂ 2 )

où tα/2,n−3 est le quantile de la distribution de Student avec n − 3 degrés de liberté tel que P(Tn−3 >
tα/2,n−3 ) = α/2.

# Données
# Remplacez les valeurs manquantes avec les données appropriées
m <- c(0, 0)
Gamma <- matrix(c(1, 0.5, 0.5, 1), nrow = 2, ncol = 2)
sigma2 <- 1
a <- c(1, 2, 3)
n <- 100
alpha <- 0.05

# Simulation des données

library(MASS)
[Link](123)
X <- mvrnorm(n, m, Gamma)
epsilon <- rnorm(n, mean = 0, sd = sqrt(sigma2))
Y <- a[1] + a[2] * X[,1] + a[3] * X[,2] + epsilon

# Modèle de régression
fit <- lm(Y ~ X[,1] + X[,2])

26
2.10 Question 10 : Calculez la région de confiance simultanée des paramètres de la regression
a = (a0 , a1 , a2 )T

# Intervalles de confiance pour les paramètres a

confint(fit, level = 1 - alpha)

## 2.5 % 97.5 %
## (Intercept) 0.944254 1.325877
## X[, 1] 1.674125 2.124479
## X[, 2] 2.711590 3.182260

# Intervalle de confiance pour sigma2

sigma2_ci <- c( sum(resid(fit)ˆ2) / qchisq(1 - alpha/2, n - 3),
sum(resid(fit)ˆ2) / qchisq(alpha/2, n - 3))
sigma2_ci

## [1] 0.6959003 1.2252932

L’intervalle de confiance à 95 pour les paramètres a0 , a1 et a2 sont respectivement [0.944, 1.326], [1.674, 2.124]
et [2.712, 3.182]. Cela signifie qu’on estime avec 95 de confiance que la vraie valeur de chaque paramètre se
trouve dans l’intervalle correspondant.
L’intervalle de confiance à 95 pour la variance σ 2 est [0.696, 1.225]. Cela signifie qu’on estime avec 95 de
confiance que la vraie valeur de la variance se trouve dans cet intervalle. Plus précisément, on peut dire que
la variance des erreurs de la régression est probablement comprise entre ces deux valeurs.

2.10 Question 10 : Calculez la région de confiance simultanée des paramètres

de la regression a = (a0 , a1 , a2 )T

La région de confiance simultanée pour les paramètres de la régression est donnée par l’intersection des inter-
valles de confiance individuels pour chaque paramètre. Pour cela, on peut utiliser la méthode de Bonferroni
qui consiste à diviser le niveau de confiance α par le nombre de paramètres à estimer.
La formule pour l’intervalle de confiance Bonferroni pour le i-ème paramètre est :

h q q i
aˆi − t 2p
α
,n−p s2ii , aˆi − t 2p
α
,n−p s2ii

où aˆi est l’estimateur du i-ème paramètre, s2ii est l’i-ème élément de la diagonale de la matrice de covariance
des estimateurs, p est le nombre de paramètres et n − p est le nombre de degrés de liberté de l’erreur.
Le code pour calculer la région de confiance simultanée pour les paramètres de la régression à l’aide de la
méthode de Bonferroni est :

# Région de confiance simultanée pour les paramètres de la régression

alpha <- 0.05
n <- 100
p <- 3
se_a <- sqrt(diag(vcov(fit)))
F_value <- qf(1 - alpha, p, n - p)
Ct <- sqrt((p*(n - p)*F_value)/(n*(n - p + p*F_value)))
delta_a <- Ct*se_a
conf_region_a <- cbind(coef(fit) - delta_a, coef(fit) + delta_a)
conf_region_a

27
2.11 Question 11 :

## [,1] [,2]
## (Intercept) 1.108786 1.161345
## X[, 1] 1.868290 1.930314
## X[, 2] 2.914514 2.979336

2.11 Question 11 :

En supposant que X1 et X2 sont des variables normales suivant respectivement N (1 , 0.5) et N (2 , 1.5),
générer Y avec a = (12, 6.5, 4).
Pour générer Y avec les valeurs données pour a et les distributions normales données pour X1 et X2 , nous
pouvons utiliser la formule de régression linéaire multiple :

Y = Xβ + ϵ

où X est une matrice de taille n×p contenant les variables explicatives (dans ce cas X1 et X2 ) et une colonne
de 1 pour l’ordonnée à l’origine, β est un vecteur de taille p + 1 contenant les coefficients de régression (dans
ce cas a0 , a1 et a2 ) et ϵ est un vecteur de bruit blanc de taille n.
Dans ce cas, n = 1000 pour générer 1000 échantillons de Y , p = 2 car nous avons deux variables explicatives
(X1 et X2 ), et X est donc une matrice de taille 1000 × 3 avec une première colonne de 1 pour l’ordonnée à
l’origine.
 
γ11 γ12
La matrice de covariance Γ est donnée par Γ = , où γ11 est la variance de X1 , γ22 est la variance
γ21 γ22
de X2 et γ12 = γ21 est la covariance entre X1 et X2 . Dans ce cas, nous avons γ11 = 0.5, γ22 = 1.5 et
γ12 = γ21 = 0 car X1 et X2 sont indépendants.
La matrice de variance covariance de ϵ est donnée par σ 2 I, où I est la matrice identité de taille n. Dans ce
cas, nous avons σ 2 = 1.

# Définir les paramètres

a <- c(12, 6.5, 4) # vecteur des coefficients de régression
sigma2 <- 2 # variance de bruit blanc

# Générer les variables X1 et X2

X1 <- rnorm(1, 1, 0.5)
X2 <- rnorm(1, 2, 1.5)

# Calculer la matrice de covariance

Sigma <- diag(rep(sigma2, 3))

# Générer le vecteur de bruit blanc centré

epsilon <- mvrnorm(n = 1, mu = rep(0, 3), Sigma = Sigma)

# Calculer le vecteur Y
X <- c(1, X1, X2)
Y <- [Link](X %*% a) + [Link](epsilon)

# Afficher le résultat
cat("Y =", Y)

## Y = 19.04109 18.17085 18.33203

28
2.12 Question 12 : Effectuer la représentation géométrique de la relation.

2.12 Question 12 : Effectuer la représentation géométrique de la relation.

La relation de régression linéaire entre Y et les variables explicatives X1 et X2 peut être représentée
géométriquement en utilisant un plan cartésien tridimensionnel, où Y est la variable réponse et X1 et X2
sont les variables explicatives.
La forme générale de l’équation de régression linéaire multiple est donnée par :

Y = a0 + a1 X1 + a2 X2 + ϵ

En utilisant les coefficients de régression a = (12, 6.5, 4) de la question 5, nous pouvons écrire l’équation
spécifique de notre modèle de régression comme suit :

Y = 12 + 6.5X1 + 4X2 + ϵ

Nous pouvons maintenant représenter cette équation sur un plan cartésien tridimensionnel en utilisant les
coordonnées (X1 , X2 , Y ). Pour cela, nous pouvons créer une surface de régression qui représente la relation
entre Y et les variables explicatives X1 et X2 .

library(plotly)

# Définir les paramètres

a <- c(12, 6.5, 4) # vecteur des coefficients de régression
sigma2 <- 2 # variance de bruit blanc

# Générer les variables X1 et X2

X1 <- rnorm(100, 1, 0.5)
X2 <- rnorm(100, 2, 1.5)

# Calculer la matrice de covariance

Sigma <- diag(rep(sigma2, 3))

# Générer le vecteur de bruit blanc centré

epsilon <- mvrnorm(n = 100, mu = rep(0, 3), Sigma = Sigma)

# Calculer la variable Y
X <- cbind(rep(1, 100), X1, X2)
Y_lin <- X %*% a # Partie linéaire de Y
Y <- Y_lin + matrix(epsilon[,1],
nrow = nrow(X), ncol = 1) # Ajouter le bruit blanc centré

# Transformer Y en une matrice de dimension 100 x 1

Z <- matrix(rep(Y[,1], 3), ncol = 3)

# Créer le nuage de points interactif

plot_ly(x = X1, y = X2, z = Y[,1], type = "scatter3d", mode = "markers") %>%
add_surface(x = X1, y = X2, z = matrix(Y_lin[,1],
ncol = 1), showscale = FALSE, opacity = 2)

29
2.13 Question 13 : Effectuer les estimations des paramètres, des intervalles et de la région de confiance
sur les données simulées.

2.13 Question 13 : Effectuer les estimations des paramètres, des intervalles et

de la région de confiance sur les données simulées.

Pour estimer les paramètres du modèle de régression linéaire multiple, on peut utiliser la fonction lm() de
R. Cette fonction permet de spécifier le modèle à estimer en utilisant une formule qui décrit la relation entre
les variables dépendantes et indépendantes. Dans notre cas, la formule sera Y ∼ X1 + X2 .

# Estimation des paramètres

fit <- lm(Y ~ X1 + X2)
summary(fit)

##
## Call:
## lm(formula = Y ~ X1 + X2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0208 -0.6831 0.1023 0.8970 3.6131
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.13663 0.36693 33.08 <2e-16 ***
## X1 6.87981 0.26905 25.57 <2e-16 ***
## X2 3.87145 0.09232 41.93 <2e-16 ***
## ---
## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
##
## Residual standard error: 1.383 on 97 degrees of freedom
## Multiple R-squared: 0.9601, Adjusted R-squared: 0.9593
## F-statistic: 1166 on 2 and 97 DF, p-value: < 2.2e-16

30
2.13 Question 13 : Effectuer les estimations des paramètres, des intervalles et de la région de confiance
sur les données simulées.

La fonction summary() appliquée à l’objet fit permet d’afficher les résultats de l’estimation, y compris les
estimations des coefficients, les intervalles de confiance et les statistiques de test.
Pour obtenir les intervalles de confiance et la région de confiance, on peut utiliser la fonction confint():

# Intervalles de confiance et région de confiance

confint(fit)

## 2.5 % 97.5 %
## (Intercept) 11.408382 12.864880
## X1 6.345824 7.413787
## X2 3.688211 4.054684

Cette fonction renvoie les intervalles de confiance pour chaque coefficient du modèle, ainsi que la région de
confiance pour l’ensemble des coefficients.

31
3 Exercice 3:
On considère une variable aléatoire suivant une loi normale

N (µ, σ 2 )

3.1 Question 1 : Densité et fonction de répartition de N (µ = 5, σ 2 = 25)

La densité de probabilité d’une variable aléatoire suivant une loi normale N (µ, σ 2 ) est donnée par :

1 (x−µ)2
f (x) = √ e− 2σ 2
2πσ 2

où µ est la moyenne de la distribution et σ est l’écart-type.

La fonction de répartition de la variable aléatoire suivant une loi normale N (µ, σ 2 ) est donnée par :

1
  
x−µ
F (x) = 1 + erf √
2 σ 2

où erf (x) est la fonction d’erreur.

mu <- 5
sigma <- sqrt(25)

# Création d'une séquence de valeurs pour la densité

x <- seq(mu - 3 * sigma, mu + 3 * sigma, length = 100)

# Densité de la loi normale

pdf <- dnorm(x, mean = mu, sd = sigma)

# Fonction de répartition de la loi normale

cdf <- pnorm(x, mean = mu, sd = sigma)

# Tracer la densité et la fonction de répartition

par(mfrow = c(1,2)) # diviser la fenêtre de tracé en deux
plot(x, pdf, type="l", main="Densité", xlab="x", ylab="f(x)")
plot(x, cdf, type="l",
main="Fonction rép",
xlab="x", ylab="F(x)")

32
3.2 Question 2 : Densité et fonction de répartition de la forme centrée et réduite de N (µ = 5, σ 2 = 25)

Densité Fonction rép

0.08

1.0
0.8
0.06

0.6
0.04

F(x)
f(x)

0.4
0.02

0.2
0.00

−10 0 5 10 20 0.0 −10 0 5 10 20

x x

3.2 Question 2 : Densité et fonction de répartition de la forme centrée et réduite

de N (µ = 5, σ 2 = 25)

La forme centrée et réduite de la loi normale est une distribution normale standardisée avec une moyenne
de 0 et un écart-type de 1. Sa densité de probabilité et sa fonction de répartition sont données par :
Densité de probabilité :

1 x2
f (x) = √ e− 2
2π

Fonction de répartition :

1
  
x
F (x) = 1 + erf √
2 2
mu <- 5
sigma <- sqrt(25)

# Création d'une séquence de valeurs pour la densité

x_std <- seq(-3, 3, length = 100)

# Densité de la loi normale centrée et réduite

pdf_std <- dnorm(x_std, mean = 0, sd = 1)

33
3.3 Question 3 : Effectuez une simulations d’un échantillon de taille n = 1000 issue de cette loi normale
N (µ = 2, σ 2 = 1.5)

# Fonction de répartition de la loi normale centrée et réduite

cdf_std <- pnorm(x_std, mean = 0, sd = 1)

# Tracer la densité et la fonction de répartition

par(mfrow = c(1,2)) # diviser la fenêtre de tracé en deux
plot(x_std, pdf_std, type="l",
main="Densité",
xlab="x", ylab="f(x)")
plot(x_std, cdf_std, type="l",
main="Fonction rép",
xlab="x", ylab="F(x)")

Densité Fonction rép

0.4

1.0
0.8
0.3

0.6
F(x)
f(x)

0.2

0.4
0.1

0.2
0.0
0.0

−3 −1 0 1 2 3 −3 −1 0 1 2 3

x x

3.3 Question 3 : Effectuez une simulations d’un échantillon de taille n = 1000

issue de cette loi normale N (µ = 2, σ 2 = 1.5)

# Simuler un échantillon de taille 1000 issu de la loi normale N(mu=2, sigmaˆ2=1.5)

n <- 1000
mu <- 2
sigma <- sqrt(1.5)
[Link](2)
echantillon <- rnorm(n, mean=mu, sd=sigma)

34
3.4 Question 4 :

# Afficher les premières valeurs de l'échantillon

head(echantillon)

## [1] 0.9015085 2.2263931 3.9447054 0.6155782 1.9017121 2.1621811

3.4 Question 4 :
Créer une fonction appelée LV de la variable θ = (µ, σ 2 ) qui calcule la log vraisemblance en θ de l’échantillon
précédant et calculer LV (θ) pour θ = (2, 1) et pour θ = (2.5, 1.5)
La log-vraisemblance est une fonction qui mesure la plausibilité d’un échantillon de données, vu un modèle
statistique paramétrique. Elle est souvent utilisée pour estimer les paramètres du modèle à partir des données
observées. Plus précisément, la log-vraisemblance d’un échantillon est la somme des logarithmes des densités
de probabilité (ou des fonctions de masse de probabilité) des observations, évaluées pour les valeurs de
paramètres considérées. En d’autres termes, elle mesure la probabilité que les observations soient générées
à partir du modèle avec les paramètres donnés.
Pour un échantillon de n observations, la log-vraisemblance en un vecteur de paramètres θ est donnée par :

n
X
LV (θ) = log f (Xi |θ)
i=1

où f (Xi |θ) est la densité de probabilité (ou la fonction de masse de probabilité) de l’observation Xi pour les
paramètres θ.
Maintenant, pour créer une fonction LV pour la variable θ = (µ, σ 2 ), on peut utiliser la densité de probabilité
de la loi normale pour chaque observation.

3.4.1 Logvraisemblance

LV <- function(x, mu, sigma2) {

n <- length(x)
log_likelihood <- -n/2*log(2*pi) - n/2*log(sigma2) - sum((x-mu)ˆ2)/(2*sigma2)
return(log_likelihood)
}

3.4.2 Calcul des valeurs

# Calcul de LV pour theta = (2,1)

theta1 <- c(2, 1)
LV1 <- LV(echantillon, mu=theta1[1], sigma2=theta1[2])
cat("LV(theta1) = ", LV1, "\n")

## LV(theta1) = -1693.636

# Calcul de LV pour theta = (2.5, 1.5)

theta2 <- c(2.5, 1.5)
LV2 <- LV(echantillon, mu=theta2[1], sigma2=theta2[2])
cat("LV(theta2) = ", LV2, "\n")

## LV(theta2) = -1696.159

35
3.5 Question 5 : Représentez la fonction de vraisemblance

3.4.3 Lequel choisir des deux

# Comparaison des log-vraisemblances

if (LV1 > LV2) {
cat("Le paramètre theta1 = (2, 1) est plus vraisemblable que theta2 = (2.5, 1.5)\n")
} else {
cat("Le paramètre theta2 = (2.5, 1.5) est plus vraisemblable que theta1 = (2, 1)\n")
}

## Le paramètre theta1 = (2, 1) est plus vraisemblable que theta2 = (2.5, 1.5)

3.5 Question 5 : Représentez la fonction de vraisemblance

Pour représenter la log vraisemblance, nous allons exécuter le code suivant :

# Créer une grille de valeurs pour mu et sigma

mu_grid <- seq(1, 3, [Link]=100)
sigma2_grid <- seq(1, 2, [Link]=100)

# Fonction pour calculer et tracer la fonction de vraisemblance pour une loi normale avec un échantillon
likelihood_normal <- function(mu, sigma2) {
n <- 100
x <- rnorm(n, mu, sqrt(sigma2))
likelihoods <- rep(NA, length(mu)*length(sigma2))
k <- 1
for (i in 1:length(mu)) {
for (j in 1:length(sigma2)) {
likelihoods[k] <- LV(x, mu[i], sigma2[j])
k <- k + 1
}
}
likelihoods <- exp(likelihoods - max(likelihoods))
likelihoods <- matrix(likelihoods, nrow=length(mu))
colnames(likelihoods) <- round(sigma2, 2)
rownames(likelihoods) <- round(mu, 2)
plot_ly(x=mu, y=sqrt(sigma2), z=likelihoods, type="surface") %>%
add_surface(contours=list(z=list(show=TRUE, usecolormap=TRUE,
highlightcolor="#ff0000", project=list(z=TRUE))))%>%
layout(scene = list(xaxis = list(title = "Mu"),
yaxis = list(title = "Sigma"),
zaxis = list(title = "Log-Likelihood")))
}

likelihood_normal(mu_grid,sigma2_grid)

36
3.6 Question 6. Cherchez l’estimateur du maximum de vraisemblance, θV M pour θ

3.6 Question 6. Cherchez l’estimateur du maximum de vraisemblance, θV M

pour θ

L’estimateur du maximum de vraisemblance (EMV) est la valeur de θ qui maximise la fonction de vraisem-
blance L(θ|x1 , . . . , xn ), où x1 , . . . , xn sont les réalisations de l’échantillon.
Pour une loi normale N (µ, σ 2 ), les estimateurs EMV de µ et σ 2 sont donnés respectivement par :

n
1X
µ̂EM V = xi
n i=1

et

n
1X
V = (xi − µ̂EM V )2
2
σ̂EM
n i=1

où x1 , . . . , xn sont les réalisations de l’échantillon.

# Fonction de log-vraisemblance pour maximisation

# Définir la fonction de log-vraisemblance pour l'échantillon de taille 100
LV_n <- function(theta, x){
mu <- theta[1]
sigma2 <- theta[2]
n <- length(x)
log_likelihood <- -n/2*log(2*pi) - n/2*log(sigma2) -
sum((x-mu)ˆ2)/(2*sigma2)
return(-log_likelihood)
}

# Trouver l'EMV de theta en utilisant la fonction optim

37
3.7 Question 7 : Effectuez une estimation par intervalle de confiance à 95% de la moyenne µ et de la
variance σ 2 par calcul et à l’aide de fonction prédéfinies

theta_EMV <- optim(c(0, 1), LV_n, x=echantillon, method="BFGS")$par

# Afficher l'estimateur de l'EMV de theta

theta_EMV

## [1] 2.075932 1.543638

3.7 Question 7 : Effectuez une estimation par intervalle de confiance à 95% de la

moyenne µ et de la variance σ 2 par calcul et à l’aide de fonction prédéfinies

Pour calculer l’intervalle de confiance à 95% de la moyenne et de la variance de l’échantillon, on peut utiliser
les formules suivantes :
Pour l’intervalle de confiance de la moyenne :

• borneinf = moyenne − t(n−1, α2 ) √σn

• bornesup = moyenne + t(n−1, α2 ) √σn

où t(n−1, α2 ) est le quantile de la distribution de Student à n − 1 degrés de liberté pour une probabilité α2 ,
ecart_type est l’écart-type de l’échantillon, n est la taille de l’échantillon et moyenne est la moyenne de
l’échantillon.
Pour l’intervalle de confiance de la variance :

2
• borneinf = (n − 1) χ(1−α,n−1)
σ

2
• borneinf = (n − 1) χ(α,n−1)
σ

où
χ(p,k)
est le quantile de la distribution du χ2 à k degrés de liberté pour une probabilité p, variance est la variance
de l’échantillon, n est la taille de l’échantillon.

# Echantillon
# Estimation par intervalle de confiance de la moyenne
mu_ci <- [Link](echantillon, [Link] = 0.95)$[Link]

# Fit a normal distribution to the data

fit <- fitdistr(echantillon, "normal")

# Calculate the confidence interval for the variance

n <- length(echantillon)
df <- n - 1
alpha <- 0.05
var_est <- var(echantillon)
chi2_left <- qchisq(1-alpha/2, df)
chi2_right <- qchisq(alpha/2, df)
var_ci <- c((df*var_est)/chi2_left, (df*var_est)/chi2_right)

mu_ci

38
3.8 Question 8 : Calculez la région de confiance simultanée des deux paramètres inconnus de θ.

## [1] 1.998796 2.153070

## attr(,"[Link]")
## [1] 0.95

var_ci

## [1] 1.418142 1.690165

3.8 Question 8 : Calculez la région de confiance simultanée des deux paramètres

inconnus de θ.

La région de confiance simultanée pour les deux paramètres inconnus µ et σ 2 de θ est l’ensemble des paires
(µ, σ 2 ) pour lesquelles on rejette l’hypothèse nulle que θ = (µ0 , σ02 ), avec une probabilité globale de 1 − α où
α est le niveau de significativité.
Pour calculer la région de confiance simultanée pour les paramètres µ et σ 2 , on peut utiliser la méthode
de l’hypersphère de Bonferroni, qui consiste à construire des régions de confiance pour chaque paramètre
individuellement, en utilisant le test de Wald, et à prendre l’intersection de ces régions.
Le test de Wald pour la moyenne µ est basé sur la statistique de test :

µ̂ − µ0
q ∼ N (0, 1)
V ar(µ̂)
d

où µ̂ est l’estimateur de la moyenne, µ0 est la valeur nulle à tester, et Vd

ar(µ̂) est l’estimateur de la variance
de µ̂.
De même, le test de Wald pour la variance σ 2 est basé sur la statistique de test :

(n − 1)σ̂ 2
∼ χ2n−1
σ02

où σ̂ 2 est l’estimateur de la variance, σ02 est la valeur nulle à tester, et χ2n−1 est la distribution chi-carré avec
n − 1 degrés de liberté.
Les régions de confiance individuelles pour µ et σ 2 sont données par :

(µ̂ − µ0 )2 t2n−1, α
≤ 2

Vd ar(µ̂) n−1

et

(n − 1)σ̂ 2 (n − 1)σ̂ 2
2 ≤ σ02 ≤ 2
χn−1, α χn−1,1− α
2 2

où tn−1, α2 et χ2n−1, α et χ2n−1,1− α sont les quantiles des distributions de Student t et chi-carré, respectivement,
2 2
avec n − 1 degrés de liberté et un niveau de significativité α2 .
L’intersection de ces deux régions de confiance donne la région de confiance simultanée pour µ et σ 2 .

39
3.8 Question 8 : Calculez la région de confiance simultanée des deux paramètres inconnus de θ.

mu_EMV <- mean(echantillon)

sigma2_EMV <- var(echantillon)
V <- matrix(c(sigma2_EMV/n, 0, 0, 2*sigma2_EMV/nˆ2), nrow = 2)

## Création de la matrice de variance-covariance de theta_EMV

V <- matrix(c(sigma2_EMV/n, 0, 0, 2*sigma2_EMV/nˆ2), nrow = 2)

## Détermination des valeurs propres et vecteurs propres de V

eigen_decomp <- eigen(V)

## Extraction des valeurs propres et vecteurs propres

eigen_values <- eigen_decomp$values
eigen_vectors <- eigen_decomp$vectors

## Détermination des axes de la région de confiance

F_quantile <- qf(0.95, 2, n - 2)
theta_EMV_centered <- theta_EMV - c(mu_EMV, 1)
theta_EMV_rotated <- eigen_vectors %*% theta_EMV_centered
theta_EMV_rotated_x <- sqrt(eigen_values[1] * F_quantile) *
cos(seq(0, 2*pi, [Link] = 100))
theta_EMV_rotated_y <- sqrt(eigen_values[2] * F_quantile) *
sin(seq(0, 2*pi, [Link] = 100))
theta_EMV_ellipse <- cbind(theta_EMV_rotated_x, theta_EMV_rotated_y)
theta_EMV_ellipse <- t(eigen_vectors %*% t(theta_EMV_ellipse) + c(mu_EMV, 1))

## Tracé de la région de confiance

plot(echantillon, main = "Région de confiance simultanée pour theta", xlab = "Index de l'échantillon", y
points(theta_EMV[1], theta_EMV[2], col = "red", pch = 19)
lines(theta_EMV_ellipse, col = "red", lwd = 2)

40
3.9 Question 9

Région de confiance simultanée pour theta

5
Valeur de l'échantillon

4
3
2
1
0
−1

0 200 400 600 800 1000

Index de l'échantillon

3.9 Question 9

On suppose désormais que l’on dispose d’une information a priori sur le paramètre θ, réduit à sa première
composante θ = µ sachant que sa deuxième composante initiale σ = 1, résumée dans la loi inverse-gamma
de densité

βα 1 −β
f (θ) = e θ 1R+ (θ)
Γ(α) θ α+1

3.10 Question 10 :

Déterminer la loi a posteriori et générer une série de taille n = 1000 issue de cette loi pour des paramètres
que vous proposerez.
La loi a posteriori pour µ sachant l’information a priori dans la loi inverse-gamma peut être calculée en util-
isant le théorème de Bayes. Supposons que nous avons un échantillon X1 , X2 , . . . , Xn tiré de la distribution
normale N (µ, σ 2 ) avec σ 2 connu. La densité de probabilité conjointe de l’échantillon et du paramètre µ est :

f (µ|x) ∝ f (x|µ)f (µ)

où f (x|µ) est la vraisemblance de l’échantillon et f (µ) est la densité de probabilité a priori pour µ dans la
loi inverse-gamma.
La vraisemblance de l’échantillon est donnée par :

41
3.10 Question 10 :

n n
!
Y 1 1 X
f (x|µ) = f (xi |µ) = √ exp − 2 (xi − µ)2
i=1 ( 2πσ 2 )n 2σ i=1

La densité de probabilité a posteriori pour µ sachant l’information a priori dans la loi inverse-gamma est
alors donnée par :

n
!
1 X βα 1 − µβ
f (µ|x) ∝ exp − 2 (xi − µ)2 × e
2σ i=1 Γ(α) µα+1

En normalisant cette densité de probabilité, on obtient :

( σβ2 + nx̄)α 1 − β+σ

2 nx̄2

f (µ|x) = e σ2 µ
Γ(α) µα+n+1

où x̄ est la moyenne de l’échantillon.

Pour générer une série de taille n = 1000 issue de cette loi, nous pouvons utiliser la méthode de rejet. Tout
d’abord, nous choisissons des paramètres α et β pour la loi inverse-gamma a priori. Ensuite, nous pouvons
générer une série de taille n = 1000 en utilisant la densité de probabilité a posteriori pour µ calculée ci-dessus.
Voici le code pour générer une série de taille n = 1000 issue de la loi a posteriori pour µ sachant l’information
a priori dans la loi inverse-gamma avec α = 2 et β = 2 :

# Définition de la loi a posteriori

mu_posteriori <- function(mu, data, alpha, beta) {
n <- length(data)
s2 <- var(data)
alpha_post <- alpha + n/2
beta_post <- beta + (n-1)*s2/2 + (n*(mu-mean(data))ˆ2)/(2*(n+1))
dgamma(mu, shape = alpha_post, rate = beta_post)
}

# Paramètres de la loi a priori

alpha <- 2
beta <- 2

# Données simulées de la loi normale

[Link](123) # pour rendre les résultats reproductibles
mu_true <- 2.5
sigma_true <- 1.2
data <- rnorm(n = 1000, mean = mu_true, sd = sigma_true)

# Génération d'un échantillon de taille n=1000 de la loi a posteriori

mu_post <- seq(from = 1, to = 4, [Link] = 1000)
posteriori_density <- mu_posteriori(mu_post, data, alpha, beta)
posteriori_density <- posteriori_density /
integrate(mu_posteriori, lower = -Inf, upper = Inf,
data = data, alpha = alpha, beta = beta)$value
mu_sample <- sample(mu_post, size = 1000,
replace = TRUE, prob = posteriori_density)

42
3.11 Question 11

3.11 Question 11

Déterminer les estimateur Bayésien d’erreur quadratique moyen minimale θM M SE et du maximum a pos-
teriori θM AP du paramètre θ et effectuer l’estimation avec la serie simulée.
Pour déterminer l’estimateur Bayésien d’erreur quadratique moyenne minimale θ̂M M SE , nous devons calculer
l’espérance conditionnelle de θ sachant les observations. Avec la loi a posteriori que nous avons déterminée
précédemment, l’espérance conditionnelle est donnée par :

α+n
θ̂M M SE = E[θ|x1 , . . . , xn ] = Pn
β + i=1 (xi − x̄)2

où x̄ est la moyenne des observations x1 , . . . , xn .

Pour déterminer l’estimateur du maximum a posteriori (MAP) θ̂M AP , nous devons trouver la valeur de θ
qui maximise la loi a posteriori. La loi a posteriori que nous avons déterminée est une loi inverse-gamma,
donc nous avons :

n
βα 1 −β Y 1 (xi −θ)2
θ̂M AP = argmax θf (θ|x1 , . . . , xn ) = argmax θ e θ √ e− 2θ
Γ(α) θ α+1
i=1
2πθ

En utilisant le logarithme pour simplifier les calculs, nous pouvons écrire :

n
βα β n 1X (xi − θ)2
log f (θ|x1 , . . . , xn ) = log − (α + 1) log θ − − log(2π) − log θ −
Γ(α) θ 2 2 i=1 2θ

Nous pouvons maximiser cette expression en trouvant la valeur de θ qui annule sa dérivée par rapport à θ.
Après quelques calculs, nous obtenons :

β + n/2
θ̂M AP = Pn
α+ i=1 (xi − x̄) /2
2

où x̄ est la moyenne des observations x1 , . . . , xn .

Maintenant, nous allons effectuer l’estimation avec la série simulée. Nous allons générer une série de taille
n = 1000 pour les paramètres α = 2 et β = 1, comme demandé dans la question précédente.

# Estimation de theta_MMSE
alpha_n <- alpha + n/2
beta_n <- beta + 1/2 * sum((mu_sample - mu_EMV)ˆ2)
theta_MMSE <- beta_n / (alpha_n + 1)

# Estimation de theta_MAP
theta_MAP <- beta_n / (alpha_n + 1)

# Affichage des résultats

cat("Estimation Bayésienne:\n")

## Estimation Bayésienne:

cat("theta_MMSE =", theta_MMSE, "\n")

## theta_MMSE = 1.147631

43
3.12 Question 12 : Créer une fonction qui permet de calculer θ̂M M SE et θ̂M AP .

cat("theta_MAP =", theta_MAP, "\n")

## theta_MAP = 1.147631

Dans ce code, on calcule d’abord les nouvelles valeurs des paramètres de la loi inverse-gamma a posteriori,
à partir de l’information a priori et des données simulées. Ensuite, on estime l’erreur quadratique moyenne
minimale et le maximum a posteriori à partir des nouvelles valeurs des paramètres de la loi a posteriori.
Enfin, on affiche les résultats.
Notons que pour ce code, on suppose que les données simulées sont stockées dans le vecteur x et que
l’estimateur de la moyenne est stocké dans la variable µ̂. Les valeurs des paramètres de la loi inverse-gamma
a priori sont stockées dans les variables alpha et beta.

3.12 Question 12 : Créer une fonction qui permet de calculer θ̂M M SE et θ̂M AP .

Ayant les paramètres de la loi inverse-gamma a priori α et β, la fonction qui permet de calculer les estimateurs
Bayésiens d’erreur quadratique moyen minimale et du maximum a posteriori pour un échantillon de données
x est la suivante :

bayesian_estimation <- function(x, alpha, beta) {

n <- length(x)
mu_hat <- mean(x)

alpha_n <- alpha + n/2

beta_n <- beta + 1/2 * sum((x - mu_hat)ˆ2)

theta_MMSE <- beta_n / (alpha_n + 1)

theta_MAP <- beta_n / (alpha_n + 1)

list(theta_MMSE = theta_MMSE, theta_MAP = theta_MAP)

}

Cette fonction utilise les mêmes expressions pour calculer les estimateurs que celles que nous avons vues
précédemment. Il suffit de fournir les paramètres alpha et beta de la loi inverse-gamma a priori, ainsi que les
données simulées et l’estimateur de la moyenne pour obtenir les estimateurs Bayésiens d’erreur quadratique
moyen minimale et du maximum a posteriori du paramètre theta.
On peut faire un test avec α = 2 et β = 2

bayesian_estimation(echantillon,2,2)

## $theta_MMSE
## [1] 1.538399
##
## $theta_MAP
## [1] 1.538399

44