Cours d'introduction à la statistique et à l'analyse de données : Notes de cours

Commençons par télécharger les données (avec la fonction wget) :

wget http://www.uow.edu.au/~mwand/webspr/janka.txt

janka.txt

Nous chargeons les données dans l'espace de travail :

janka <- read.table("janka.txt",header=TRUE)

Le jeu de données n'étant pas trop grand, nous pouvons l'afficher dans son ensemble :

janka

   dens hardness
1  24,7      484
2  24,8      427
3  27,3      413
4  28,4      517
5  28,4      549
6  29,0      648
7  30,3      587
8  32,7      704
9  35,6      979
10 38,5      914
11 38,8     1070
12 39,3     1020
13 39,4     1210
14 39,9      989
15 40,3     1160
16 40,6     1010
17 40,7     1100
18 40,7     1130
19 42,9     1270
20 45,8     1180
21 46,9     1400
22 48,2     1760
23 51,5     1710
24 51,5     2010
25 53,4     1880
26 56,0     1980
27 56,5     1820
28 57,3     2020
29 57,6     1980
30 59,2     2310
31 59,8     1940
32 66,0     3260
33 67,4     2700
34 68,8     2890
35 69,1     2740
36 69,1     3140

Comme toujours, la première chose à faire est une figure de la variable à expliquer en fonction de la (des) variable(s) explicatives (ou prédictives) :

par(cex=2)
with(janka,plot(dens,hardness,col="blue",xlab="Densité",ylab="Dureté"))
grid(10,10,lwd=2)

On peut remarquer plusieurs choses :

• la densité apparaît comme un bon prédicteur de la dureté ;
• la variance de la dureté semble augmenter avec la densité ;
• la « fonction de régression » semble être légèrement incurvée vers le haut pour les plus grandes valeurs de la densité.

J'utilise ici GNU octave un logiciel libre qui ressemble beaucoup à Matlab.

Je commence par charger les données :

janka = dlmread("janka.txt")
janka = janka(2:37,1:2)

Je construis ensuite la matrice du plan d'expérience :

X = [ones(36,1) janka(:,1) janka(:,1).^2]

J'obtiens la décomposition QR de celle-ci :

[Q,R] = qr(X)

Puis \(\hat{β}\) :

