Exercise 1 - Les deux composants
électroniques.
Question 1. Que signifient les variables aléatoires
Zi et Wi ?
# Paramètres
lambda <- 2
mu <- 3
n <- 1000 # Taille de l'échantillon
# Simulation
set.seed(123)
X <- rexp(n, rate = lambda) # Composant X
Y <- rexp(n, rate = mu) # Composant Y
# Calcul des variables Z et W
Z <- pmin(X, Y) # Temps jusqu'à la défaillance (minimum entre X et Y)
W <- ifelse(Z == X, 1, 0) # Identifie quel composant a échoué en premier
# Résultats de base
cat("Premières valeurs de Z (temps jusqu'à la défaillance) :\n")
## Premières valeurs de Z (temps jusqu'à la défaillance) :
print(head(Z))
## [1] 0.07449894 0.03920849 0.66452743 0.01578868 0.02810549 0.10433064
cat("\nProportion de W (échec du premier composant) :\n")
##
## Proportion de W (échec du premier composant) :
print(table(W) / n)
## W
## 0 1
## 0.607 0.393
Question 2. Donner les lois des deux variables
aléatoires Zi et Wi. Vérifier vos affirmations avec R.
# Paramètres
lambda <- 2 # Taux du composant X
mu <- 3 # Taux du composant Y
n <- 10000 # Taille de l'échantillon
# Simulation
set.seed(123)
X <- rexp(n, rate = lambda) # Temps de défaillance du composant X
Y <- rexp(n, rate = mu) # Temps de défaillance du composant Y
# Calcul de Z et W
Z <- pmin(X, Y) # Temps minimum jusqu'à la défaillance
W <- ifelse(Z == X, 1, 0) # Quel composant a échoué en premier
# Validation de la distribution de Z
dens_Z <- density(Z)
plot(dens_Z, main = "Densité simulée de Z", xlab = "Z", ylab = "Densité")
curve((lambda + mu) * exp(-(lambda + mu) * x), add = TRUE, col = "red", lty = 2)
legend("topright", legend = c("Simulée", "Théorique"), col = c("black", "red"), lty = 1:2)
# Validation de la distribution de W
prop_W <- table(W) / n
cat("\nProportion simulée de W :\n")
##
## Proportion simulée de W :
print(prop_W)
## W
## 0 1
## 0.5984 0.4016
# Probabilités théoriques
cat("\nProbabilités théoriques de W :\n")
##
## Probabilités théoriques de W :
cat("P(W = 1) :", lambda / (lambda + mu), "\n")
## P(W = 1) : 0.4
cat("P(W = 0) :", mu / (lambda + mu), "\n")
## P(W = 0) : 0.6
Question 3. Déterminer l’estimateur du maximum de
vraisemblance λˆn de λ. Vérifier votre résultat avec R.
# Paramètres connus
mu <- 1 # Paramètre du second composant
lambda_true <- 2 # Valeur réelle de lambda pour la simulation
n <- 1000 # Taille de l'échantillon
# Simulation
set.seed(123)
X <- rexp(n, rate = lambda_true)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
# Estimateur de maximum de vraisemblance
lambda_hat <- sum(W) / sum(Z)
cat("Estimateur de maximum de vraisemblance (lambda) :", lambda_hat, "\n")
## Estimateur de maximum de vraisemblance (lambda) : 1.900831
# Comparaison avec la valeur réelle
cat("Valeur réelle de lambda :", lambda_true, "\n")
## Valeur réelle de lambda : 2
Question 4. À l’aide de R, représenter, pour une
valeur de λ=2, le biais, la variance et l’EQM de l’estimateur du maximum
de vraisemblance λˆn, en fonction de la taille de l’échantillon
n. L’estimateur semble-t-il convergent ? L’estimateur semble-t-il
absolument convergent ?
# Paramètres
lambda <- 2 # Valeur réelle de lambda
mu <- 3 # Paramètre fixe du second composant
n_sim <- 100 # Nombre de simulations par taille d'échantillon
n_values <- seq(10, 1000, by = 10) # Tailles d'échantillon
# Fonctions auxiliaires pour le calcul
simulate_lambda_hat <- function(n, lambda, mu) {
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
sum(W) / sum(Z)
}
# Résultats stockés
bias <- numeric(length(n_values))
variance <- numeric(length(n_values))
mse <- numeric(length(n_values))
means <- numeric(length(n_values))
# Simulations
set.seed(123)
for (i in seq_along(n_values)) {
n <- n_values[i]
estimates <- replicate(n_sim, simulate_lambda_hat(n, lambda, mu))
means[i] <- mean(estimates)
bias[i] <- mean(estimates) - lambda
variance[i] <- var(estimates)
mse[i] <- variance[i] + bias[i]^2
}
# Graphiques
par(mfrow = c(3, 1)) # Trois graphiques dans une seule fenêtre
# Tendance
plot(n_values, means, type = "l", col = "blue", lwd = 2,
xlab = "Taille de l'échantillon (n)", ylab = "Tendance",
main = "Tendance de l'estimateur")
abline(h = lambda, col = "red", lty = 2) # Valeur réelle de lambda
legend("topright", legend = c("Estimation", "Valeur réelle"),
col = c("blue", "red"), lty = c(1, 2))
# Variance
plot(n_values, variance, type = "l", col = "green", lwd = 2,
xlab = "Taille de l'échantillon (n)", ylab = "Variance",
main = "Variance de l'estimateur")
# EQM
plot(n_values, mse, type = "l", col = "purple", lwd = 2,
xlab = "Taille de l'échantillon (n)", ylab = "Erreur Quadratique Moyenne (EQM)",
main = "Erreur Quadratique Moyenne de l'estimateur")
Question 5. À l’aide de R, représenter, pour une
valeur de n=30, le biais, la variance et l’EQM de l’estimateur du
maximum de vraisemblance λˆn, en fonction du paramètre λ.
# Paramètres fixes
mu <- 3 # Paramètre fixe du composant Y
n <- 30 # Taille de l'échantillon
n_sim <- 100 # Nombre de simulations par valeur de lambda
lambda_values <- seq(0.5, 5, by = 0.1) # Valeurs de lambda à tester
# Fonctions auxiliaires
simulate_lambda_hat <- function(lambda, mu, n) {
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
sum(W) / sum(Z)
}
# Résultats stockés
means <- numeric(length(lambda_values))
variance <- numeric(length(lambda_values))
rmse <- numeric(length(lambda_values))
# Simulations
set.seed(123)
for (i in seq_along(lambda_values)) {
lambda <- lambda_values[i]
estimates <- replicate(n_sim, simulate_lambda_hat(lambda, mu, n))
means[i] <- mean(estimates)
variance[i] <- var(estimates)
bias <- means[i] - lambda
rmse[i] <- sqrt(variance[i] + bias^2)
}
# Graphiques
par(mfrow = c(3, 1)) # Trois graphiques dans une seule fenêtre
# Tendance
plot(lambda_values, means, type = "l", col = "blue", lwd = 2,
xlab = expression(lambda), ylab = "Tendance",
main = "Tendance de l'estimateur en fonction de λ")
abline(a = 0, b = 1, col = "red", lty = 2) # Ligne de référence (y = x)
legend("topright", legend = c("Estimation", "Ligne de référence"),
col = c("blue", "red"), lty = c(1, 2))
# Variance
plot(lambda_values, variance, type = "l", col = "green", lwd = 2,
xlab = expression(lambda), ylab = "Variance",
main = "Variance de l'estimateur en fonction de λ")
# RMSE
plot(lambda_values, rmse, type = "l", col = "purple", lwd = 2,
xlab = expression(lambda), ylab = "RMSE",
main = "Erreur Quadratique Moyenne (RMSE) en fonction de λ")
Question 6. Calculer l’information de Fisher
contenue dans l’échantillon (Z1,W1),…,(Zn,Wn) pour le paramètre λ. En
utilisant le résultat du cours, donner la distribution asymptotique de
l’estimateur λˆn . Vérifier votre résultat avec R.
# Paramètres
lambda <- 2 # Valeur réelle de lambda
mu <- 3 # Paramètre fixe de mu
n <- 30 # Taille de l'échantillon
n_sim <- 1000 # Nombre de simulations
# Fonction pour calculer l'information de Fisher
fisher_information <- function(lambda, mu, n) {
term1 <- n / (lambda + mu)^2
term2 <- n / (lambda * (lambda + mu))
term1 + term2
}
# Calcul de l'information de Fisher
I_n <- fisher_information(lambda, mu, n)
cat("Information de Fisher :", I_n, "\n")
## Information de Fisher : 4.2
# Variance asymptotique
var_assintotica <- 1 / I_n
cat("Variance asymptotique :", var_assintotica, "\n")
## Variance asymptotique : 0.2380952
# Simuler l'estimateur de lambda
simulate_lambda_hat <- function(n, lambda, mu) {
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
sum(W) / sum(Z)
}
set.seed(123)
lambda_hats <- replicate(n_sim, simulate_lambda_hat(n, lambda, mu))
# Comparaison avec la distribution normale
hist(lambda_hats, breaks = 30, probability = TRUE, main = "Distribution de λ̂n",
xlab = expression(hat(lambda)), col = "lightblue")
curve(dnorm(x, mean = lambda, sd = sqrt(var_assintotica)), col = "red", lwd = 2, add = TRUE)
legend("topright", legend = c("Simulée", "Théorique (Normale)"),
col = c("lightblue", "red"), lty = c(1, 1), lwd = c(1, 2))
Question 7. Déduire de la distribution asymptotique
obtenue à la question 6. précédente un intervalle de confiance bilatéral
symétrique asymptotique de niveau 1−α pour le paramètre λ. Vérifier
votre résultat avec R.
# Paramètres
lambda <- 2 # Valeur réelle de lambda
mu <- 3 # Paramètre fixe de mu
n <- 30 # Taille de l'échantillon
alpha <- 0.05 # Niveau de signification
n_sim <- 1000 # Nombre de simulations
# Fonction pour calculer l'information de Fisher
fisher_information <- function(lambda, mu, n) {
term1 <- n / (lambda + mu)^2
term2 <- n / (lambda * (lambda + mu))
term1 + term2
}
# Information de Fisher et variance asymptotique
I_n <- fisher_information(lambda, mu, n)
var_assintotica <- 1 / I_n
z_alpha <- qnorm(1 - alpha / 2)
# Simuler l'estimateur de lambda
simulate_lambda_hat <- function(n, lambda, mu) {
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
sum(W) / sum(Z)
}
set.seed(123)
lambda_hats <- replicate(n_sim, simulate_lambda_hat(n, lambda, mu))
# Calcul des intervalles de confiance
lower_bounds <- lambda_hats - z_alpha * sqrt(var_assintotica)
upper_bounds <- lambda_hats + z_alpha * sqrt(var_assintotica)
# Proportion de couverture de l'intervalle
coverage <- mean(lower_bounds <= lambda & upper_bounds >= lambda)
# Résultats
cat("Information de Fisher :", I_n, "\n")
## Information de Fisher : 4.2
cat("Variance asymptotique :", var_assintotica, "\n")
## Variance asymptotique : 0.2380952
cat("Proportion de couverture de l'intervalle de confiance :", coverage, "\n")
## Proportion de couverture de l'intervalle de confiance : 0.9
# Graphique des intervalles de confiance
plot(lambda_hats, ylim = c(min(lower_bounds), max(upper_bounds)), pch = 16,
col = "blue", xlab = "Simulations", ylab = "Estimation et Intervalle")
arrows(1:n_sim, lower_bounds, 1:n_sim, upper_bounds, angle = 90, code = 3, length = 0.05, col = "blue")
abline(h = lambda, col = "red", lwd = 2, lty = 2) # Valeur réelle de lambda
legend("topright", legend = c("Intervalles", "Valeur réelle"),
col = c("blue", "red"), lty = c(1, 2), lwd = c(1, 2))
Question 8. Le taux de couverture est la
probabilité avec laquelle un intervalle de confiance contient la vraie
valeur inconnue du paramètre. Vous calculerez, par simulation avec R, le
taux de couverture réel de l’intervalle de confiance construit à la
question 7. pour des valeurs valeurs de α=1%, α=5% et α=10%, des valeurs
de λ=1, λ=0,1 et λ=0,01 et des effectifs d’échantillons n=10 n=20 et
n=30.
# Paramètres généraux
alpha_values <- c(0.01, 0.05, 0.10) # Niveaux de signification
lambda_values <- c(1, 0.1, 0.01) # Valeurs de lambda
n_values <- c(10, 20, 30) # Tailles d'échantillon
mu <- 3 # Valeur fixe pour mu
n_sim <- 1000 # Nombre de simulations
# Fonction pour calculer l'information de Fisher
fisher_information <- function(lambda, mu, n) {
term1 <- n / (lambda + mu)^2
term2 <- n / (lambda * (lambda + mu))
term1 + term2
}
# Fonction pour simuler et calculer le taux de couverture
simulate_coverage <- function(n, lambda, mu, alpha, n_sim) {
# Information de Fisher et variance
I_n <- fisher_information(lambda, mu, n)
var_assintotica <- 1 / I_n
z_alpha <- qnorm(1 - alpha / 2)
# Simuler l'estimateur
lambda_hats <- replicate(n_sim, {
X <- rexp(n, rate = lambda)
Y <- rexp(n, rate = mu)
Z <- pmin(X, Y)
W <- ifelse(Z == X, 1, 0)
sum(W) / sum(Z)
})
# Calculer les intervalles de confiance
lower_bounds <- lambda_hats - z_alpha * sqrt(var_assintotica)
upper_bounds <- lambda_hats + z_alpha * sqrt(var_assintotica)
# Taux de couverture
mean(lower_bounds <= lambda & upper_bounds >= lambda)
}
# Réaliser les simulations
results <- expand.grid(Alpha = alpha_values, Lambda = lambda_values, N = n_values)
results$Coverage <- NA
for (i in 1:nrow(results)) {
alpha <- results$Alpha[i]
lambda <- results$Lambda[i]
n <- results$N[i]
# Calculer le taux de couverture
results$Coverage[i] <- simulate_coverage(n, lambda, mu, alpha, n_sim)
}
# Afficher les résultats
print(results)
## Alpha Lambda N Coverage
## 1 0.01 1.00 10 0.946
## 2 0.05 1.00 10 0.905
## 3 0.10 1.00 10 0.809
## 4 0.01 0.10 10 0.956
## 5 0.05 0.10 10 0.940
## 6 0.10 0.10 10 0.895
## 7 0.01 0.01 10 0.969
## 8 0.05 0.01 10 0.968
## 9 0.10 0.01 10 0.964
## 10 0.01 1.00 20 0.952
## 11 0.05 1.00 20 0.902
## 12 0.10 1.00 20 0.844
## 13 0.01 0.10 20 0.975
## 14 0.05 0.10 20 0.940
## 15 0.10 0.10 20 0.896
## 16 0.01 0.01 20 0.937
## 17 0.05 0.01 20 0.933
## 18 0.10 0.01 20 0.920
## 19 0.01 1.00 30 0.965
## 20 0.05 1.00 30 0.925
## 21 0.10 1.00 30 0.862
## 22 0.01 0.10 30 0.975
## 23 0.05 0.10 30 0.952
## 24 0.10 0.10 30 0.926
## 25 0.01 0.01 30 0.935
## 26 0.05 0.01 30 0.906
## 27 0.10 0.01 30 0.914
# Graphiques
par(mfrow = c(1, 3)) # Diviser en trois graphiques
for (alpha in alpha_values) {
subset_data <- subset(results, Alpha == alpha)
plot(subset_data$Lambda, subset_data$Coverage, type = "b", pch = 16,
xlab = expression(lambda), ylab = "Taux de couverture",
main = paste("Alpha =", alpha), ylim = c(0, 1),
col = as.factor(subset_data$N))
legend("bottomright", legend = paste("N =", n_values),
col = 1:length(n_values), lty = 1, pch = 16)
}
Question 9. Déterminer l’estimateur du maximum de
vraisemblance aˆn du paramètre réel a > 0.
# Paramètres
set.seed(123)
a_true <- 3 # Valeur réelle de a
n <- 1000 # Taille de l'échantillon
simulations <- 1000 # Nombre de simulations
# Fonction pour simuler des échantillons de Rayleigh
simulate_rayleigh <- function(n, a) {
sqrt(-2 * a * log(runif(n)))
}
# Fonction pour calculer l'EMV
estimate_a <- function(x) {
sum(x^2) / (2 * length(x))
}
# Simuler et calculer l'EMV
mle_values <- replicate(simulations, {
x <- simulate_rayleigh(n, a_true)
estimate_a(x)
})
# Résultats
cat("Moyenne de l'EMV :", mean(mle_values), "\n")
## Moyenne de l'EMV : 3.005436
cat("Valeur réelle de a :", a_true, "\n")
## Valeur réelle de a : 3
cat("Variance de l'EMV :", var(mle_values), "\n")
## Variance de l'EMV : 0.008952662
# Graphique de convergence
hist(mle_values, probability = TRUE, main = "Distribution de l'EMV pour a",
xlab = expression(hat(a)), col = "lightblue", border = "white")
abline(v = a_true, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Valeur réelle de a"), col = "red", lty = 2, lwd = 2)
Question 10. L’estimateur du maximum de
vraisemblance aˆn est-il sans biais ?
# Paramètres
set.seed(123)
a_true <- 3 # Valeur réelle de a
n <- 1000 # Taille de l'échantillon
simulations <- 1000 # Nombre de simulations
# Fonction pour simuler la distribution de Rayleigh
simulate_rayleigh <- function(n, a) {
sqrt(-2 * a * log(runif(n)))
}
# Fonction pour calculer l'estimateur de maximum de vraisemblance (MLE)
estimate_a <- function(x) {
sum(x^2) / (2 * length(x))
}
# Simuler et calculer le MLE
mle_values <- replicate(simulations, {
x <- simulate_rayleigh(n, a_true)
estimate_a(x)
})
# Vérifier les propriétés
mean_mle <- mean(mle_values) # Moyenne de l'EMV
var_mle <- var(mle_values) # Variance de l'EMV
cat("Moyenne de l'EMV :", mean_mle, "\n")
## Moyenne de l'EMV : 3.005436
cat("Variance de l'EMV :", var_mle, "\n")
## Variance de l'EMV : 0.008952662
cat("Valeur réelle de a :", a_true, "\n")
## Valeur réelle de a : 3
# Graphique de convergence
hist(mle_values, probability = TRUE, main = "Distribution de l'EMV pour a",
xlab = expression(hat(a)), col = "lightblue", border = "white")
abline(v = a_true, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Valeur réelle de a"), col = "red", lty = 2, lwd = 2)
Question 11. L’estimateur du maximum de
vraisemblance aˆn est-il convergent ? Est-il absolument convergent
?
# Paramètres
set.seed(123)
a_true <- 3 # Valeur réelle de a
n <- c(10, 50, 100, 500) # Tailles d'échantillon
simulations <- 1000 # Nombre de simulations
# Fonction pour simuler la distribution de Rayleigh
simulate_rayleigh <- function(n, a) {
sqrt(-2 * a * log(runif(n)))
}
# Fonction pour calculer le MLE
estimate_a <- function(x) {
sum(x^2) / (2 * length(x))
}
# Simuler et vérifier la convergence
results <- data.frame(SampleSize = integer(), Bias = numeric(), Variance = numeric(), MSE = numeric())
for (sample_size in n) {
mle_values <- replicate(simulations, {
x <- simulate_rayleigh(sample_size, a_true)
estimate_a(x)
})
bias <- mean(mle_values) - a_true
variance <- var(mle_values)
mse <- bias^2 + variance
results <- rbind(results, data.frame(SampleSize = sample_size, Bias = bias, Variance = variance, MSE = mse))
}
# Afficher les résultats
print(results)
## SampleSize Bias Variance MSE
## 1 10 0.012168787 0.85043728 0.85058536
## 2 50 0.009022152 0.18584951 0.18593091
## 3 100 -0.003439051 0.09283169 0.09284352
## 4 500 0.004711651 0.01776622 0.01778842
# Graphique de convergence de l'estimateur
plot(n, results$MSE, type = "b", col = "blue", lwd = 2,
xlab = "Taille de l'échantillon (n)", ylab = "Erreur Quadratique Moyenne (EQM)",
main = "Convergence de l'estimateur de maximum de vraisemblance")
abline(h = 0, col = "red", lty = 2)
legend("topright", legend = c("EQM", "Valeur réelle"), col = c("blue", "red"), lty = c(1, 2))
Question 12. Calculer l’information de Fisher
contenue dans l’échantillon (X1,…,Xn) pour le paramètre a. En utilisant
le résultat du cours, donner la distribution asymptotique de
l’estimateur aˆn. Vérifier votre résultat avec R.
# Paramètres
set.seed(123)
a_true <- 3 # Valeur réelle de a
n <- 50 # Taille de l'échantillon
simulations <- 1000 # Nombre de simulations
# Fonction pour calculer l'information de Fisher
fisher_information <- function(a, n) {
2 * n / a^2
}
# Calcul de l'information de Fisher
I_n <- fisher_information(a_true, n)
cat("Information de Fisher :", I_n, "\n")
## Information de Fisher : 11.11111
# Simulations pour valider la distribution asymptotique
simulate_mle <- function(n, a) {
x <- sqrt(-2 * a * log(runif(n)))
sum(x^2) / (2 * n)
}
mle_values <- replicate(simulations, simulate_mle(n, a_true))
# Comparaison avec la distribution normale
mean_mle <- mean(mle_values)
var_mle <- var(mle_values)
theoretical_var <- 1 / I_n
cat("Moyenne de l'EMV :", mean_mle, "\n")
## Moyenne de l'EMV : 3.017395
cat("Variance de l'EMV (simulée) :", var_mle, "\n")
## Variance de l'EMV (simulée) : 0.1771729
cat("Variance théorique :", theoretical_var, "\n")
## Variance théorique : 0.09
# Graphique pour comparer la distribution simulée et théorique
hist(mle_values, probability = TRUE, main = "Distribution de l'EMV pour a",
xlab = expression(hat(a)), col = "lightblue", border = "white")
curve(dnorm(x, mean = a_true, sd = sqrt(theoretical_var)),
add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("Distribution simulée", "Distribution théorique"),
col = c("lightblue", "red"), lty = c(1, 1), lwd = c(2, 2))
Question 13. Déduire de la distribution
asymptotique précédente un intervalle de confiance bilatéral symétrique
asymptotique de niveau 1−α pour le paramètre a. Vérifier votre résultat
avec R.
# Paramètres
set.seed(123)
a_true <- 3 # Valeur réelle de a
n <- 50 # Taille de l'échantillon
alpha <- 0.05 # Niveau de signification
z_alpha <- qnorm(1 - alpha / 2) # Quantile de la normale standard
# Fonction pour calculer l'information de Fisher
fisher_information <- function(a, n) {
2 * n / a^2
}
# Calculer l'information de Fisher et la variance
I_n <- fisher_information(a_true, n)
var_a_hat <- 1 / I_n
# Simuler l'estimateur a_hat
simulate_a_hat <- function(n, a) {
x <- sqrt(-2 * a * log(runif(n)))
sum(x^2) / (2 * n)
}
set.seed(123)
a_hats <- replicate(1000, simulate_a_hat(n, a_true))
# Calcul des intervalles de confiance
lower_bounds <- a_hats - z_alpha * sqrt(var_a_hat)
upper_bounds <- a_hats + z_alpha * sqrt(var_a_hat)
# Proportion de couverture
coverage <- mean(lower_bounds <= a_true & upper_bounds >= a_true)
# Résultats
cat("Information de Fisher :", I_n, "\n")
## Information de Fisher : 11.11111
cat("Variance asymptotique :", var_a_hat, "\n")
## Variance asymptotique : 0.09
cat("Proportion de couverture :", coverage, "\n")
## Proportion de couverture : 0.827
# Graphique des intervalles de confiance
hist(a_hats, probability = TRUE, main = "Distribution de l'EMV et intervalles de confiance",
xlab = expression(hat(a)), col = "lightblue", border = "white")
abline(v = a_true, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Valeur réelle de a"), col = "red", lty = 2, lwd = 2)
Question 14. Le taux de couverture est la
probabilité avec laquelle un intervalle de confiance contient la vraie
valeur inconnue du paramètre. Calculer, par simulation, le taux de
couverture réel de l’intervalle de confiance construit à la question 13.
précédente pour des valeurs de α=1%, α=5% et α=10%, des valeurs de a=10,
a=1 et a=0,1 et des effectifs d’échantillons n=8 n=16 et n=32.
# Paramètres généraux
alpha_values <- c(0.01, 0.05, 0.10) # Niveaux de signification
a_values <- c(10, 1, 0.1) # Valeurs de a
n_values <- c(8, 16, 32) # Tailles d'échantillon
n_sim <- 1000 # Nombre de simulations
# Fonction pour calculer l'information de Fisher
fisher_information <- function(a, n) {
2 * n / a^2
}
# Simuler et calculer le taux de couverture
simulate_coverage <- function(n, a, alpha, n_sim) {
# Information de Fisher
I_n <- fisher_information(a, n)
var_a_hat <- 1 / I_n
z_alpha <- qnorm(1 - alpha / 2)
# Simuler l'estimateur
a_hats <- replicate(n_sim, {
x <- sqrt(-2 * a * log(runif(n)))
sum(x^2) / (2 * n)
})
# Calculer les intervalles de confiance
lower_bounds <- a_hats - z_alpha * sqrt(var_a_hat)
upper_bounds <- a_hats + z_alpha * sqrt(var_a_hat)
# Taux de couverture
mean(lower_bounds <= a & upper_bounds >= a)
}
# Résultats
results <- expand.grid(Alpha = alpha_values, A = a_values, N = n_values)
results$Coverage <- NA
for (i in 1:nrow(results)) {
alpha <- results$Alpha[i]
a <- results$A[i]
n <- results$N[i]
# Calculer le taux de couverture
results$Coverage[i] <- simulate_coverage(n, a, alpha, n_sim)
}
# Afficher les résultats
print(results)
## Alpha A N Coverage
## 1 0.01 10.0 8 0.942
## 2 0.05 10.0 8 0.866
## 3 0.10 10.0 8 0.773
## 4 0.01 1.0 8 0.935
## 5 0.05 1.0 8 0.870
## 6 0.10 1.0 8 0.783
## 7 0.01 0.1 8 0.952
## 8 0.05 0.1 8 0.852
## 9 0.10 0.1 8 0.746
## 10 0.01 10.0 16 0.934
## 11 0.05 10.0 16 0.836
## 12 0.10 10.0 16 0.770
## 13 0.01 1.0 16 0.929
## 14 0.05 1.0 16 0.824
## 15 0.10 1.0 16 0.756
## 16 0.01 0.1 16 0.933
## 17 0.05 0.1 16 0.862
## 18 0.10 0.1 16 0.757
## 19 0.01 10.0 32 0.931
## 20 0.05 10.0 32 0.815
## 21 0.10 10.0 32 0.776
## 22 0.01 1.0 32 0.921
## 23 0.05 1.0 32 0.854
## 24 0.10 1.0 32 0.734
## 25 0.01 0.1 32 0.944
## 26 0.05 0.1 32 0.852
## 27 0.10 0.1 32 0.748
# Graphiques
par(mfrow = c(1, 3))
for (alpha in alpha_values) {
subset_data <- subset(results, Alpha == alpha)
plot(subset_data$A, subset_data$Coverage, type = "b", pch = 16,
xlab = expression(a), ylab = "Taux de couverture",
main = paste("Alpha =", alpha), ylim = c(0, 1),
col = as.factor(subset_data$N))
legend("bottomright", legend = paste("N =", n_values),
col = 1:length(n_values), lty = 1, pch = 16)
}
Question 15. Pendant une période de huit ans, nous
avons observé les hauteurs maximales en mètres suivantes pour le fleuve
:(x1,…,x8)=(2,27;3,75;0,29;6,88;5,17;4,29;4,07;4,47). Une compagnie
d’assurance estime qu’une crue catastrophique avec une hauteur de 6
mètres au moins n’arrive au plus qu’une fois tous les mille ans. Est-ce
justifié? Quel est le niveau de risque associé à cette décision ?
# Données fournies
altitudes <- c(2.27, 3.75, 0.29, 6.88, 5.17, 4.29, 4.07, 4.47)
# Statistiques descriptives
summary(altitudes)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.290 3.380 4.180 3.899 4.645 6.880
sd_altitudes <- sd(altitudes)
cat("Écart type :", sd_altitudes, "\n")
## Écart type : 1.95341
# Vérifier si les données suivent une distribution de Rayleigh
library(MASS)
# Ajustement du modèle de Rayleigh
fit <- fitdistr(altitudes, densfun = "weibull")
cat("Estimation des paramètres du modèle de Weibull :\n")
## Estimation des paramètres du modèle de Weibull :
print(fit)
## shape scale
## 2.0004250 4.3058384
## (0.6192116) (0.7848329)
# Graphique de la distribution observée
hist(altitudes, probability = TRUE, main = "Distribution des hauteurs maximales",
xlab = "Hauteur (m)", ylab = "Densité", col = "lightblue")
curve(dweibull(x, shape = fit$estimate[1], scale = fit$estimate[2]),
add = TRUE, col = "red", lwd = 2)
legend("topright", legend = c("Histogramme", "Ajustement Weibull"),
col = c("lightblue", "red"), lty = c(1, 1), lwd = c(2, 2))
# Test d'ajustement
ks_test <- ks.test(altitudes, "pweibull", fit$estimate[1], fit$estimate[2])
cat("Résultat du test de Kolmogorov-Smirnov :\n")
## Résultat du test de Kolmogorov-Smirnov :
print(ks_test)
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: altitudes
## D = 0.2816, p-value = 0.4677
## alternative hypothesis: two-sided