Université de Caen M1

SOLUTION TP no 7

Solution 1. On considère le jeu de données "rayons" disponible ici :

https://chesneau.users.lmno.cnrs.fr/rayons.txt

Dans celui-ci, on mesure les temps en secondes entre les arrivées de n = 3935 photons.
L’histogramme associé à ces valeurs est :

Vu son allure, il est raisonnable de penser que le temps en secondes entre deux arrivées de photons
est une var X suivant la loi gamma Γ(α, β), i.e. de densité :
 α
 β xα−1 e−βx si x ≥ 0,
f (x) = Γ(α)

0 sinon,
R∞
où Γ :]0, ∞[→ [0, ∞[ désigne la fonction gamma d’Euler définie par : Γ(α) = 0 tα−1 e−t dt.
Ici, α et β sont des paramètres inconnus que l’on souhaite estimer ponctuellement à l’aide des
données : x1 , . . . , xn (observation / réalisation d’un n-échantillon (X1 , . . . , Xn ) de X).

1. Montrer que
α α
E(X) = , V(X) = :
β β2
En faisant le changement de variable y = βx et l’égalité : Γ(α + 1) = αΓ(α), on obtient
Z ∞ Z ∞ Z ∞
β α α−1 −βx βα
E(X) = xf (x)dx = x x e dx = xα e−βx dx
−∞ Γ(α) Γ(α)
α Z ∞
 α 0 Z ∞
0
β y 1 1
= e−y dy = y (α+1)−1 e−y dy
Γ(α) 0 β β βΓ(α) 0
1 α
= Γ(α + 1) = .
βΓ(α) β

C. Chesneau 1 TP no 7 : solution
Université de Caen M1

On a V(X) = E(X 2 ) − (E(X))2 . Il reste à calculer E(X 2 ). En faisant de nouveau le

changement de variable y = βx et l’égalité : Γ(α + 2) = (α + 1)αΓ(α), on obtient
Z ∞ Z ∞ α Z ∞
2 β α−1 −βx βα
2
E(X ) = 2
x f (x)dx = x x e dx = xα+1 e−βx dx
−∞ 0 Γ(α) Γ(α) 0
α Z ∞
 α+1 Z ∞
β y 1 1
= e−y dy = 2 y (α+2)−1 e−y dy
Γ(α) 0 β β β Γ(α) 0
1 α(α + 1)
= 2
Γ(α + 2) = .
β Γ(α) β2
D’où
 2
α(α + 1) α α
V(X) = − = 2.
β2 β β

2. Mettre les données dans un vecteur x :

x = scan("https://chesneau.users.lmno.cnrs.fr/rayons.txt")

3. En estimant ponctuellement E(X) par x et V(X) par s2 , donner des estimations ponctuelles
de α et β par la méthode des moments : La question 1 entraîne α = βE(X), donc on peut
écrire :
α βE(X) E(X)
V(X) = 2
= 2
= ,
β β β
soit
α βE(X) E(X)
β= 2
= 2
= .
β β V(X)
Ainsi, une estimation de β par la méthode des moments est
x
b= .
s2
D’autre part, comme α = βE(X), une estimation ponctuelle de α est

x x2
a = bx = x = .
s2 s2
On considère les commandes :

b = mean(x) / sd(x)^2
a = mean(x)^2 / sd(x)^2
c(b, a)

Cela renvoie : [1] 0.01266144 1.01209495. Par la méthode des moments, une estimation
ponctuelle de α est 1.01209495 et une estimation ponctuelle de β est b = 0.01266144.

4. Consulter les aides suivantes de R : help(dgamma) et help(nlm).

C. Chesneau 2 TP no 7 : solution
Université de Caen M1

5. On cherche maintenant à estimer ponctuellement α et β par la méthode du maximum de

vraisemblance. Reproduire et comprendre les enjeux des commandes suivantes :

mlog = function(theta, x){

sum(-dgamma(x, shape = theta[1], rate = theta[2], log = TRUE))
}
nlm(mlog, c(1, 1), x = x)$estimate

Cela renvoie : [1] 1.02634117 0.01283969. En utilisant la méthode du maximum de

vraisemblance, on a ainsi une estimation ponctuelle de α qui est 1.02634117 et une estimation
ponctuelle de β qui est 0.01283969.

6. Quel est l’objet mathématique renvoyé par les commandes suivantes ?

mv = nlm(mlog, c(1, 1), x = x, hessian = TRUE)

mv$hessian

Cela renvoie :

[,1] [,2]
[1,] 6231.295 -305284.4
[2,] -305284.399 24121374.2

On obtient une estimation ponctuelle de la matrice Hessienne de la fonction de log-vraisemblance.

7. Quel sont les objets mathématiques renvoyés par les commandes suivantes ?

hat.alpha = mv$estimate[1]
hat.beta = mv$estimate[2]
z = qnorm((1 + 0.95) / 2)
inv.fish = solve(mv$hessian)
hat.alpha + c(-1, 1) * z * sqrt(inv.fish[1, 1])
hat.beta + c(-1, 1) * z * sqrt(inv.fish[2, 2])

Les commandes renvoient deux intervalles de confiance pour α et β au niveau 95% en utilisant
la convergence asymptotique de l’estimateur du maxium de vraisemblance.

Solution 2. Soient θ ∈ N∗ , X une var suivant la loi uniforme U({1, . . . , θ}) :

1
P(X = x) = , x ∈ {1, . . . , θ},
θ
n ∈ N∗ et (X1 , . . . , Xn ) un n-échantillon de X. Ici, θ est un paramètre inconnu que l’on souhaite
estimer à l’aide de (X1 , . . . , Xn ).

C. Chesneau 3 TP no 7 : solution
Université de Caen M1

1. Déterminer un estimateur obtenu par la méthode des moments, et l’estimateur du maximum

de vraisemblance : On a
θ
X X 1 1 θ(θ + 1) θ+1
E(X) = xP(X = x) = x× = × = .
x=1
θ θ 2 2
x∈X(Ω)

Il vient θ = 2E(X) − 1. Ainsi, la méthode des moments nous donne l’estimateur

n
2X
Un = Xi − 1.
n i=1

La fonction de vraisemblance est

n
Y 1
L(θ; x1 , . . . , xn ) = P(Xi = xi ) = , sup(x1 , . . . , xn ) ∈ {1, . . . , θ}.
i=1
θn

Comme cette fonction est décroissante en θ, on a

Vn = Argmax L(θ; X1 , . . . , Xn ) = sup(X1 , . . . , Xn )

θ∈N∗

Ainsi, la méthode du maximum de vraisemblance nous donne l’estimateur

Vn = sup(X1 , . . . , Xn ).

2. Reproduire et comprendre les enjeux des commandes suivantes :

simunif = function(teta, n, k) {
est1 = numeric(k)
est2 = numeric(k)
for (i in 1 : k) {
echant = sample(1 : teta, n, replace = T)
est1[i] = 2 * mean(echant) - 1
est2[i] = max(echant)
}
yrange = range(est1, est2)
plot(1:k, est1, type = "p", col = "blue", ylim = yrange, main = "", sub = "",
xlab = "", ylab = "")
points(1:k, est2, type = "p", col = "red", ylim = yrange, main = "", sub = "",
xlab = "", ylab = "")
abline(h = teta, col = "black")
abline(h = mean(est1), col = "blue")
abline(h = mean(est2), col = "red")
r1 = round(sum(est1 - teta)^2 / k, 2)
r2 = round(sum(est2 - teta)^2 / k, 2)
text(1, min(yrange) + 10, paste("r1 = ",r1," r2 = ",r2), pos = 4)
}
simunif(1000, 20, 200)
simunif(1000, 2000, 200)

C. Chesneau 4 TP no 7 : solution
Université de Caen M1

On constate que, quand n est petit, l’estimateur Un est meilleur, et quand n est grand, Vn
est meilleur.

C. Chesneau 5 TP no 7 : solution