beta_chapeau = R(1:3,1:3)\(Q(:,1:3)'*janka(:,2))
ans = beta_chapeau

GNU octave ne semble pas avoir l'équivalent de backsolve.

Je répète l'opération en travaillant avec la log(dureté) et en centrant la densité :

Xb = [ones(36,1) janka(:,1)-mean(janka(:,1)) (janka(:,1)-mean(janka(:,1))).^2]
[Qb,Rb] = qr(Xb)
beta_chapeau_b = Rb(1:3,1:3)\(Qb(:,1:3)'*log(janka(:,2)))
ans = beta_chapeau_b

Il suffit alors de continuer de la sorte pour obtenir des résultats équivalents, pour ne pas dire identiques, à ceux de R…

Les plus aventureux pourons jetter un oeil sur une manière de reproduire cette analyse en utilisant le Common Lisp.

Nous allons construire des intervalles de confiances suivant deux méthodes : celle de Wald et celle des rapports de vraisemblance, en réduisant encore la taille de l'échantillon considéré, c.-à-d. en prenant une taille de 5 au lieu de la plus petite taille de 10. Commençons donc par construire un échantillon de taille 5 :

echantillon5 <- sapply(1:nombreDeRepliques,
                       function(i) rexp(5,1/tau.vrai))

D'après la théorie ci-dessus, si la limite asymptotique est atteinte, l'EMV de τ, c'est-à-dire dans ce cas la moyenne empirique, doit suivre une loi normale centrée sur la vraie valeur (1) avec une variance égale à \(1^2/n\), c'est-à-dire, 0.2. Nous pouvons ici construire une estimation de la densité de probabilité de \(\hat{τ}\) avec un histogramme et lui superposer la densité théorique :

par(cex=2)
hist(apply(echantillon5,2,mean),breaks=50,
     main="",xlab=expression(hat(tau)),
     ylab="Densité empirique",probability=TRUE)
toto <- seq(-0.5,2.5,len=501)
lines(toto,dnorm(toto,1,sqrt(0.2)),col=2,lwd=2)

On voit ici deux problèmes :

• la « théorie » prise au pieds de la lettre nous prévoit des valeurs possibles de \(\hat{τ}\) que nous savons être impossibles car négatives (puisque l'EMV est par définition dans ce cas une moyenne empirique d'observations strictement positives) ;
• la densité empirique (l'histogramme) est bien plus asymétrique que ne le prévoit la loi asymptotique, le mode (pic) est un peu à gauche de 1 et d'assez grandes déviations positives sont observées.

Une façon plus quantitative de vérifier l'adéquation entre un histogramme – construit avec une taille de bin suffisamment petite – est une densité théorique est le hanging rootogram. L'idée est la suivante : soit \(p_i\) la probabilité d'observer un événement dans le bin i – théoriquement on a \(p_i = \int_{g_i}^{d_i} f(x) \, dx\), où f est la densité théorique et où le bin i est défini comme l'intervalle \([g_i,d_i]\) et soit n le nombre total d'observations IIDs effectuées. Le nombre d'observations dans le bin i devrait suivre, si la théorie est exacte, une loi binomiale de paramètres \(p_i\) et n (c'est-à-dire avoir une espérance de \(n \, p_i\) et une variance \(n \, p_i \, (1-p_i)\)). Maintenant si \(n\) devient très grand et si \(p_i\) devient très petit avec un produit \(λ_i \equiv n \, p_i\) constant, la loi binomiale tend vers une loi de Poisson de paramètre \(λ_i\). Si nos hypothèses sur f (ici une loi normale centré sur 1 de variance 0.2) sont correctes nos comptes par bin devraient suivre approximativement une loi de Poisson. Cela implique que si le bin est dans une région où f est grande l'amplitude des écarts entre prédiction théorique et observations devraient être typiquement plus grands que si le bin est dans une région où f est petite (puisque que la variance d'une variable aléatoire de Poisson est égale à sa moyenne) cela signifie que nous ne devrions pas avoir des écarts homogènes entre densité empirique (histogramme) et densité théorique; dit autrement, juger à l'œil l'écart entre la courbe rouge et l'histogramme de la figure précédente pourrait nous induire en erreur. Pour rendre une telle comparaison plus objective, nous allons chercher une transformation de nos observations qui stabilise leur variance. Pour ce faire nous partons de la convergence en loi d'une loi de Poisson paramètre \(λ_i\) vers une loi normale d'espérance et de variance \(λ_i\) – ça se monte bien avec les fonctions caractéristiques, les transformées de Fourier, pour vous en convaincre, faites quelques figures avec \(λ_i \ge 25\). Si nous notons \(B_i\) la variable aléatoire modélisant notre nombre d'observations dans le bin i nous avons approximativement : \[ B_i \sim \mathcal{N}(λ_i,λ_i) \] et nous pouvons écrire pour chacune de nos observations : \[ b_i = λ_i + \sqrt{λ_i} \, ε_i \] où \(ε_i\) est la réalisation d'une variable aléatoire normale centrée réduite (de moyenne 0 et de variance 1). Si nous définissons une nouvelle variable aléatoire \(Z_i \equiv g(B_i)\), où g est une fonction inversible (on peut démontrer ce qui suit avec des conditions moins restrictives sur g) on a alors : \[ z_i \approx g(λ_i) + \sqrt{λ_i} \, g'(λ_i) \, ε_i \; .\] Dit autrement, \(Z_i\) est approximativement normale de moyenne \(g(λ_i)\) et de variance \(λ_i \, g'(λ_i)^2\) (en Physique vous appelez ça la propagation des incertitudes ou des erreurs, en Stats cela s'appelle la méthode delta). Ici il nous reste à trouver g telle que \(\sqrt{x} \, g'(x) = Cste\) soit \(g'(x) α 1/\sqrt{x}\). On voit donc que \(g(x) \equiv \sqrt{x}\) fait ce qu'on veut. Ainsi, avec \(Z_i = \sqrt{B_i}\) on a : \[ z_i \approx \sqrt{λ_i} + 1/2 \, ε_i \] et \(Z_i\) est approximativement normale de moyenne \(\sqrt{λ_i}\) et de variance 1/4. On pourra s'en convaincre facilement avec une petite simulation :

λV <- c(25,50,100,500,1000)
rbind(λV,
      sapply(λV,
             function(λ) var(sqrt(rpois(10000,λ)))))

Pour faire notre test visuel nous prenons la racine des comptes par bin et nous soustrayons la racine carrée de la densité théorique multipliée par la taille du bin et par le nombre total d'événements dans l'échantillon. Nous pouvons rajouter au graphe des lignes horizontales correspondant à plus ou moins 2 fois l'écart type théorique (0.5) :

par(cex=2)
mesBornes <- seq(0,3.5,0.1)
histTaille5 <- hist(apply(echantillon5,2,mean),
                    breaks=mesBornes,plot=FALSE)
theoTaille5 <- diff(pnorm(histTaille5$breaks,1,sqrt(0.2)))*dim(echantillon5)[2]
plot(histTaille5$mids,
     sqrt(histTaille5$counts)-sqrt(theoTaille5),
     type="l",lwd=2,
     main="",xlab=expression(hat(tau)),
     ylab="")
abline(h=c(-1,1),col=2,lty=2)
abline(h=0,col="grey80")

On peut comparer au même graphe construit avec la plus grande taille d'échantillon simulé (1000) :

par(cex=2)
σ <- sqrt(1/1000)
mesBornes <- seq(-5,5,0.2)*σ+1
histTaille1000 <- hist(dureeDeVie[,100],
                       breaks=mesBornes,plot=FALSE)
theoTaille1000 <- diff(pnorm(histTaille1000$breaks,1,σ))*dim(dureeDeVie)[1]
plot(histTaille1000$mids,
     sqrt(histTaille1000$counts)-sqrt(theoTaille1000),
     type="l",lwd=2,
     main="",xlab=expression(hat(tau)),
     ylab="",ylim=range(sqrt(histTaille5$counts)-sqrt(theoTaille5)))
abline(h=c(-1,1),col=2,lty=2)
abline(h=0,col="grey80")

Nous venons de construire des hanging rootograms pour la statistique de Wald comparée à sa loi asymptotique avec deux tailles d'échantillons : 5 et 1000. On a vu clairement qu'avec la plus petite taille d'échantillon, la limite asymptotique n'a pas été atteinte. Nous pouvons répéter l'opération avec une autre statistique, celle du rapport de vraisemblance. Pour la cas qui nous occupe, comme l'EMV, \(\hat{τ}\) est la moyenne empirique, on peut écrire la log-vraisemblance : \[l(τ) = -n \, \left(\log τ + \hat{τ}/τ\right)\] et la statistique du rapport de vraisemblance prend la forme – en notant τ₀ la vraie valeur de τ : \[\hat{r} = - 2 \, n \, \left( \log (\hat{τ}/τ_0) + 1 - \hat{τ}/τ_0 \right) \; .\] Dans le cadre de nos simulation nous avons de plus : τ₀ = 1. Il ne nous reste plus qu'à calculer la valeur de la statistique pour chaque simulation, construire une estimation de la densité (c.-à-d. un histogramme) à partir de ses valeurs, calculer la valeur de \(n \, \int_{g_i}^{d_i} χ^2_1(x) \, dx\) – ici la loi asymptotique est un χ² à un degré de liberté – puis prendre les racines carrées avant de former la différence ¹.

par(cex=2)
τ.chapeau <- apply(echantillon5,2,mean)
rV <- -5*(log(τ.chapeau/tau.vrai) + 1- τ.chapeau/tau.vrai)*2
yl <- range(rV)
mesBornes <- seq(yl[1]-0.5,yl[2]+0.5,len=51)
histTaille5 <- hist(rV,
                    breaks=mesBornes,plot=FALSE)
chi2Taille5 <- diff(pchisq(histTaille5$breaks,1))*dim(echantillon5)[2]
plot(histTaille5$mids,
     sqrt(histTaille5$counts)-sqrt(chi2Taille5),
     type="l",lwd=2,
     main="",xlab=expression(hat(r)),
     ylab="")
abline(h=c(-1,1),col=2,lty=2)
abline(h=0,col="grey80")

On constate que même pour une aussi petite taille d'échantillon, la statistique semble assez proche de sa valeur asymptotique. On est donc pas surpris par le comportement observé avec la plus grande taille d'échantillon (1000) :

par(cex=2)
τ.chapeau <- dureeDeVie[,100]
rV <- -1000*(log(τ.chapeau/tau.vrai) + 1- τ.chapeau/tau.vrai)*2
yl <- range(rV)
mesBornes <- seq(yl[1]-0.1,yl[2]+0.1,len=51)
histTaille1000 <- hist(rV,
                       breaks=mesBornes,plot=FALSE)
chi2Taille1000 <- diff(pchisq(histTaille1000$breaks,1))*dim(dureeDeVie)[1]
plot(histTaille1000$mids,
     sqrt(histTaille1000$counts)-sqrt(chi2Taille1000),
     type="l",lwd=2,
     main="",xlab=expression(hat(r)),
     ylab="",ylim=range(sqrt(histTaille5$counts)-sqrt(chi2Taille5)))
abline(h=c(-1,1),col=2,lty=2)
abline(h=0,col="grey80")

On construit immédiatement des intervalles de confiances avec la méthode de Wald : \[ [\hat{τ} + φ^{-1}\left(α/2\right) \, \sqrt{\hat{τ}/n} , \hat{τ} - φ^{-1}\left(α/2\right) \, \sqrt{\hat{τ}/n} ] \; ,\] où \(φ\) désigne la fonction de répartition d'une loi normale centrée réduite et \(φ^{-1}\) est son inverse ; α est le niveau du test. On peut donc, pour chacune de nos 1000 expériences simulées, tracer un intervalle de confiance à 95% est regarder si la vraie valeur (1) s'y trouve bien 95 fois sur 100, c'est-à-dire à peu près 950 fois sur 1000 dans notre cas. Clairement la validité des intervalles de confiance dépend de la validité de l'approximation asymptotique de la statistique utilisée. Compte tenu de ce que nous avons vu dans les deux sections précédentes, nous ne nous attendons pas à ce que les intervalles de confiance à 95% basés sur la statistique de Wald aient une probabilité de recouvrement empirique identique à leur probabilité nominale (ou théorique).

par(cex=1.5)
τ.chapeau <- apply(echantillon5,2,mean)
τ.IC <- τ.chapeau + qnorm(0.975) * (c(-1,1) %o% sqrt(τ.chapeau/5))
τ.vrai.dedans <- apply(τ.IC,2, function(x) x[1] <= 1 && 1 <= x[2])
plot(c(0,1000),
     range(τ.IC),
     type="n",
     xlab="Indice de réplique",
     ylab="",
     main="Intervalles basés sur la méthode de Wald")
abline(h=1,col="orange",lwd=2)
invisible(sapply(1:1000,
                 function(i) {
                   x <- τ.IC[,i]
                   dedans <- τ.vrai.dedans[i]
                   segments(i,x[1],i,x[2],col=ifelse(dedans,1,2))}))

On a ici seulement :

sum(τ.vrai.dedans)

914

qui contiennent la vraie valeur.

Pour la méthode du rapport de vraisemblance, nous avons, lorsque la limite asymptotique est atteinte : \[2 \, \left( l(\hat{τ}) - l(τ_0) \right) \sim χ^2_1 \; .\] Ainsi : \[ \mathrm{Pr} \{ 2 \, \left( l(\hat{τ}) - l(τ_0) \right) \le z \} = \mathrm{pchisq}(z,1) \; ,\] où \(\mathrm{pchisq}(z,1)\) désigne la valeur de la fonction de répartition d'un loi \(χ^2_1\) en z. Nous allons donc construire un intervalle de niveau α en cherchant τ_g ≤ \(\hat{τ}\) ≤ τ_d, tels que : \[ 2 \, \left( l(\hat{τ}) - l(τ_{g,d}) \right) = \mathrm{qchisq}(1-α,1) \; \] où \(\mathrm{qchisq}(1-α,1)\) désigne la valeur de l'inverse de la fonction de répartition d'un χ²₁ en 1-α. Pour trouver les valeurs \(τ_{g,d}\) nous devons procéder numériquement :

q95 <- qchisq(0.95,1)
τ.IC <- sapply(1:1000,
               function(i) {
                 τ.emv <- τ.chapeau[i]
                 lFct <- function(τ) -5*(log(τ) + τ.emv/τ)
                 lAτ.chapeau <- lFct(τ.emv)
                 cible <- function(τ) 2*(lAτ.chapeau - lFct(τ)) - q95 
                 icBas <- uniroot(cible,c(0.001,τ.emv))$root
                 icHaut <- uniroot(cible,c(τ.emv,100))$root
                 c(icBas,icHaut)
               })
τ.vrai.dedans <- apply(τ.IC,2, function(x) x[1] <= 1 && 1 <= x[2])
plot(c(0,1000),
     range(τ.IC),
     type="n",
     xlab="Indice de réplique",
     ylab="",
     main="Intervalles basés sur le rapport de vraisemblance")
abline(h=1,col="orange",lwd=2)
invisible(sapply(1:1000,
                 function(i) {
                   x <- τ.IC[,i]
                   dedans <- τ.vrai.dedans[i]
                   segments(i,x[1],i,x[2],col=ifelse(dedans,1,2))}))

On a cette fois :

sum(τ.vrai.dedans)

946

intervalles contenant la bonne valeur…

Nous considérons les données de fertilité de la communauté Huttérite étudiées par R. Ihaka dans son 4e cours. En nous basant sur l'analyse d'Ihaka, nous « adoptons comme vraie » la fonction de répartition du temps (en mois) entre marriage et première grossesse suivante : \[F(t) = 1 - \frac{h^a}{(h+t)^a} \quad \textrm{avec} \quad a \, = \, 3,1 \quad \textrm{et} \quad h \, = \, 9,6 \; . \]

L'échantillon considéré par Ihaka comprend 342 femmes. Les temps entre mariage et grossesse avaient de plus été tabulé suivant une partition dont les limites étaient :

(breaks <- c(0:12,15,18,24,Inf))

[1]   0   1   2   3   4   5   6   7   8   9  10  11  12  15  18  24 Inf

Pour mémoire, les données sont :

Nous allons supposer que nous disposons de données non discrétisées et nous souhaitons comparer les performances de trois estimateurs du couple (a,h) de paramètres du modèles :

• l'estimateur du maximum de vraisemblance appliqué aux données brutes ;
• l'estimateur du maximum de vraisemblance appliqué aux données tabulées selon la même partition que les données originale ;
• l'estimateur du moindre χ², c'est-à-dire celui utilisé par Ihaka.

Nous allons faire cette étude par simulation, c'est-à-dire que nous allons simuler, disons 1000 fois, une échantillon de 342 variables IID dont la fonction de répartition est donnée par notre première équation. Pour chacun de nos 1000 échantillons, nous obtiendrons 3 estimations du couple (a,b) de paramètres et nous comparerons nos estimations en construisant leurs histogrammes.

La Méthode de la transformée inverse, toujours applicable mais souvent coûteuse consiste en générer un nombre, U, suivant une loi uniforme sur [0,1] avant de trouver X tel que F(X) = U, où F est la fonction de répartition de la loi cible. X a alors F pour fonction de répartition. On prouve cela facilement lorsque F est strictement croissante, on a en effet : \[X \sim \mathrm{F}^{-1}(U)\, , \quad \mathrm{où} \quad U \sim \mathcal{U}(0,1) \; .\] Par définition: \[\mathrm{F}_X(x) \equiv \mathrm{Pr}\{X \le x\} \; ,\] mais \[\begin{array}{l l l} \mathrm{Pr}\{X \le x\} & = & \mathrm{Pr}\{\mathrm{F}(X) \le \mathrm{F}(x)\} \quad \mathrm{(F} \; \textrm{strictement croissante)} \\ \mathrm{Pr}\{X \le x\} & = & \mathrm{Pr}\{U \le \mathrm{F}(x)\} \quad \textrm{(définition de} \; X \mathrm{)} \\ \mathrm{Pr}\{X \le x\} & = & \mathrm{F}(x) \quad \textrm{(car} \; U \; \textrm{est uniforme sur} \; [0,1) \textrm{)} \\ \mathrm{F}_X(x) & = & \mathrm{F}(x) \; . \end{array}\]

La mise en oeuvre « brutale » de cette méthode est numérique et utilise la recherche de la racine de la fonction U-F(X), avec uniroot dans R. On pourrait ainsi générer n observations suivant notre loi avec la fonction :

rHutterite <- function(n)
  sapply(runif(n), function(u) uniroot(function(x) u-1+1/(1+x/9.6)^3.1,c(0,100))$root)

Une façon plus élégante et plus efficace de faire la même chose et d'inverser analytiquement notre fonction de répartition. Nous avons en effet : \[\begin{array}{l l l} F(t) & = & 1 - \frac{h^a}{(h+t)^a} \\ F(t) & = & 1 - \frac{1}{(1+t/h)^a} \\ 1 - F(t) & = & \frac{1}{(1+t/h)^a} \\ \frac{1}{1 - F(t)} & = & (1+t/h)^a \\ \frac{1}{\sqrt[a]{1 - F(t)}} & = & (1+t/h) \\ \frac{h}{\sqrt[a]{1 - F(t)}} - h & = & t \end{array}\]

Nous pouvons ainsi définir une suite de fonctions liées à la loi qui nous intéresse en suivant les conventions (implicites) de R :

dHutterite <- function(x, a=3.1, h=9.6) a/h/(1+x/h)^(a+1)
pHutterite  <- function(q, a=3.1, h=9.6) 1-1/(1+q/h)^a
qHutterite <- function(p, a=3.1, h=9.6) h*((1-p)^(-1/a)-1)
rHutterite <- function(n, a=3.1, h=9.6) qHutterite(runif(n),a,h)

On vérifie facilement que de grosses erreurs n'ont pas été commises dans nos définitions de la fonction de répartition (pHutterite) et de la fonction quantile (qHutterite) avec :

pp <- seq(0,1,len=101)
max(abs(pHutterite(qHutterite(pp))-pp))

2.98372437868011e-16

Nous avons désormais tout ce qu'il nous faut pour répondre à notre question. Un truc à utiliser pour gagner du temps et de commencer nos optimisations avec la vraie valeur des paramètre du modèle :

nrep <- 1000
set.seed(20061001,"Mersenne-Twister")
N <- 342
sim <- sapply(1:nrep,
              function(i) {
                obs <- rHutterite(N)
                fit1 <- optim(log(c(3.1,9.6)),function(p) moinsLogVraiBrute(exp(p),obs),method="BFGS")
                fit2 <- optim(log(c(3.1,9.6)),function(p) moinsLogVraiTab(exp(p),obs),method="BFGS")
                fit3 <- optim(log(c(3.1,9.6)),function(p) chi2Tab(exp(p),obs),method="BFGS")
                c(fit1conv = fit1$convergence,
                  fit1a = fit1$par[1],
                  fit1h = fit1$par[2],
                  fit2conv = fit2$convergence,
                  fit2a = fit2$par[1],
                  fit2h = fit2$par[2],
                  fit3conv = fit3$convergence,
                  fit3a = fit3$par[1],
                  fit3h = fit3$par[2]) } )
c(sum(sim[1,] > 0),sum(sim[4,] > 0),sum(sim[7,] > 0))

Soit N le nombre total d'enfants de la famille (de trois enfants) infectés avant la fin de l'épidémie. Montrez que : \[ \mathrm{Pr}\left(N=1\right) = (1-\theta)^2 \; ; \quad \mathrm{Pr}\left(N=2\right) = 2\, \theta \, (1-\theta)^2 \; ; \quad \mathrm{Pr}\left(N=3\right) = \theta^2 \, (3-2\,\theta) \; ; \] Vous supposerez que les infections ont lieu indépendamment et vous pourrez vous rappeler que la somme des probabilités de tous les cas possibles doit être égale à 1.

• La probabilité pour qu'aucun des deux autres enfants ne soit infecté est, du fait de l'indépendance, le produit des probabilités pour qu'ils ne le soient pas individuellement d'où : \(\mathrm{Pr}\left(N=1\right) = (1-\theta)^2\) ;
• Pour qu'un seul enfant soit infecté il faut qu'il soit infecté (probabilité θ) et que l'autre ne le soit ni par le premier enfant (probabilité 1-θ) ni par celui nouvellement infecté (probabilité 1-θ), d'où la probabilité pour que un enfant donné soit infecté alors que l'autre ne l'est pas: \(\theta \, (1-\theta)^2\) et comme il y a deux enfants, la probabilité pour qu'on est au total 2 enfants infectés : \(2 \, \theta \, (1-\theta)^2\) ;
• on doit avoir : \(\mathrm{Pr}\left(N=1\right) + \mathrm{Pr}\left(N=2\right) + \mathrm{Pr}\left(N=3\right) = 1\), soit \(\mathrm{Pr}\left(N=3\right) = 1 - \mathrm{Pr}\left(N=1\right) - \mathrm{Pr}\left(N=2\right)\), \(\mathrm{Pr}\left(N=3\right) = 1 - (1-\theta)^2 - 2 \, \theta \, (1-\theta)^2\). En développant le deuxième terme du second membre on obtient : \(\mathrm{Pr}\left(N=3\right) = 2\, \theta-\theta^2 - 2 \, \theta \, (1-\theta)^2\), en développant tout il vient : \(\mathrm{Pr}\left(N=3\right) = 2\, \theta-\theta^2 - 2 \, \theta + 4 \, \theta^2 - 2 \, \theta^3\) et en simplifiant : \(\mathrm{Pr}\left(N=3\right) = \theta^2 \, (3-2\,\theta)\).

Vous écrirez les expressions de la log-vraisemblance, de la fonction score et de la fonction d'information observée lorsqu'un échantillon de K familles de trois enfants est observé et que k₁ familles restent avec un seul enfant infecté, k₂ familles en ont deux et k₃ familles en ont trois. Les familles sont bien sûr supposées indépendantes (au sens statistique du terme). Inutile de vous « alourdir » avec des termes ne faisant pas intervenir le paramètre θ.

En supposant l'indépendance des familles on a : \[\mathrm{Pr}\left(K_1=k_1,K_2=k_2,K_3=k_3 | \theta \right) = C \, \mathrm{Pr}\left(N=1\right)^{k_1} \, \mathrm{Pr}\left(N=2\right)^{k_2} \, \mathrm{Pr}\left(N=3\right)^{k_3} \; ,\] où C est une constante de normalisation ne dépendant que de K, k₁, k₂ et k₃ – on a ici affaire à une loi multinomiale et \(C=K!/(k_1 ! \, k_2 ! \, k_3 !)\).

La log vraisemblance devient : \[l(\theta) = 2 \, (k_1 + k_2) \, \log (1-\theta) + (k_2 + 2 \, k_3) \, \log \theta + k_3 \, \log (3-2 \, \theta ) \; .\]

La fonction score : \[ S(\theta) = - \frac{2 \, (k_1 + k_2)}{1 - \theta} + \frac{k_2+2\, k_3}{\theta} - \frac{2 \, k_3}{3 - 2\, \theta} \; .\]

L'information observée : \[\mathcal{J}(\theta) = \frac{2 \, (k_1 + k_2)}{(1 - \theta)^2} + \frac{k_2+2\, k_3}{\theta^2} + \frac{4 \, k_3}{(3 - 2\, \theta)^2} \; .\]

On voit que comme 0 ≤ k₁, k₂ et k₃ ≤ K et k₁ + k₂ + k₃ = K, l'information observée est toujours positive, ce qui implique que \(l(\theta)\) est concave ce qui veut dire que nous n'aurons pas de problème de convergence des méthodes d'optimisation numérique.

Par curiosité, je donne les commandes permettant d'effectuer ces petits calculs avec le logiciel Maxima (les commandes seraient très proches avec Mapple). Définition de la log vraisemblance:

2*(k1+k2)*log(1-θ) + (k2+2*k3)*log(θ)+k3*log(3-2*θ)

print(<<log-vraisemblance()>>);

(2 k3 + k2) log(θ) + 2 (k2 + k1) log(1 - θ) + k3 log(3 - 2 θ) 

La fonction score est alors simplement obtenue :

diff(2*(k1+k2)*log(1-θ)+(k2+2*k3)*log(θ)+k3*log(3-2*θ),θ,1);

diff(<<log-vraisemblance()>>,θ,1)

display2d:false;print(<<score>>);

(2*k3+k2)/θ-2*(k2+k1)/(1-θ)-2*k3/(3-2*θ) 

De même que la fonction d'information observée:

-diff(2*(k1+k2)*log(1-θ)+(k2+2*k3)*log(θ)+k3*log(3-2*θ),θ,2);

-diff(<<log-vraisemblance()>>,θ,2)

display2d:false;print(<<information-observée>>);

(2*k3+k2)/θ^2+2*(k2+k1)/(1-θ)^2+4*k3/(3-2*θ)^2 

Une étude portant sur K = 100 familles de trois enfants donne les résultats suivants :

Vous construirez une fonction R renvoyant la log-vraisemblance pour ce jeu de données.

Je procède comme à mon habitude en deux temps en définissant d'abord une fonction qui construit la fonction de log-vraisemblance² :

faitFLV <- function(k1=48,k2=32,k3=20)
  function(θ) 2*(k1+k2)*log(1-θ) + (k2+2*k3)*log(θ) + k3 * log(3-2*θ)

Puis je construis la fonction de log vraisemblance :

FLVvraie <- faitFLV(k1=48,k2=32,k3=20)

Vous pouvez aussi définir FLVvraie de la façon suivante:

FLVvraie <- function(θ,k1=48,k2=32,k3=20)
  2*(k1+k2)*log(1-θ) + (k2+2*k3)*log(θ) + k3 * log(3-2*θ)

Comme il est plutôt difficile de trouver la bonne racine de la fonction score analytiquement (avec maxima c'est faisable), nous allons le faire numériquement avec optimise. Comme cette dernière fonction demande à son utilisateur de fournir des bornes encadrant l'optimum, nous allons commencer par tracer un graphe de la log vraisemblance sur le domaine de définition de θ : (0,1) pour trouver une limite raisonnable au domaine fourni à optimise. Vous prendrez soin de mettre des labels à vos axes. Vous pourrez néanmoins regarder de près votre information observée et commenter sur la nécessité de définir un domaine « raisonnable » plus court que (0,1).

Le code suivant nous donne directement la figure (\ref{fig:Stat2013reponse4}) demandée. (J'ai ajusté la taille des labels, par(cex=2) pour que la figure soit plus lisible ; je ne vous demandais pas de le faire dans vos fichiers de script. J'ai aussi utilisé des commandes un peu sophistiquées pour avoir des labels avec des lettres grecs, ce que je ne vous demandais pas de faire non plus.)

θθ <- seq(0,1,len=501)
par(cex=2)
plot(θθ, FLVvraie(θθ),
     type="l",lwd=3,
     xlab=expression(theta),
     ylab=expression(l(theta)))

Comme mentionné dans l'énonce de la question, maxima nous permet d'obtenir analytiquement les racines de la fonction score:

solve(diff(2*(k1+k2)*log(1-θ)+(k2+2*k3)*log(θ)+k3*log(3-2*θ),θ,1),θ);

solve(diff(<<log-vraisemblance()>>,θ,1),θ)

display2d:false;print(<<racines-de-score>>);

[θ = -(sqrt((48*k2+48*k1)*k3+49*k2^2+84*k1*k2+36*k1^2)-12*k3-11*k2-6*k1)
   /(12*k3+12*k2+8*k1),
 θ = (sqrt((48*k2+48*k1)*k3+49*k2^2+84*k1*k2+36*k1^2)+12*k3+11*k2+6*k1)
   /(12*k3+12*k2+8*k1)]
  

Nous pouvons directement utiliser ces expression pour construire une fonction R

racineScore <- function(k1=48,k2=32,k3=20)
  c(-(sqrt((48*k2+48*k1)*k3+
           49*k2^2+
           84*k1*k2+
           36*k1^2)-
      12*k3-
      11*k2-
      6*k1)/(12*k3+12*k2+8*k1),
    (sqrt((48*k2+48*k1)*k3+
          49*k2^2+
          84*k1*k2+
          36*k1^2)+
     12*k3+
     11*k2+
     6*k1)/(12*k3+12*k2+8*k1))

Vous trouverez la valeur de votre estimateur du maximum de vraisemblance numériquement (avec, par exemple, optimise).

(θ.chapeau <- optimize(FLVvraie,lower=0.2,
                       upper=0.4,maximum=TRUE)$maximum)

[1] 0,2954604

Les expressions analytiques des racines de la fonction score nous donnent :

racineScore(k1=48,k2=32,k3=20)

[1] 0,2954474 1,4505843

Vous préparerez un nouveau graphe avec la log vraisemblance en trait continu noir, le point maximum figuré par un disque rouge (par exemple) et son abscisse figurée par une ligne verticale pointillée. Vous êtes invités à faire un choix « judicieux » d'échelle pour vos abscisses et vos ordonnées.

Là encore, la courte séquence de commandes suivante génère la figure demandée (\ref{fig:Stat2013reponse6}).

par(cex=2)
plot(θθ, FLVvraie(θθ),
     type="l",lwd=3,
     xlab=expression(theta),
     ylab=expression(l(theta)),
     xlim=c(0.2,0.4),
     ylim=c(FLVvraie(0.2),FLVvraie(θ.chapeau)))
abline(v=θ.chapeau,col="red",lty=3)
points(θ.chapeau,FLVvraie(θ.chapeau),col="red",pch=16)

Vous utiliserez votre expression de l'information observée que vous combinerez avec votre estimateur \(\hat{\theta}\) de θ pour obtenir un intervalle de confiance à 95% avec la méthode de Wald.

Je commence par définir une fonction qui renvoie une fonction qui renvoie l'information observée

faitInfObsF <- function(k1=48,k2=32,k3=20)
  function(θ) 2*(k1+k2)/(1-θ)^2 + (k2+2*k3)/θ^2 + 4*k3/(3-2*θ)^2

Puis j'appelle cette fonction pour obtenir ma fonction d'information observée :

InfObsFvraie <- faitInfObsF(k1=48,k2=32,k3=20)

La valeur de l'information observée au maximum de vraisemblance est donc :

(J.en.θ.chapeau <- InfObsFvraie(θ.chapeau))

[1] 1160,893

D'où l'erreur type sur \(\hat{\theta}\) (la racine carrée de l'inverse de l'information observée au maximum de vraisemblance) :

(σ.θ.chapeau <- 1/sqrt(J.en.θ.chapeau))

[1] 0,02934972

L'intervalle de confiance à 95% avec la méthode de Wald est donc :

θ.chapeau + c(-1,1)*qnorm(0.975)*σ.θ.chapeau

[1] 0,2379360 0,3529848

Vous simulerez 1000 ou 10000 répliques de l'expérience ci-dessus en utilisant pour θ, dans les formules donnant les probabilités d'observer un, deux ou trois enfants malades, le \(\hat{\theta}\) que vous venez d'obtenir. Vous utiliserez la fonction rmultinom de R pour vos simulations; vous commencerez par lire la documentation associée et vous verrez que dans votre cas, l'argument prob doit être un vecteur de trois éléments : les probabilités d'avoir 1, 2 ou 3 enfants malades. Vous n'oublierez pas de spécifier la graine de votre générateur de nombres (pseudo)aléatoires.

Je commence par définir une fonction de θ qui renvoie le vecteur de probabilité :

faitVecProb <- function(θ)
  c(un=(1-θ)^2,
    deux=2*θ*(1-θ)^2,
    trois=θ^2*(3-2*θ))

Je vérifie que les éléments du vecteur se somment à 1 :

sum(faitVecProb(θ.chapeau))

1

Je fais la simulation :

set.seed(20061001)
n.rep <- 10000
nb.familles <- 100
sim.exp <- rmultinom(n.rep,
                     nb.familles,
                     prob=faitVecProb(θ.chapeau))

Pour chacune des expériences simulées vous obtiendrez l'estimateur du maximum de vraisemblance comme vous l'avez fait avec les vraies données, vous construirez alors un histogramme (avec entre 50 et 100 bins) des \(\hat{\theta}\) auquel vous superposerez la loi asymptotique (loi normale centrée sur la vraie valeur dont la variance est l'inverse de l'information de Fisher ; pour avoir l'info de Fisher vous remplacez dans l'info observée les k₁, k₂, k₃ par leur valeur théorique).

Je commence par obtenir l'estimateur du maximum de vraisemblance pour chaque réplique :

θ.chapeau.sim <- apply(sim.exp,2,
                       function(K) {
                         FLV <- faitFLV(k1=K[1],k2=K[2],k3=K[3])
                         optimize(FLV,lower=0,
                                  upper=1,
                                  maximum=TRUE)$maximum
                       })

J'obtiens l'information de Fisher :

Kideal <- nb.familles * faitVecProb(θ.chapeau)
InfFisher <- faitInfObsF(k1=Kideal[1],k2=Kideal[2],k3=Kideal[3])

Je construis l'histogramme et je superpose la densité asymptotique comme explicité dans le code ci-dessous qui donne la figure \ref{fig:Stat2013reponse9}.

par(cex=2)
hist(θ.chapeau.sim,
     breaks=100,
     prob=TRUE,
     xlab=expression(hat(theta)),
     ylab="Densité estimée",
     main="")
tt <- seq(min(θ.chapeau.sim),
          max(θ.chapeau.sim),
          len=501)
lines(tt,dnorm(tt,θ.chapeau,1/sqrt(InfFisher(θ.chapeau))),
      col="red",lwd=2)

Vous construirez le hanging rootogramme vous permettant de comparer votre histogramme avec la densité théorique normale.

On adapte les exemples du cours comme ci-dessous pour obtenir la figure \ref{fig:reponse10}.

mon.histo <- hist(θ.chapeau.sim,breaks=100,plot=FALSE)
ideal <- diff(pnorm(mon.histo$breaks,
                    θ.chapeau,
                    1/sqrt(InfFisher(θ.chapeau))))*n.rep
par(cex=2)
plot(mon.histo$mids,
     sqrt(mon.histo$counts) - sqrt(ideal),
     type="l",lwd=2,
     xlab=expression(hat(theta)),
     ylab="",main="")
abline(h=c(-1,1),col=2,lty=2)
abline(h=0,col="grey80")

On constate que trop de points (de l'ordre de 20) sont en dehors de la bande centrale (entre les lignes pointillées rouges) qui devrait contenir en moyenne 95% des points – le rootogramme étant construit avec 100 bins. La converge de la loi de \(\hat{\theta}\) vers la loi normale n'est donc pas encore parfaite.

Pour chaque expérience simulée, vous construirez un intervalle de confiance à 95% avec la méthode de Wald et vous obtiendrez la probabilité de recouvrement empirique de vos intervalles (c.-à-d. que vous obtiendrez la fraction de vos intervalles qui contiennent la vraie valeur). Commentez (c.-à-d. la méthode de Wald vous fournit-elle des intervalles de confiance à 95% fiables dans ce cas).

dansIC95.Wald <- sapply(1:n.rep,
                        function(i) {
                          K <- sim.exp[,i]
                          θ <- θ.chapeau.sim[i]
                          JF <- faitInfObsF(k1=K[1],
                                            k2=K[2],
                                            k3=K[3])
                          J <- JF(θ)
                          σ <- sqrt(1/J)
                          IC <- θ + c(-1,1)*qnorm(0.975)*σ
                          IC[1] <= θ.chapeau &&
                            θ.chapeau <= IC[2]
                        })
sum(dansIC95.Wald)/length(dansIC95.Wald)

[1] 0,9484

On constate que même si le test du « rootogramme » n'était pas excellent, la probabilité de recouvrement empirique de nos intervalles de confiance à 95% obtenus avec la méthode de Wald est bonne.

La densité exponentielle « complète » est \(p(t) = \nu \exp ( - t \nu)\), la fonction de répartition est donc \[ P(t) = 1 - \exp ( - t \nu) \] et \[ \textrm{Pr}\{T > \delta\} = \exp ( - \delta \nu) \; . \] La densité de probabilité supposée pour le problème est donc : \[ p(t \mid t \ge \delta) = \nu \exp \left( - (t-\delta) \nu\right)\; . \] La log vraisemblance du problème est ainsi : \[ l(\nu) = N \, \log \nu - \nu \sum_i (t_i-\delta) \; , \] où \(N=20\) et \(i=1,\ldots,20\). La fonction score (dérivée première) est : \[ l'(\nu) = N/\nu - \sum_i (t_i - \delta) \] et la dérivée seconde est \[ l''(\nu) = -N/\nu^2\; . \] Cette dernière est toujours négative (tant que \(N>0\), ce qui est le cas ici) et \(l'\) est donc concave et admet un maximum unique donné par : \[ \hat{\nu} = N / \sum_i (t_i - \delta) \; . \] Numériquement il vient :

nu.chapeau <- 1/(mean(intervalles)-delta)

0.815361408944515

D'après le développement de la réponse 1 on a : \[ \hat{\sigma}_{\hat{\nu}} = \hat{\nu} / \sqrt{N} \; , \] soit :

se.nu.chapeau <- nu.chapeau/sqrt(n.intervalles)

0.182320353662994

D'où l'intervalle à 95 % demandé :

round(nu.chapeau+1.96*c(-1,1)*se.nu.chapeau,digits=2)

compte.ll <- function(nu) n.intervalles*((delta-mean(intervalles))*nu+log(nu))
nn <- seq(-3,3,len=101) * se.nu.chapeau + nu.chapeau
par(cex=2)
plot(nn,compte.ll(nn),type="l",lwd=2,
     xlab=expression(nu),
     ylab="log vraisemblance")
lines(nn,compte.ll(nu.chapeau)-(nn-nu.chapeau)^2/2/se.nu.chapeau^2,col=2)
abline(v=-1.96*se.nu.chapeau + nu.chapeau,lty=2)
abline(v=1.96*se.nu.chapeau + nu.chapeau,lty=2)

On voit que l'approximation quadratique commence à être « assez loin » de la fonction approximée. Cela suggère que l'intervalle de confiance que nous venons d'obtenir n'est pas très fiable.

set.seed(225246) 
nu.chapeau.boot.echantillon <- sapply(1:5000,
                                      function(i) 1/(mean(qexp(runif(n.intervalles,
                                                                     pexp(delta,nu.chapeau),1),
                                                               nu.chapeau)
                                                          )-delta
                                                     )
                                      )
round(quantile(nu.chapeau.boot.echantillon,probs=c(0.025,0.975)),digits=2)

On voit que l'intervalle de bootstrap est « décalé » vers la droite par-rapport à l'intervalle obtenu avec l'approximation normale. Ceci est cohérent avec la différence entre l'approximation quadratique et la fonction de log vraisemblance.

La moyenne de l'estimateur de bootstrap est :

mean(nu.chapeau.boot.echantillon)

0.860996385694117

Cet estimateur est :

nu.chapeau.boot.moy <- mean(nu.chapeau.boot.echantillon)
nu.chapeau.boot.se <- sd(nu.chapeau.boot.echantillon)/sqrt(length(nu.chapeau.boot.echantillon))
pnorm((nu.chapeau-nu.chapeau.boot.moy)/nu.chapeau.boot.se)

1.1862179628973e-56

significativement différent de \(\hat{\nu}\). Cela signifie que l'estimateur du maximum de vraisemblance est biaisé pour une taille d'échantillon aussi petite. Avec un échantillon de taille 100 fois plus grance nous aurions :

set.seed(231815)
nu.chapeau.boot.echantillon2 <- sapply(1:1000,
                                       function(i) 1/(mean(round(qexp(runif(n.intervalles*100,
                                                                            pexp(delta,nu.chapeau),
                                                                            1),
                                                                      nu.chapeau
                                                                      ),
                                                                 digits=3)
                                                           )-delta
                                                      )
                                       )
nu.chapeau.boot.moy2 <- mean(nu.chapeau.boot.echantillon2)
nu.chapeau.boot.se2 <- sd(nu.chapeau.boot.echantillon2)/sqrt(length(nu.chapeau.boot.echantillon2))
pnorm((nu.chapeau-nu.chapeau.boot.moy2)/nu.chapeau.boot.se2)

0.183803649654853

Le biais disparaît.

Notons \(N_T=10\) la taille de l'échantillon et \(0 \le m \le N_T\) le nombre d'observations censurées, c.-à-d. le cardinal de l'ensemble \(C\) de l'équation (dans le premier cas, \(m=0\)). Comme les observations sont (supposées) IID nous pouvons réarranger les indices sans perte de généralité, de sorte que les \(N_T-m\) premiers indices correspondent aux observations non censurées. La log vraisemblance peut alors s'écrire : \[ l(\nu) = (N_T-m) \, \log \nu - \nu \, \sum_{i=1}^{N_T-m} t_i - \nu \, m \, T \; . \] La fonction score est donc : \[ l'(\nu) = \frac{N_T-m}{\nu} - \left( \sum_{i=1}^{N_T-m} t_i + m \, T\right) \] et la dérivée seconde est : \[ l''(\nu) = - \frac{N_T-m}{\nu^2} \; . \] Tant que toutes les observations ne sont pas censurées \(l'\) est donc concave et admet un unique maximum. Dans le cas contraire la vraisemblance est maximale en \(\nu = 0\). Comme dans chacun de nos trois cas, l'échantillon contient des éléments non censurés, nous avons : \[ \hat{\nu} = \frac{N_T-m}{N_T \, \bar{t}} \; , \] où \[ \bar{t} = \frac{1}{N_T} \left( \sum_{i=1}^{N_T-m} t_i + m \, T\right) \; . \]

Pour chacun des trois cas nous avons donc :

durées1 <- c(70,11,66,5,20,4,35,40,29,8)
m.cas.1 <- 0
t.bar.cas.1 <- mean(durées1)
durées2 <- c(50,11,50,5,20,4,35,40,29,8)
m.cas.2 <- sum(durées2 == 50)
t.bar.cas.2 <- mean(durées2)
durées3 <- c(25,11,25,5,20,4,25,25,25,8)
m.cas.3 <- sum(durées3 == 25)
t.bar.cas.3 <- mean(durées3)
N.T <- length(durées1)
ν.cas.1 <- (N.T-m.cas.1)/N.T/t.bar.cas.1
ν.cas.2 <- (N.T-m.cas.2)/N.T/t.bar.cas.2
ν.cas.3 <- (N.T-m.cas.3)/N.T/t.bar.cas.3
ν <- c(ν.cas.1,ν.cas.2,ν.cas.3)
data.frame(cas=paste("cas",1:3),
           ν=round(ν,digits=4)
           )

Le développement précédent nous donne : \[ \hat{\sigma}_{\hat{\nu}} = \frac{\nu}{\sqrt{N_T-m}} \; , \] d'où :

σ <- ν / sqrt(N.T-c(m.cas.1,m.cas.2,m.cas.3))
ic.ν <- cbind(-1.96 * σ,+1.96 * σ) + ν
round(ic.ν,digits=4)

mkLLFct <- function(durées,
                      m
                      ) {
    t.bar <- mean(durées)
    N.T <- length(durées)
    facteur <- N.T-m
    function(ν) facteur * log(ν) - ν*N.T*t.bar
  }

  llFct.cas1 <- mkLLFct(durées1,m.cas.1)
  llFct.cas2 <- mkLLFct(durées2,m.cas.2)
  llFct.cas3 <- mkLLFct(durées3,m.cas.3)

  ν.lim <- c(min(ic.ν[,1]),max(ic.ν[,2]))
  νν <- seq(ν.lim[1],ν.lim[2],len=101)
  par(cex=2)
  plot(νν,
       llFct.cas1(νν)-llFct.cas1(ν[1]),
       type="l",
       xlab="ν",
       ylab="log vraisemblance - log vraisemblance maximale",
       lwd=2,
       col=1,
       ylim=c(-6,0))
  lines(νν, - (νν-ν[1])^2/2/σ[1]^2,
        lwd=2,lty=2,col=1)

  lines(νν, llFct.cas2(νν)-llFct.cas2(ν[2]),
        lwd=2,lty=1,col=2)
  lines(νν, - (νν-ν[2])^2/2/σ[2]^2,
        lwd=2,lty=2,col=2)

  lines(νν, llFct.cas3(νν) - llFct.cas3(ν[3]),
        lwd=2,lty=1,col=4)
  lines(νν, - (νν-ν[3])^2/2/σ[3]^2,
        lwd=2,lty=2,col=4)

On constate que l'approximation quadratique n'est jamais bonne. De plus les 3 courbes diffèrent le plus pour les petites valeurs de ν, ce qui est normal car ce sont les grandes valeurs de temps de pannes qui donnent de l'information « de ce coté là » et ce sont les valeurs censurées dans les cas 2 et 3.

On procède comme suit :

set.seed(141718)
n.boot <- 1000
ν.boot.cas.1 <- sapply(1:n.boot,
                       function(i) {
                         durées <- rexp(N.T,ν[1])
                         1/mean(durées)
                       }
                       )
μ.ν.boot.cas.1 <- mean(ν.boot.cas.1)
σ.ν.boot.cas.1 <- sd(ν.boot.cas.1)/sqrt(n.boot)
pnorm((ν[1]-μ.ν.boot.cas.1)/σ.ν.boot.cas.1)

ν.boot.cas.2 <- sapply(1:n.boot,
                       function(i) {
                         durées <- rexp(N.T,ν[2])
                         durées[durées>=50] <- 50
                         m <- sum(durées==50)
                         (N.T-m)/N.T/mean(durées)
                       }
                       )
μ.ν.boot.cas.2 <- mean(ν.boot.cas.2)
σ.ν.boot.cas.2 <- sd(ν.boot.cas.2)/sqrt(n.boot)
pnorm((ν[2]-μ.ν.boot.cas.2)/σ.ν.boot.cas.2)

ν.boot.cas.3 <- sapply(1:n.boot,
                       function(i) {
                         durées <- rexp(N.T,ν[3])
                         durées[durées>=25] <- 25
                         m <- sum(durées==25)
                         (N.T-m)/N.T/mean(durées)
                       }
                       )
μ.ν.boot.cas.3 <- mean(ν.boot.cas.3)
σ.ν.boot.cas.3 <- sd(ν.boot.cas.3)/sqrt(n.boot)
pnorm((ν[3]-μ.ν.boot.cas.3)/σ.ν.boot.cas.3)

ic.ν.boot <- rbind(quantile(ν.boot.cas.1,probs=c(0.025,0.975)),
                   quantile(ν.boot.cas.2,probs=c(0.025,0.975)),
                   quantile(ν.boot.cas.3,probs=c(0.025,0.975))
                   )

round(ic.ν.boot,digits=4)

Une fois les données chargées dans l'espace de travail, il suffit de faire :

par(cex=2)
with(cancer,plot(V2,V1,
                 xlab="Population féminine adulte (1960)",
                 ylab="Nombre de décès par cancer du sein (1950-1960)"))

La façon la plus directe de précéder dans R est d'utiliser la fonction lm. La seule petite subtilité vient du fait qu'on ne veut pas d'ordonnée à l'origine dans le modèle, ce qui nous donne l'appel suivant :

cancer.fit <- lm(V1 ~ V2 - 1, data=cancer)
summary(cancer.fit)

Call:
lm(formula = V1 ~ V2 - 1, data = cancer)

Residuals:
    Min      1Q  Median      3Q     Max 
-65.654  -4.428   0.245   4.779  87.901 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
V2 3.559e-03  4.205e-05   84.62   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.98 on 300 degrees of freedom
Multiple R-squared:  0.9598,	Adjusted R-squared:  0.9597 
F-statistic:  7161 on 1 and 300 DF,  p-value: < 2.2e-16

Cela nous donne directement notre estimation \(\hat{β}\) de β : \(3.56 \times 10^{-3}\).

La méthode plus générale passe par la construction de la matrice du plan d'expérience :

X <- cbind(cancer$V2)

suivie par la décomposition QR et la « récupération » des matrices Q et R :

XdecompQR <- qr(X)
Q <- qr.Q(XdecompQR)
R <- qr.R(XdecompQR)

et l'estimateur est obtenu par :

(beta.chapeau <- backsolve(R, t(Q) %*% cancer$V1))

0.00355876748882502

Si la fonction lm a été utilisée, il suffit d'appeler la fonction plot sur l'objet cancer.fit pour générer les graphes de diagnostique :

plot(cancer.fit)

Ce sont les deux premiers graphes générés par cette commande que je vous ai dit de toujours faire.

Si la « voie QR » a été choisie il fallait faire quelque chose du type :

layout(matrix(c(1:3,3),nc=2,byrow=TRUE))
par(cex=2)
with(cancer,plot(V2,V1,col="blue",xlab="Population",ylab="Décès"))
domaine.pop <- range(cancer$V2)
dd <- seq(domaine.pop[1],domaine.pop[2],len=501)
lines(dd,dd %*% beta.chapeau,col="red",lwd=2)
Y.chapeau <- X %*% beta.chapeau
Z.chapeau <- cancer$V1 - Y.chapeau
qqnorm(Z.chapeau, col="blue",main="",
       xlab="Quantiles théoriques (loi normale)",ylab="Quantiles empiriques")
qqline(Z.chapeau,col="red",lwd=2)
plot(Y.chapeau, Z.chapeau,
     col="blue",xlab="Valeurs prédites",ylab="Résidus")
abline(h=0,col="red",lwd=2)

Ici il y a une dépendance nette de la variance des résidus vis-à-vis de la valeur prédite. Les résidus ne suivent pas non plus une loi normale. Ceci signifie que les hypothèses requises pour construire un intervalle de confiance ayant un sens ne sont pas vérifiées. C'est pourquoi nous allons chercher une transformation stabilisant la variance.

Quand les données sont susceptibles de venir d'une loi de Poisson, la transformation à essayer en priorité est la transformation racine carrée (précisément celle qui est employée dans la construction d'un hanging rootogram). Je donne maintenant la version « longue » passant par la décomposition QR :

Xs <- cbind(sqrt(cancer$V2))
XdecompQRs <- qr(Xs)
Qs <- qr.Q(XdecompQRs)
Rs <- qr.R(XdecompQRs)
(gamma.chapeau <- backsolve(Rs, t(Qs) %*% sqrt(cancer$V1)))    

0.058918594589966

Et les graphes suivent :

layout(matrix(c(1:3,3),nc=2,byrow=TRUE))
par(cex=2)
with(cancer,plot(sqrt(V2),sqrt(V1),col="blue",xlab="Population^1/2",ylab="Décès^1/2"))
domaine.pops <- range(sqrt(cancer$V2))
dds <- seq(domaine.pops[1],domaine.pops[2],len=501)
lines(dds,dds %*% gamma.chapeau,col="red",lwd=2)
Ys.chapeau <- Xs %*% gamma.chapeau
Zs.chapeau <- sqrt(cancer$V1) - Ys.chapeau
qqnorm(Zs.chapeau, col="blue",main="",
       xlab="Quantiles théoriques (loi normale)",ylab="Quantiles empiriques")
qqline(Zs.chapeau,col="red",lwd=2)
plot(Ys.chapeau, Zs.chapeau,
     col="blue",xlab="Valeurs prédites",ylab="Résidus")
abline(h=0,col="red",lwd=2)

La transformation racine carrée a nettement atténué la non homogénéité de la variance. Les résidus sont aussi nettement plus proches d'une loi normale. On peut donc calculer une erreur type reposant sur l'hypothèse d'un bruit homogène suivant une loi normale.

On suit la section 2.2.2.6 du cours pour obtenir :

sigma2.chapeau <- sum(Zs.chapeau^2)/-diff(dim(Xs))
gamma.chapeau.se <- sqrt(sigma2.chapeau*diag(solve(t(Rs) %*% Rs)))
(gamma.chapeau.ci <- qt(c(0.025,0.975),-diff(dim(Xs)))*gamma.chapeau.se+gamma.chapeau)

Soient \(R_1,\ldots,R_N\) N variables aléatoires IID suivant une loi de Rayleigh de paramètre θ. La densité de probabilité conjointe s'écrit : \[p(R_1=r_1,\ldots,R_N=r_n) = \prod_{i=1}^N \frac{r_i}{θ^2} \, \exp\left(\frac{-r_i^2}{2\, θ^2}\right) \; .\] La log-vraisemblance s'écrit donc : \[\mathcal{l}(θ) = -2 N \log θ - \frac{1}{2\, θ^2} \sum_{i=1}^N r_i^2 \, ,\] où le terme \(\sum_{i=1}^N \log r_i\) a été éliminé puisqu'il ne dépend pas de θ. La fonction score est donc : \[S(θ) = -\frac{2 N}{θ} + \frac{1}{θ^3} \sum_{i=1}^N r_i^2 \] et l'information observée : \[\mathcal{J}(θ) = -\frac{2 N}{θ^2} + \frac{3}{θ^4} \sum_{i=1}^N r_i^2 \; .\] L'estimateur du maximum de vraisemblance (EMV), dans les bons cas au moins, est la (ou une des) racine(s) de la fonction score, ce qui nous donne ici : \[\hat{θ} = \sqrt{\frac{\sum_{i=1}^N r_i^2}{2 N}} \; .\] On voit de plus que : \[\mathcal{J}(\hat{θ}) = \frac{8 N^2}{\sum_{i=1}^N r_i^2} > 0\] on a donc bien un maximum puisque l'information observée est l'opposée de la dérivée seconde.

Il suffit d'écrire une fonction qui renvoie par exemple l'EMV, l'erreur type (racine carrée de l'inverse de l'information observée évaluée à l'EMV) et les bornes gauches et droites de l'intervalle de confiance :

toutEnUn <- function(distances) {
    N <- length(distances)
    r2 <- sum(distances^2)
    theta.chapeau <- sqrt(r2/2/N)
    sigma.chapeau <- sqrt(r2/2)/2/N
    gauche <- theta.chapeau - 1.96*sigma.chapeau
    droite <- theta.chapeau + 1.96*sigma.chapeau
    c(emv=theta.chapeau,
      erreurType=sigma.chapeau,
      icGauche=gauche,
      icDroite=droite)}

Nous obtenons donc pour les trois jeux de données :

court.res <- toutEnUn(court)
moyen.res <- toutEnUn(moyen)
long.res <- toutEnUn(long)
cbind(court.res,
      moyen.res,
      long.res)

Le plus simple est de définir une fonction qui génère un graphe comme demandé pour un jeu de données passé comme paramètre :

grapheApproxQuad <- function(distances) {
  stats <- toutEnUn(distances)
  N <- length(distances)
  r2 <- sum(distances^2)
  ll <- function(theta) -2*N*log(theta) - r2/2/theta^2
  tt <- (-200:200)/200*stats[2]*6+stats[1]
  plot(tt,ll(tt),type="l",
       xlab="theta",ylab="log-vraisemblance")
  lines(tt,ll(stats[1])-(tt-stats[1])^2/2/stats[2]^2,
        col=2,lty=2)
}

Et on obtient les figures avec :

layout(matrix(1:3,nc=3))
par(cex=2)
grapheApproxQuad(court)
grapheApproxQuad(moyen)
grapheApproxQuad(long)

On constate que nos fonctions de log-vraisemblance sont légèrement asymétriques et on s'attend donc à ce que nos intervalles basés sur la méthode de Wald soient légèrement faux.

Je vais faire les hanging rootograms en passant, encore une fois, par la définition d'une fonction. Pour faire les choses proprement, j'ai besoin de l'expression de la fonction de répartition d'une loi de Rayleigh : \[F(r|θ) = \int_0^r f(u|θ) du = \int_0^r \frac{u}{θ^2} \exp \left(\frac{-u^2}{2θ^2}\right) du = \int_0^r d\left[-\exp \left(\frac{-u^2}{2θ^2}\right)\right] = 1 - \exp \left(\frac{-r^2}{2θ^2}\right) \; .\]

D'où :

hangingRooto <- function(distances, br=20) {
    stats <- toutEnUn(distances)
    N <- length(distances)
    F <- function(r) 1-exp(-r^2/2/stats[1]^2)
    dist.hist <- hist(distances,breaks=br,plot=FALSE)
    Y <- dist.hist$counts
    X <- dist.hist$mids
    pred <- diff(F(dist.hist$breaks)*N)
    plot(X,sqrt(Y)-sqrt(pred),type="l",
         xlab="Distance (μm)",ylab="")
    abline(h=1.96/2,col=2,lty=2)
    abline(h=-1.96/2,col=2,lty=2)
}

Et on construit les trois figures :

layout(matrix(1:3,nc=3))
par(cex=2)
hangingRooto(court)
hangingRooto(moyen)
hangingRooto(long)

Les ajustements sont plutôt bons !

par(cex=2)
plot(sqrt(c(10,64,192)),
     c(court.res[1],
       moyen.res[1],
       long.res[1]),
     xlab="Racine carrée de la séparation entre marqueurs (Mbp)",
     ylab="Paramètre de la loi de Rayleigh (μm)")
lines(sqrt(c(10,192)),
      c(court.res[1],
       long.res[1]),
      col=2,lty=2)

Le modèle est satisfaisant jusque là !

L'expression de la fonction de répartition établie dans la question précédente nous donne rapidement l'expression de la fonction inverse : \[F^{-1}(u) = θ \sqrt{- 2 \log u} \; ,\]

d'où la réponse à la question en définissant d'abord une fonction générale :

paraBoot <- function(distances, nrep=10000) {
     stats <- toutEnUn(distances)
     N <- length(distances)
     theta <- stats[1]
     boot.rep <- sapply(1:nrep,
                        function(i) {
                            dist.boot <- theta*sqrt(-2*log(runif(N)))
                            stats.boot <- toutEnUn(dist.boot)
                            stats.boot[3] <= theta & theta <= stats.boot[4]})
     c(sum(boot.rep),nrep,sum(boot.rep)/nrep)}

On peut alors obtenir les probabilités de recouvrement pour les trois cas avec :

set.seed(20100928)
cbind(paraBoot(court),
      paraBoot(moyen),
      paraBoot(long))

On voit que les probabilités de recouvrement empiriques sont très proches des probabilités nominales (0,95), malgré la légère asymétrie de la fonction de log-vraisemblance que nous avions constaté précédemment.