Projet SY02
Eduardo KRIK ANTUNES DE LIMA
Larissa CABIDO RUFFO
Matheus AVELLAR RODRIGUES
2025-01-10
Projet avec R SY02
Exercice 1 - Les deux composants électroniques
Introduction
Nous étudions les durées de vie de deux composants électroniques. Les durées de vie sont modélisées par deux variables aléatoires indépendantes \(X \sim \text{Exp}(\lambda)\) et \(Y \sim \text{Exp}(\mu)\). Le système doit être remplacé dès qu’un des deux composants tombe en panne. Nous observons les variables \(Z = \min(X, Y)\) et \(W\) définies comme suit : \[ W = \begin{cases} 1 & \text{si } Z = X, \\ 0 & \text{si } Z = Y. \end{cases} \]
Question 1 : Interprétation des variables \(Z\) et \(W\)
Variable \(Z\) : Elle représente la durée de vie minimale entre les deux composants \(X\) et \(Y\), car le système tombe en panne dès qu’un composant échoue.
Variable \(W\) : Elle indique quel composant est tombé en panne en premier :
- \(W = 1\) si le composant \(X\) tombe en panne en premier.
- \(W = 0\) si le composant \(Y\) tombe en panne en premier.
Question 2 : Lois de \(Z\) et \(W\)
Partie 1 : Lois des variables aléatoires
Loi de \(Z\)
\(Z = \min(X, Y)\). Puisque \(X \sim \text{Exp}(\lambda)\) et \(Y \sim \text{Exp}(\mu)\), la fonction de répartition de \(Z\) est donnée par : \[ F_Z(z) = P(Z \leq z) = P(\min(X, Y) \leq z) = 1 - P(X > z \cap Y > z). \] Les deux événements \(X > z\) et \(Y > z\) sont indépendants : \[ P(X > z \cap Y > z) = P(X > z) P(Y > z). \] En utilisant les propriétés des lois exponentielles : \[ P(X > z) = e^{-\lambda z}, \quad P(Y > z) = e^{-\mu z}. \] Donc : \[ F_Z(z) = 1 - e^{-\lambda z} e^{-\mu z} = 1 - e^{-(\lambda + \mu)z}. \] Cela montre que \(Z \sim \text{Exp}(\lambda + \mu)\).
Loi de \(W\)
La probabilité que \(W = 1\) est donnée par : \[ P(W = 1) = P(Z = X) = P(X < Y) = \int_0^\infty P(X < Y | X = x) f_X(x) dx. \] En utilisant les densités des lois exponentielles, il est possible de montrer que : \[ P(W = 1) = \frac{\lambda}{\lambda + \mu}, \quad P(W = 0) = \frac{\mu}{\lambda + \mu}. \] Ainsi, \(W\) suit une loi Bernoulli de paramètre \(p = \frac{\lambda}{\lambda + \mu}\).
Partie 2 : Simulation en R
set.seed(120)
# Parametres pour les lois choisies arbitrairement
lambda <- 3
mu <- 5
n <- 1000
# Simulation pour les lois de la question
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y) # Duree de vie minimale entre les deux composants
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu)) # Composant est tombe en panne en premier# Histogramme de Z
hist(Z, breaks = 50, freq = FALSE, col = "lightblue", main = "Histogramme de Z",
xlab = "Z (duree de vie du systeme)")
curve(dexp(x, rate = lambda + mu), add = TRUE, col = "red", lwd = 2)
# Legende pour le graphique
legend("topright", legend = c("Histogramme de Z", "Courbe theorique (loi exponentielle)"),
col = c("lightblue", "red"), lwd = 2, bty = "n")
text(x = 2, y = 0.4, labels = "La courbe rouge = loi exponentielle de parametre (lambda + mu)", col = "red")
# Probabilites theoriques et empiriques
cat("Probabilites pour W = 1 : Theorique =", (lambda / (lambda + mu)), ", Empirique =", (mean(W == 1)), "\n")## Probabilites pour W = 1 : Theorique = 0.375 , Empirique = 0.378cat("Probabilites pour W = 0 : Theorique =", (mu / (lambda + mu)), ", Empirique =", (mean(W == 0)), "\n")## Probabilites pour W = 0 : Theorique = 0.625 , Empirique = 0.622Les résultants ont été validés, et le graphique démontre également l’exponentielle.
Question 3 : Estimation par maximum de vraisemblance (EMV) de λ
Partie 1 : Theorie
Définition de la vraisemblance
La vraisemblance totale des variables \((Z_1, W_1), \dots, (Z_n, W_n)\) est donnée par le produit des vraisemblances pour \(Z_i\) et \(W_i\), en supposant leur indépendance :
\[ L(\lambda; Z_1, \dots, Z_n, W_1, \dots, W_n) = L(\lambda; Z_1, \dots, Z_n) \cdot L(\lambda; W_1, \dots, W_n), \]
où :
- \(Z_i \sim \text{Exponentielle}(\lambda)\), avec \(f_{Z_i}(z) = \lambda e^{-\lambda z}\) pour \(z \geq 0\),
- \(W_i \sim \text{Bernoulli}\left(\frac{\lambda}{\lambda + \mu}\right)\).
Vraisemblance pour \(Z_i\)
La vraisemblance totale des variables \((Z_1, W_1), \dots, (Z_n, W_n)\) est donnée par le produit des vraisemblances pour \(Z_i\) et \(W_i\), en supposant leur indépendance :
\[ L(\lambda; Z_1, \dots, Z_n, W_1, \dots, W_n) = L(\lambda; Z_1, \dots, Z_n) \cdot L(\lambda; W_1, \dots, W_n), \]
où :
- \(Z_i \sim \text{Exponentielle}(\lambda)\), avec \(f_{Z_i}(z) = \lambda e^{-\lambda z}\) pour \(z \geq 0\)
- \(W_i \sim \text{Bernoulli}\left(\frac{\lambda}{\lambda + \mu}\right)\)
Vraisemblance pour \(Z_i\)
La fonction de vraisemblance pour \(Z_1, \dots, Z_n\) est donnée par :
\[ L(\lambda; Z_1, \dots, Z_n) = \prod_{i=1}^n \lambda e^{-\lambda Z_i} \]
Vraisemblance pour \(W_i\)
Pour \(W_1, \dots, W_n\), la vraisemblance est :
\[ L(\lambda; W_1, \dots, W_n) = \prod_{i=1}^n \left(\frac{\lambda}{\lambda + \mu}\right)^{W_i} \cdot \left(\frac{\mu}{\lambda + \mu}\right)^{1 - W_i} \]
En combinant les deux :
\[ L(\lambda) = \left(\lambda^n e^{-\lambda \sum_{i=1}^n Z_i}\right) \cdot \prod_{i=1}^n \left(\frac{\lambda}{\lambda + \mu}\right)^{W_i} \cdot \left(\frac{\mu}{\lambda + \mu}\right)^{1 - W_i} \]
Fonction de log-vraisemblance
Prenons le logarithme de la fonction de vraisemblance pour simplifier :
\[ \ell(\lambda) = n \ln(\lambda) - \lambda \sum_{i=1}^n Z_i + \sum_{i=1}^n \left(W_i \ln\left(\frac{\lambda}{\lambda + \mu}\right) + (1 - W_i) \ln\left(\frac{\mu}{\lambda + \mu}\right)\right) \]
Simplifions :
\[ \ell(\lambda) = n \ln(\lambda) - \lambda \sum_{i=1}^n Z_i + \sum_{i=1}^n W_i \ln(\lambda) - n \ln(\lambda + \mu) + \sum_{i=1}^n (1 - W_i) \ln(\mu) \]
Regroupons les termes :
\[ \ell(\lambda) = (n + \sum_{i=1}^n W_i) \ln(\lambda) - \lambda \sum_{i=1}^n Z_i - n \ln(\lambda + \mu) + \sum_{i=1}^n W_i \ln(\mu) \]
Estimation par Maximum de Vraisemblance
Pour maximiser la log-vraisemblance, nous dérivons \(\ell(\lambda)\) par rapport à \(\lambda\) et résolvons \(\frac{\partial \ell(\lambda)}{\partial \lambda} = 0\)
Dérivée de \(\ell(\lambda)\)
\[ \frac{\partial \ell(\lambda)}{\partial \lambda} = \frac{n + \sum_{i=1}^n W_i}{\lambda} - \sum_{i=1}^n Z_i - \frac{n}{\lambda + \mu} \]
En posant \(\frac{\partial \ell(\lambda)}{\partial \lambda} = 0\), nous obtenons :
\[ \hat{\lambda} = \frac{\sum_{i=1}^n W_i}{\sum_{i=1}^n Z_i} \]
Partie 2 : Simulation en R
set.seed(120)
# Parametres pour les lois
n <- 1000
mu <- 1
lambda <- 2
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
# Estimation par Maximum de Vraisemblance
lambda_hat <- sum(W) / sum(Z)## Valeur estimée de λ : 2.127941 , Valeur vraie de λ : 2L’estimateur du maximum de vraisemblance fournit une estimation proche de la vraie valeur de 𝜆, et cette estimation devient plus précise à mesure que 𝑛 augmente. Comme démontré ci-dessous :
set.seed(120)
# Parametres pour les lois
n <- 1000000
mu <- 1
lambda <- 2
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
# Estimation par Maximum de Vraisemblance
lambda_hat <- sum(W) / sum(Z)# Resultat
cat("Valeur estimée de λ : ", lambda_hat, ", Valeur vraie de λ : ", lambda, "\n")## Valeur estimée de λ : 1.99576 , Valeur vraie de λ : 2Question 4 : Biais, Variance et EQM de l’estimateur \(\hat{\lambda}_n\)
Partie 1: Theorie
Pour \(\lambda = 2\), nous allons calculer et représenter le biais, la variance et l’erreur quadratique moyenne (EQM) de l’estimateur du maximum de vraisemblance \(\hat{\lambda}_n\), en fonction de la taille de l’échantillon \(n\).
Définitions
Biais : La différence entre la valeur attendue de l’estimateur et la vraie valeur de \(\lambda\), c’est-à-dire \(\text{Biais}(\hat{\lambda}_n) = E[\hat{\lambda}_n] - \lambda\).
Variance : La variance de l’estimateur \(\hat{\lambda}_n\), c’est-à-dire \(\text{Var}(\hat{\lambda}_n) = E(\hat{\lambda}_n^2) - E(\hat{\lambda}_n)^2\)
Erreur quadratique moyenne (EQM) : L’Erreur quadratique moyenne est définie comme \(\text{EQM}(\hat{\lambda}_n) = \text{Biais}^2(\hat{\lambda}_n) + \text{Var}(\hat{\lambda}_n)\)
Partie 2 : Simulation en R
set.seed(120)
# Parametres pour les lois
lambda <- 2
mu <- 1
n_values <- seq(10, 500, by = 10) # Analyse avec differentes tailles d'échantillons
taille_n <- length(n_values) # Nombre d'analyses a faire
# Les vecteurs sont initialises avec le meme longueur
bias <- taille_n
variance <- taille_n
eqm <- taille_n
# Simulation pour differentes tailles d'echantillons
for (i in 1:taille_n) {
n <- n_values[i]
estimates <- 1000 # Nombre de simulations
for (j in 1:estimates) {
# Simuler le Z et le W
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
# Estimation de lambda
estimates[j] <- sum(W) / sum(Z)
}
# Calcul du biais, variance et EQM
bias[i] <- mean(estimates) - lambda
variance[i] <- var(estimates)
eqm[i] <- bias[i]^2 + variance[i]
}
# Representation graphique
plot(n_values, bias, type = "b", col = "blue", pch = 19,
main = "Biais, Variance et EQM de λ",
xlab = "Taille de l'echantillon (n)", ylab = "Valeurs",
ylim = range(c(bias, variance, eqm)))
# Ajouter la variance
lines(n_values, variance, type = "b", col = "green", pch = 19)
# Ajouter l'EQM
lines(n_values, eqm, type = "b", col = "purple", pch = 19)
# Ajouter une legende
legend("topright", legend = c("Biais", "Variance", "EQM"),
col = c("blue", "green", "purple"), pch = 19)
Le graphique précédent représente le biais, la variance et l’EQM de l’estimateur \(\hat{\lambda}_n\) pour différentes tailles d’échantillons \(n\).
- Biais (en bleu) : Le biais tend vers zéro à mesure que \(n\) augmente.
- Variance (en vert) : La variance diminue avec l’augmentation de \(n\), donc l’estimateur est plus précis.
- EQM (en violet) : L’EQM diminue également, indiquant une meilleure performance de l’estimateur avec un plus grand échantillon.
Donc, l’estimateur \(\hat{\lambda}_n\) est convergent, car le biais disparaît avec un grand \(n\) et la variance ainsi que l’EQM diminuent. Cela indique que l’estimateur devient de plus en plus précis à mesure que la taille de l’échantillon augmente, et l’estimateur est sans biais lorsque \(n\) est suffisamment grand.
Question 5 : Représenter, le biais, la variance et l’EQM de l’estimateur du maximum de vraisemblance \(\hat{\lambda}_n\)
Partie 1: Theorie
Nous souhaitons représenter le biais, la variance et l’erreur quadratique moyenne (EQM) de l’estimateur du maximum de vraisemblance \(\hat{\lambda}_n\), en fonction du paramètre \(\lambda\).
Biais : La différence entre la valeur attendue de l’estimateur et la vraie valeur de \(\lambda\), c’est-à-dire \(\text{Biais}(\hat{\lambda}_n) = E(\hat{\lambda}_n) - \lambda\)
Variance : La variance de l’estimateur \(\hat{\lambda}_n\), c’est-à-dire \(\text{Var}(\hat{\lambda}_n) = E(\hat{\lambda}_n^2) - E(\hat{\lambda}_n)^2\)
Erreur quadratique moyenne (EQM) : L’Erreur quadratique moyenne est définie comme \(\text{EQM}(\hat{\lambda}_n) = \text{Biais}^2(\hat{\lambda}_n) + \text{Var}(\hat{\lambda}_n)\)
Partie 2 : Simulation en R
set.seed(120)
# Parametres pour les lois
n <- 30 # Taille de l'echantillon
lambda_values <- seq(0.5, 5, by = 0.1) # Analyse avec differentes valeurs de λ
taille_lambda <- length(lambda_values) # Nombre d'analyses a faire
# Les vecteurs sont initialises avec le meme longueur
bias <- taille_lambda
variance <- taille_lambda
eqm <- taille_lambda
# Simulation pour chaque valeur de λ
for (i in 1:taille_lambda) {
lambda <- lambda_values[i]
estimates <- 1000 # Nombre de simulations
for (j in 1:1000) {
# Simuler le Z et le W
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
# Estimation de lambda
estimates[j] <- sum(W) / sum(Z)
}
# Calcul du biais, variance et EQM
bias[i] <- mean(estimates) - lambda
variance[i] <- var(estimates)
eqm[i] <- bias[i]^2 + variance[i]
}
# Representation graphique
plot(lambda_values, bias, type = "b", col = "blue", pch = 19,
main = "Biais, Variance et EQM de l'estimateur λ̂n",
xlab = "Valeurs de λ", ylab = "Valeurs",
ylim = range(c(bias, variance, eqm)))
# Ajouter la variance
lines(lambda_values, variance, type = "b", col = "green", pch = 19)
# Ajouter l'EQM
lines(lambda_values, eqm, type = "b", col = "purple", pch = 19)
# Ajouter une legende
legend("topright", legend = c("Biais", "Variance", "EQM"),
col = c("blue", "green", "purple"), pch = 19)
Le graphique précédent représente le biais, la variance et l’EQM de l’estimateur \(\hat{\lambda}_n\) pour différentes tailles d’échantillons \(n\).
- Biais (en bleu) : Le biais est proche de zéro pour toutes les valeurs de \(\lambda\), indiquant que l’estimateur est presque non biaisé.
- Variance (en vert) : La variance augmente avec \(\lambda\), donc ça reflète une précision réduite de l’estimation pour des valeurs plus élevées. Cela montre que l’incertitude de l’estimateur augmente pour des paramètres \(\lambda\) plus importants.
- EQM (en violet) : L’EQM suit principalement la courbe de la variance, car le biais est très faible. Cela confirme que la variance domine l’erreur quadratique totale pour cet estimateur.
Question 6 : Calcul de l’information de Fisher et distribution asymptotique de \(\hat{\lambda}_n\)
Partie 1: Theorie
Dans cette analyse, nous calculons l’information de Fisher pour un échantillon \((Z_1, W_1), \dots, (Z_n, W_n)\) et déterminons la distribution asymptotique de l’estimateur \(\hat{\lambda}_n\).
Information de Fisher
L’information de Fisher pour le paramètre \(\lambda\), notée \(I(\lambda)\), est définie comme : \[ I(\lambda) = -\mathbb{E}\left[\frac{\partial^2}{\partial \lambda^2} \log L(\lambda)\right], \] où \(L(\lambda)\) est la fonction de vraisemblance.
Dans ce cas particulier : \[ I(\lambda) = \frac{n}{\lambda^2}, \] car \(Z_i \sim \text{Exp}(\lambda)\) et \(W_i \sim \text{Binomiale}(\lambda / (\lambda + \mu))\).
Distribution asymptotique
Sous des conditions régulières, l’estimateur du maximum de vraisemblance \(\hat{\lambda}_n\) suit asymptotiquement une distribution normale centrée en \(\lambda\) avec une variance égale à l’inverse de l’information de Fisher :
\[ \hat{\lambda}_n \sim \mathcal{N}\left(\lambda, \frac{1}{I(\lambda)}\right). \]
Qui, une fois développé, sera :
\[ \hat{\lambda}_n \sim \mathcal{N}\left(\lambda, \frac{\lambda^2}{n}\right) \]
Partie 2 : Simulation en R
set.seed(120)
# Parametres pour les lois
lambda_true <- 2
mu <- 1
n <- 30
simulations <- 1000 # Nombre de simulations
# Calcul de l'information de Fisher pour le vrai lambda
I_Z <- n / (lambda_true + mu)^2 # Information de Fisher pour Z
p_lambda <- lambda_true / (lambda_true + mu) # Probabilite de W = 1
I_W <- n * (p_lambda * (1 - p_lambda)) # Information de Fisher pour W
fisher_info_true <- I_Z + I_W # Information de Fisher totale
# Variance theorique asymptotique
variance_theorique <- 1 / fisher_info_true
# Simulation pour chaque valeur de λ
for (sim in 1:simulations) {
Z <- rexp(n, rate = lambda_true + mu)
W <- rbinom(n, size = 1, prob = lambda_true / (lambda_true + mu))
# Estimation de lambda
simulations[sim] <- sum(W) / sum(Z)
}
# Calcul de moyenne et variance empiriques
moyenne_empirique <- mean(simulations)
variance_empirique <- var(simulations)
# Representation graphique
hist(simulations, probability = TRUE, breaks = 30, col = "grey",
main = "Distribution de l'estimateur", xlab = expression(hat(lambda)))
curve(dnorm(x, mean = lambda_true, sd = sqrt(variance_theorique)),
col = "red", lwd = 2, add = TRUE)
# Ajouter une legende
legend("topright", legend = c("Empirique", "Théorique"), col = c("grey", "red"), lwd = 2)
Question 7 : Déduction de l’intervalle de confiance asymptotique pour \(\lambda\)
Partie 1 : Théorie
L’estimateur du maximum de vraisemblance pour le paramètre \(\lambda\), noté \(\hat{\lambda}_n\), suit une distribution asymptotique normale centrée autour de la vraie valeur de \(\lambda\), avec une variance donnée par l’inverse de l’information de Fisher. Formellement :
\[ \hat{\lambda}_n \sim \mathcal{N}\left(\lambda, \frac{\lambda^2}{n}\right) \]
Pour un niveau de confiance \(1 - \alpha\), l’intervalle de confiance bilatéral symétrique est donné par :
\[ \left[\hat{\lambda}_n - z_{\alpha/2} \cdot \sqrt{\frac{{\lambda}({\lambda} + \mu)}{n}}, \hat{\lambda}_n + z_{\alpha/2} \cdot \sqrt{\frac{{\lambda}({\lambda} + \mu)}{n}}\right] \]
Partie 2 : Implémentation en R
set.seed(120)
n <- 30
alpha <- 0.05
simulations <- 1000
lambda <- 2
mu <- 1
estimations <- numeric(simulations)
borne_inf <- numeric(simulations)
borne_sup <- numeric(simulations)
for (sim in 1:simulations) {
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
estimations[sim] <- sum(W) / sum(Z)
# Calcul des bornes
sqrt_variance <- sqrt((estimations[sim] * (estimations[sim] + mu)) / n)
borne_inf[sim] <- max(0, estimations[sim] - qnorm(1 - alpha / 2) * sqrt_variance) #Borne inferieur
borne_sup[sim] <- estimations[sim] + qnorm(1 - alpha / 2) * sqrt_variance #Borne superieur
}
cat("L'intervalle de confiance moyen est : [",
mean(borne_inf), ", ", mean(borne_sup), "]\n")## L'intervalle de confiance moyen est : [ 1.172191 , 2.979406 ]mean_lower <- mean(borne_inf)
mean_upper <- mean(borne_sup)
cat("Intervalle de confiance moyen : [", round(mean_lower, 3), ", ", round(mean_upper, 3), "]\n")## Intervalle de confiance moyen : [ 1.172 , 2.979 ]# Pour voir le graphique
library(ggplot2)
ggplot(data.frame(estimations), aes(x = estimations)) +
geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "grey", color = "black") +
geom_vline(xintercept = mean_lower, color = "blue", linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = mean_upper, color = "blue", linetype = "dashed", linewidth = 1) +
labs(title = "Histogramme des estimations de λ avec IC moyen",
x = expression(hat(lambda)[n]),
y = "Densite") +
theme_minimal()
Question 8 : Taux de couverture des intervalles de confiance
Partie 1 : Théorie
le taux de couverture représente la proportion des cas où un intervalle de confiance, construit à partir de différents échantillons, inclut réellement la vraie valeur du paramètre.
Dans cet exercice, nous simulons des données issues de deux composants électroniques, où leurs durées de vie suivent des lois exponentielles. Notre objectif est d’évaluer la performance de l’intervalle de confiance asymptotique que nous avons construit pour le paramètre 𝜆, dans différents scénarios :
En variant le niveau de confiance (𝛼) : 1%, 5%, et 10%;
En modifiant les valeurs de 𝜆: 1, 0,1 et 0,01;
En testant différentes tailles d’échantillons : 10, 20 et 30.
Partie 2 : Implémentation en R
# Parametres
alpha_valeurs <- c(0.01, 0.05, 0.10)
lambda_valeurs <- c(1, 0.1, 0.01)
n_valeurs <- c(10, 20, 30)
mu <- 1
simulations <- 1000
# Calcul du taux de couverture
simuler_taux_couverture <- function(alpha, lambda, n, simulations, mu) {
# Estimations
res <- replicate(simulations, {
Z <- rexp(n, rate = lambda + mu)
W <- rbinom(n, size = 1, prob = lambda / (lambda + mu))
estimations <- sum(W) / sum(Z)
# Calcul des bornes
erreur_type <- sqrt((estimations * (estimations + mu)) / n)
borne_inf <- max(0, estimations - qnorm(1 - alpha / 2) * erreur_type)
borne_sup <- estimations + qnorm(1 - alpha / 2) * erreur_type
# Vérification de la couverture
lambda >= borne_inf && lambda <= borne_sup
})
# Taux de couverture
mean(res)
}
# Application sur les paramètres
resultats <- expand.grid(alpha = alpha_valeurs,
lambda = lambda_valeurs,
n = n_valeurs) %>%
rowwise() %>%
mutate(taux_couverture = simuler_taux_couverture(alpha, lambda, n, simulations, mu)) %>%
ungroup()
print(resultats)## # A tibble: 27 × 4
## alpha lambda n taux_couverture
## <dbl> <dbl> <dbl> <dbl>
## 1 0.01 1 10 0.952
## 2 0.05 1 10 0.93
## 3 0.1 1 10 0.884
## 4 0.01 0.1 10 0.605
## 5 0.05 0.1 10 0.624
## 6 0.1 0.1 10 0.622
## 7 0.01 0.01 10 0.092
## 8 0.05 0.01 10 0.105
## 9 0.1 0.01 10 0.093
## 10 0.01 1 20 0.964
## # ℹ 17 more rowsExercice 2 - La hauteur d’un fleuve.
Question 9 : Estimateur du Maximum de Vraisemblance
Partie 1 : Théorie
Pour \(a > 0\) et des
observations \(x_1, \ldots, x_n \in
\mathbb{R}_+^*\), on a :
\[
L_n(x_1, \ldots, x_n; a) = \prod_{i=1}^n \mathbf{1}{\mathbb{R}+^*}(x_i)
\frac{x_i}{a} \exp\left(-\frac{x_i^2}{2a}\right).
\]
En simplifiant, comme \(x_i > 0\)
pour tout \(i\), on peut écrire :
\[
L_n(x_1, \ldots, x_n; a) = \frac{\prod_{i=1}^n x_i}{a^n}
\exp\left(-\frac{\sum_{i=1}^n x_i^2}{2a}\right).
\]
En passant au logarithme, on obtient :
\[
\log L_n(x_1, \ldots, x_n; a) = \sum_{i=1}^n \log x_i - n \log a -
\frac{\sum_{i=1}^n x_i^2}{2a}.
\]
En dérivant par rapport à \(a\), on
obtient :
\[
\frac{\partial}{\partial a} \log L_n(x_1, \ldots, x_n; a) = -\frac{n}{a}
+ \frac{\sum_{i=1}^n x_i^2}{2a^2}.
\]
En résolvant \(\frac{\partial}{\partial a}
\log L_n(x_1, \ldots, x_n; a) = 0\), on trouve :
\[
\frac{n}{a} = \frac{\sum_{i=1}^n x_i^2}{2a^2}.
\]
En isolant \(a\), on obtient :
\[
\hat{a}_{\text{MV}} = \frac{\sum_{i=1}^n x_i^2}{2n}.
\]
Ainsi, l’estimateur du maximum de vraisemblance pour \(a\) est donné par :
\[
\hat{a}_{\text{MV}} = \frac{\sum_{i=1}^n x_i^2}{2n}
\]
Partie 2 : Simulation en R
set.seed(120)
true_a <- 1 # Valeur réelle du paramètre a
n_sim <- 1000 # Nombre de simulations
n_sample <- 100 # Taille de l'échantillon observé
# Fonction pour générer des échantillons Rayleigh et calculer l'estimateur a_hat
simulate_a_hat <- function(a, n, n_sim) {
a_hats <- numeric(n_sim) # Stockage des résultats
for (i in 1:n_sim) {
sample <- sqrt(-2 * a * log(runif(n))) # Génération de la loi de Rayleigh
a_hats[i] <- (1 / (2 * n)) * sum(sample^2) # Estimation de a_hat
}
return(a_hats)
}
# Générer les estimations simulées
simulated_a_hats <- simulate_a_hat(true_a, n_sample, n_sim)
# Calculer le biais, la variance et la convergence
mean_a_hat <- mean(simulated_a_hats) # Moyenne des estimations
bias_a_hat <- mean_a_hat - true_a # Calcul du biais
variance_a_hat <- var(simulated_a_hats) # Variance des estimations
variance_a_hat## [1] 0.01029914bias_a_hat## [1] 0.00122328mean_a_hat## [1] 1.001223Question 10 : L’estimateur du maximum de vraisemblance est-il sans biais ?
Un estimateur est considéré comme sans biais si son espérance est égale à la vraie valeur du paramètre estimé, soit : \[ \mathbb{E}[\hat{a}_n] = a. \]
Calcul de l’espérance de \(\hat{a}_n\)
L’estimateur du maximum de vraisemblance de \(a\) est donné par : \[ \hat{a}n = \frac{1}{2n} \sum_{i=1}^n X_i^2, \] où \(X_i\) suit une loi de Rayleigh avec le paramètre \(a\).
L’espérance de \(\hat{a}_n\) s’écrit : \[ \mathbb{E}[\hat{a}n] = \mathbb{E}\left[\frac{1}{2n} \sum_{i=1}^n X_i^2\right]. \]
Puisque les \(X_i\) sont indépendants et identiquement distribués, on obtient : \[ \mathbb{E}[\hat{a}_n] = \frac{1}{2n} \cdot n \cdot \mathbb{E}[X_1^2]. \]
Pour une variable \(X_i\) suivant une loi de Rayleigh, l’espérance de \(X_i^2\) est connue : \[ \mathbb{E}[X_i^2] = 2a. \]
Ainsi : \[ \mathbb{E}[\hat{a}_n] = \frac{1}{2n} \cdot n \cdot 2a = a. \]
Conclusion
L’estimateur \(\hat{a}_n\) est donc sans biais, car : \[ \mathbb{biais} = a - a = 0. \]
cat("Biais de a_hat=:", bias_a_hat)## Biais de a_hat=: 0.00122328Le faible biais est conforme au fait que l’espérance est sans biais.
Question 11 : L’estimateur du maximum de vraisemblance est-il convergent ?
Partie 1 : Convergence en probabilité
Un estimateur \(\hat{a}_n\) est dit convergent en probabilité s’il converge vers le vrai paramètre \(a\), c’est-à-dire : \[ \hat{a}_n \xrightarrow{} a \quad \text{lorsque } n \to \infty. \]
L’estimateur \(\hat{a}_n\) est donné par : \[ \hat{a}n = \frac{1}{2n} \sum_{i=1}^n X_i^2, \] où les \(X_i\) sont indépendants et identiquement distribués suivant une loi de Rayleigh.
D’après la loi forte des grands nombres, on a : \[ \frac{1}{n} \sum_{i=1}^n X_i^2 \xrightarrow{\text{}} \mathbb{E}[X_i^2]. \]
Pour une loi de Rayleigh de paramètre \(a\), \(\mathbb{E}[X_i^2] = 2a\). Ainsi : \[ \frac{1}{n} \sum_{i=1}^n X_i^2 \xrightarrow{\text{}} 2a. \]
En divisant par 2, on obtient : \[ \hat{a}_n \xrightarrow{\text{}} a. \]
Conclusion : \(\hat{a}_n\) est convergent en probabilité.
Partie 2 : Absolument convergent
Pour vérifier si \(\hat{a}_n\) est absolument convergent, nous devons montrer qu’il converge presque sûrement vers \(a\), soit : \[ \hat{a}_n \xrightarrow{\text{}} a \]
En utilisant la loi forte des grands nombres, on sait que : \[ \frac{1}{n} \sum_{i=1}^n X_i^2 \xrightarrow{\text{}} \mathbb{E}[X_i^2] = 2a \]
En divisant par 2, on obtient : \[ \hat{a}n = \frac{1}{2n} \sum_{i=1}^n X_i^2 \xrightarrow{\text{}} a \]
Conclusion : \(\hat{a}_n\) est absolument convergent.
cat("Variance des estimations :", variance_a_hat)## Variance des estimations : 0.01029914La variance diminue avec n, conformément à la convergence.
Question 12: Information de Fisher
Partie 1 : Théorie
- L’information de Fisher pour un échantillon de taille \(n\) est donnée par
-
\[ I(a) = \frac{n}{a^2}. \] La variance asymptotique de l’estimateur est : \[ \text{Var}(\hat{a}_n) = \frac{1}{I(a)} = \frac{a^2}{n}. \]
Partie 2 : Simulation en R
# Calcul de l'information de Fisher
information_fisher <- function(a, n) {
return(n / (a^2)) # Formule
}
# Calcul de l'information de Fisher pour les données simulées
fisher_info <- information_fisher(true_a, n)
variance_asymptotique <- 1 / fisher_info # Variance asymptotique Résultats de la simulation :
## Information de Fisher : 30## Variance asymptotique de l'estimateur : 0.03333333Question 13: Intervalle de Confiance
Partie 1 : Théorie
L’intervalle de confiance asymptotique à un niveau de confiance \(1-\alpha\) est : \[ IC = \left[ \hat{a}_n - z_{1-\alpha/2} \sqrt{\text{Var}(\hat{a}_n)}, \; \hat{a}_n + z_{1-\alpha/2} \sqrt{\text{Var}(\hat{a}_n)} \right]. \]
Partie 2 : Calcul de l’Intervalle en R
alpha <- 0.05 # Niveau de confiance (95%)
z_alpha <- qnorm(1 - alpha / 2) # Quantile de la loi normale
# Calcul de l'intervalle de confiance
lower_bound <- mean_a_hat - z_alpha * sqrt(variance_asymptotique / n)
upper_bound <- mean_a_hat + z_alpha * sqrt(variance_asymptotique / n)Résultats de la simulation :
## Intervalle de confiance asymptotique à 95% : [ 0.9358911 , 1.066555 ]Question 14: Simulation pour vérifier le taux de couverture
Partie 1 : Fonction de Simulation
coverage_simulation <- function(true_a, n, n_sim, alpha) {
covered <- 0
for (i in 1:n_sim) {
sample <- sqrt(-2 * true_a * log(runif(n))) # Generation de Rayleigh
a_hat <- (1 / (2 * n)) * sum(sample^2) # Calcul de a_hat
fisher_info <- n / (2 * a_hat^2) # Information de Fisher pour l'estimation
variance_asymptotique <- 1 / fisher_info # Variance asymptotique
z_alpha <- qnorm(1 - alpha / 2)
lower_bound <- a_hat - z_alpha * sqrt(variance_asymptotique) # Borne inferieure
upper_bound <- a_hat + z_alpha * sqrt(variance_asymptotique) # Borne superieure
if (true_a >= lower_bound && true_a <= upper_bound) {
covered <- covered + 1
}
}
return(covered / n_sim) # Taux de couverture
}Partie 2 : Simulation pour différents niveaux de alpha
# Parametres
true_a <- 2 # Valeur du paramètre a
n_sample <- 50 # Taille de l'échantillon
n_sim <- 1000 # Nombre de simulations
# Calcul des taux de couverture
coverage_1 <- coverage_simulation(true_a, n_sample, n_sim, 0.01) # alpha = 1%
coverage_5 <- coverage_simulation(true_a, n_sample, n_sim, 0.05) # alpha = 5%
coverage_10 <- coverage_simulation(true_a, n_sample, n_sim, 0.10) # alpha = 10%## Taux de couverture pour alpha = 1% : 0.997## Taux de couverture pour alpha = 5% : 0.983## Taux de couverture pour alpha = 10% : 0.973Les résultats montrent que l’intervalle de confiance fonctionne très bien, avec des taux de couverture proches de 1.
Question 15: Analyse de la Probabilité de Dépassement
Partie 1 : Théorie
La probabilité qu’une hauteur maximale du fleuve dépasse 6 mètres est
donnée par :
\[
P(X > 6) = 1 - F(6)
\]
où \(F(x)\) est la fonction de
répartition de la loi de Rayleigh.
La fonction de répartition est définie par :
\[
F(x) = 1 - \exp\left(-\frac{x^2}{2a}\right)
\]
Cela implique que :
\[
P(X > 6) = \exp\left(-\frac{6^2}{2a}\right)
\]
Pour une période de 8 ans, la probabilité qu’au moins une hauteur
dépasse 6 mètres est donnée par :
\[
P_{\text{au moins une fois}} = 1 - \left(1 - P(X > 6)\right)^8
\]
Ce calcul repose sur l’indépendance des observations chaque année.
Partie 2 : Implémentation dans R
# Données observées
hauteurs <- c(2.27, 3.75, 0.29, 6.88, 5.17, 4.29, 4.07, 4.47)
# Estimation de a par maximum de vraisemblance
n <- length(hauteurs) # Taille de l'échantillon
a_estime <- sum(hauteurs^2) / (2 * n) # Estimation de a## Estimation de a par maximum de vraisemblance : 9.269544Partie 3 : Probabilité de dépassement du seuil de 6 mètres
# Calcul de P(X > 6)
seuil <- 6
P_X_sup_6 <- exp(-seuil^2 / (2 * a_estime))## Probabilité : 0.1434393Partie 4 : Probabilité d’au moins une crue en 8 ans
# Probabilité d'au moins 8 ans
annees <- 8
P_au_moins_une_fois <- 1 - (1 - P_X_sup_6)^annees## Probabilité : 0.7102222L’estimation de l’entreprise d’assurance n’est pas justifiée, car la probabilité d’une crue de 6 mètres est de 14,34% par an, bien plus élevée que les 0,1% estimés par l’entreprise.
Le niveau de risque est élevé, avec une probabilité de 71,02% qu’une crue de 6 mètres se produise au moins une fois en 8 ans.