# Confidence vs. Credibility Intervals

Tomorrow, for the final lecture of the Mathematical Statistics course, I will try to illustrate – using Monte Carlo simulations – the difference between classical statistics, and the Bayesien approach.

The (simple) way I see it is the following,

• for frequentists, a probability is a measure of the the frequency of repeated events, so the interpretation is that parameters are fixed (but unknown), and data are random
• for Bayesians, a probability is a measure of the degree of certainty about values, so the interpretation is that parameters are random and data are fixed

Or to quote Frequentism and Bayesianism: A Python-driven Primer,  a Bayesian statistician would say “given our observed data, there is a 95% probability that the true value of $\theta$ falls within the credible region” while a Frequentist statistician would say “there is a 95% probability that when I compute a confidence interval from data of this sort, the true value of $\theta$ will fall within it”.

To get more intuition about those quotes, consider a simple problem, with Bernoulli trials, with insurance claims. We want to derive some confidence interval for the probability to claim a loss. There were = 1047 policies. And 159 claims.

Consider the standard (frequentist) confidence interval. What does that mean that $$\overline{x}\pm\sqrt{\frac{\overline{x}(1-\overline{x})}{n}}$$is the (asymptotic) 95% confidence interval? The way I see it is very simple. Let us generate some samples, of size $n$, with the same probability as the empirical one, i.e. $\widehat{\theta}$ (which is the meaning of “from data of this sort”). For each sample, compute the confidence interval with the relationship above. It is a 95% confidence interval because in 95% of the scenarios, the empirical value lies in the confidence interval. From a computation point of view, it is the following idea,

> xbar <- 159
> n <- 1047
> ns <- 100
> M=matrix(rbinom(n*ns,size=1,prob=xbar/n),nrow=n)

I generate 100 samples of size . For each sample, I compute the mean, and the confidence interval, from the previous relationship

> fIC=function(x) mean(x)+c(-1,1)*1.96*sqrt(mean(x)*(1-mean(x)))/sqrt(n)
> IC=t(apply(M,2,fIC))
> MN=apply(M,2,mean)

Then we plot all those confidence intervals. In red when they do not contain the empirical mean

> k=(xbar/n<IC[,1])|(xbar/n>IC[,2])
> plot(MN,1:ns,xlim=range(IC),axes=FALSE,
+ xlab="",ylab="",pch=19,cex=.7,
+ col=c("blue","red")[1+k])
> axis(1)
> segments(IC[,1],1:ns,IC[,2],1:
+ ns,col=c("blue","red")[1+k])
> abline(v=xbar/n) Now, what about the Bayesian credible interval ? Assume that the prior distribution for the probability to claim a loss has a distribution. We’ve seen in the course that, since the Beta distribution is the conjugate of the Bernoulli one, the posterior distribution will also be Beta. More precisely

Based on that property, the confidence interval is based on quantiles of that (posterior) distribution

> u=seq(.1,.2,length=501)
> v=dbeta(u,1+xbar,1+n-xbar)
> plot(u,v,axes=FALSE,type="l")
> I=u<qbeta(.025,1+xbar,1+n-xbar)
> polygon(c(u[I],rev(u[I])),c(v[I],
+ rep(0,sum(I))),col="red",density=30,border=NA)
> I=u>qbeta(.975,1+xbar,1+n-xbar)
> polygon(c(u[I],rev(u[I])),c(v[I],
+ rep(0,sum(I))),col="red",density=30,border=NA)
> axis(1) What does that mean, here, that we have a 95% credible interval. Well, this time, we do not draw using the empirical mean, but some possible probability, based on that posterior distribution (given the observations)

> pk <- rbeta(ns,1+xbar,1+n-xbar)

In green, below, we can visualize the histogram of those values

> hist(pk,prob=TRUE,col="light green",
+ border="white",axes=FALSE,
+ main="",xlab="",ylab="",lwd=3,xlim=c(.12,.18))

And here again, let us generate samples, and compute the empirical probabilities,

> M=matrix(rbinom(n*ns,size=1,prob=rep(pk,
+ each=n)),nrow=n)
> MN=apply(M,2,mean)

Here, there is 95% chance that those empirical means lie in the credible interval, defined using quantiles of the posterior distribution. We can actually visualize all those means : in black the mean used to generate the sample, and then, in blue or red, the averages obtained on those simulated samples,

> abline(v=qbeta(c(.025,.975),1+xbar,1+
+ n-xbar),col="red",lty=2)
> points(pk,seq(1,40,length=ns),pch=19,cex=.7)
> k=(MN<qbeta(.025,1+xbar,1+n-xbar))|
+ (MN>qbeta(.975,1+xbar,1+n-xbar))
> points(MN,seq(1,40,length=ns),
+ pch=19,cex=.7,col=c("blue","red")[1+k])
> segments(MN,seq(1,40,length=ns),
+ pk,seq(1,40,length=ns),col="grey") More details and exemple on Bayesian statistics, seen with the eyes of a (probably) not Bayesian statistician in my slides, from my talk in London, last Summer,

# Puissance et tests statistiques (simples)

Pour comprendre la notion de puissance dans un test statistique, reprenons Statistical power analysis for the behavioral sciences de Jacob Cohen. Comme il le rappelle Fisher (dans The Design of Experiments) ne souhaite pas “prouver” qu’une hypothèse (appelée H_0) est valide, mais il espère, en réalité, rejeter l’hypothèse en question (c’est ce que rappellent d’ailleurs Hubbard et Bayarri en 2003). Autrement dit, supposons que l’on cherche à faire un test médical, et à détecter si une personne est malade (ou pas). Si une personne est généralement malade dès lors qu’un taux de quelque chose est faible (noté ), on va chercher à tester . Dans un test statistique, on espère rejeter l’hypothèse nulle .

The power of a statistical test of a null hypothesis is the probability that it will lead to the rejection of the null hypothesis.

# Tests basés sur la vraisemblance – score

Une autre grandeur intéressant est le score, qui est la dérivée de la vraisemblance. Intuitivement (c’est l’idée de la condition du premier ordre),   et  seront proches si les dérivées en ces points sont proches. En  la dérivée est nulle, donc on va se demander ici, tout simplement, si la dérivée en  est proche de 0. Ou pas. # Tests basés sur la vraisemblance – Rapport de Vraisemblance

Chose promise, chose due. J’avais dit qu’on parlerait du test de rapport de vraisemblance. L’idée – visuelle – est d’avoir une lecture dans l’autre sens : au lieu de se demander si  et  sont proches, on va se demander si   et  sont proches. Si la fonction de vraisemblance est suffisamment régulière, on se pose la même question.

# Tests basés sur la vraisemblance – Wald

Petit rappel préliminaire. Si on dispose d’un échantillon , i.id de loi , où le paramètre  est inconnu, alors on peut calculer le maximum de la vraisemblance, ce qui nous donnera un estimateur intéressant (cf. cours de stat). Illustrons avec un jeu de pile ou face. Reprenons l’échantillon de mon précédent billet,

> X=c(0, 0, 1, 1, 0, 1, 1, 1, 1, 0,
0, 0, 1, 0, 1, 0, 1, 1, 0, 1)

On peut ici tracer la fonction de vraisemblance

ou (un peu plus malin) la fonction de log-vraisemblance

> p=seq(0,1,by=.01)
> logL=function(p)
+ {sum(log(dbinom(X,size=1,prob=p)))}
> plot(p,Vectorize(logL)(p),
+ type="l",col="red",lwd=2)
> p0=.5
> abline(v=p0,col="blue")

La valeur correspondant au trait bleu correspondant au cas que l’on va chercher à tester, de pièce non pipée, soit  (de manière très générale) # Les tests et la logique, modus tollens

Quand on apprend la logique, on apprend la notion de modus tollens, correspondant à une notion de contraposition. Si on a une proposition du genre , alors la proposition contraposée est .  Et on apprend que les deux propositions sont équivalentes (je fais de la logique classique). Par exemple si  correspond à “feu” et  à “fumée“,  signifie que tout feu fait de la fumée. Si cette affirmation est vrai, alors il n’y a pas de fumée sans feu, et donc , autrement dit, s’il n’y a pas de fumée, c’est qu’il n’y a pas de feu.

# Test du chi-deux et indépendance

Considérons le tableau de contingence suivant

> N=margin.table(HairEyeColor, c(1,2))
> N
Eye
Hair    Brown Blue Hazel Green
Black    68   20    15     5
Brown   119   84    54    29
Red      26   17    14    14
Blond     7   94    10    16

avec ici les comptages, que l’on peut aussi traduire sous la forme de probabilités jointes (empiriques)

> N/n
Eye
Hair    Brown  Blue Hazel Green
Black 0.115 0.034 0.025 0.008
Brown 0.201 0.142 0.091 0.049
Red   0.044 0.029 0.024 0.024
Blond 0.012 0.159 0.017 0.027

# Loi multinomiale et loi du chi-deux

La semaine passée, en cours, j’avais rappelé que quand décrivait le compte de  variable multinomiales prenant modalités, la variable

suit asymptotiquement une loi . Et plus généralement, on peut montrer que . Le soucis est que la matrice de variance covariance n’est pas la matrice identité. Pire que ça, elle n’est pas diagonale. Encore pire, elle n’est pas inversible.

On ne peut alors pas utiliser le joli résultat qui nous dit que si , alors suit une loi du chi-deux à  degrés de liberté (dimension du vecteur ). Car il faut que  soit inversible. Or ici, je pense qu’on peut montrer que

(à un facteur près). Et comme on a une loi de probabilité, le terme de droite est nul. Bref, et son inverse n’existe alors pas. Cela dit,  n’est valeur propre qu’une fois, d’où la dimension  qui apparaît sans cesse (à condition que les soient non-nuls, ce qui était une hypothèse dans la définition de la loi multinomiale : toutes les modalités existent vraiment).

L’idée est alors de poser

suit asymptotiquement une loi . Et plus généralement, là encore, .

Ce qui est amusant, c’est que dans le cas binomiale (j’avais dit qu’il était intéressant de voir la loi multinomiale comme une généralisation de la loi binomiale),

puisque

Et en fait, un peut de calculs rapides permettent de noter que

de telle sorte que alors que . Bref, c’est clairement sur cette dernière quantité qu’il faut travailler !

Et là on fait des calculs.

L’astuce est d’utiliser la formule de Stirling

alors pour

soit

et donc

de telle sorte que

On peut alors faire des développements limités, en faisant intervenir

et non pas la version avec au dénominateur. On peut alors réécrire la somme en faisant intervenir termes (et non pas ), exactement comme dans les calculs faits dans le premier cours, avec les calculs de variances. Bref, un peu de calcul supplémentaire permet de montrer que

C’est donc

et non pas

qui sera utilisée comme statistique de test par la suite.

Pour utiliser le test du chi-deux comme test d’ajustement de loi, je pourrais reprendre l’exemple des buts lors des coupes du monde de Football

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
S98=(soccer1998$V1+soccer1998$V2)
S10=(soccer2010$V1+soccer2010$V3)

# Introduction aux tests nonparamétriques

Aujourd’hui, on va commencer la partie du cours sur les tests, avec une introduction par les tests non-paramétriques. On parlera (dans la continuité de ce qu’on a fait jusqu’à présent) de tests d’ajustement de lois, mais aussi de tests plus divers, comme des tests d’indépendance (entre observations successives).

Le test le plus classique est le test des signes. L’idée est de regarder puis de compter le nombre de fois où (différence positive) et le nombre de fois où  (différence négative).

> (X=rnorm(20))
  0.81352495  0.78076709 -0.58654359  0.52805014 -0.31101282 -0.65022034 -0.93619729 -0.16483729 -1.97396753 -0.22506985  1.14959851 -0.38895448 -0.09298421
> (S=c("N","P")[1+1*(diff(X)>0)])
 "N" "N" "P" "N" "N" "N" "P" "N" "P" "P" "N" "P" "P" "N" "P" "N" "P" "P" "N"

On verra ce test en cours. On peut par exemple compter le nombre de fois où la différence était positive,

> sim=function(a=0){
+ n=100
+ epsilon=X=rnorm(n)
+ for(t in 2:n) X[t]=a*X[t-1]+epsilon[t]
+ S=c("N","P")[1+1*(diff(X)>0)]
+ return(sum(S=="P"))}

(je fais ici une fonction un peu plus compliqué pour voir ce qui se passe si les données ne sont pas indépendantes, justement). Si on génère plein d’échantillons, on peut voir que le nombre de variations positives suit une loi normale, si les observations sont effectivement indépendantes,

> library(ggplot2)
> set.seed(1)
> vectS=unlist(lapply(rep(0,1000),sim))
> ggplot(data=data.frame(vectS = vectS),
+ aes(x = vectS)) +
+ geom_histogram(aes(y = ..density..),
+ binwidth = 1, fill = "light blue",
+ col = "black") +
+ stat_function(fun=function(x) dnorm(x,
+ (n-1)/2+.5,sqrt((n+1)/12)),
+ col = "red", size = 1) + ylab("Densite") +
+ ggtitle("Histogramme") # Cross Validation for Kernel Density Estimation

In a post publihed in July, I mentioned the so called the Goldilocks principle, in the context of kermel density estimation, and bandwidth selection. The bandwith should not be too small (the variance would be too large) and it should not be too large (the bias would be too large). Another standard method to select the bandwith, as mentioned this afternoon in class is the cross-validation technique (described in Chiu (1991)). Here, we would like to minimize

The integral can be writen

Since the third component is constant, we have to minimize the expected value of the sum of the first two.

The idea is to approximate it as

which can easily be computed. Consider here some sample, with 50 observations, from a Gaussian distribution,

> set.seed(1)
> X=rnorm(50)

From Silverman’s rule of thumb (which should be appropriate here since the sample has a Gaussian sample) the optimal bandwidth is

> 1.06*sd(X)*length(X)^(-1/5)
 0.4030127

Using the cross-validation technique mentioned above, compute

> J=function(h){
+ fhat=Vectorize(function(x) density(X,from=x,to=x,n=1,bw=h)$y) + fhati=Vectorize(function(i) density(X[-i],from=X[i],to=X[i],n=1,bw=h)$y)
+ F=fhati(1:length(X))
+ return(integrate(function(x) fhat(x)^2,-Inf,Inf)$value-2*mean(F)) + } > vx=seq(.1,1,by=.01) > vy=Vectorize(J)(vx) > df=data.frame(vx,vy) > library(ggplot2) > qplot(vx,vy,geom="line",data=df) The function has the following shape and the optimal value is > optimize(J,interval=c(.1,1))$minimum
 0.4687553

$objective  -0.3355477 Note that, indeed, it is close to Siverman’s optimal bandwidth. # PP plot, QQ plot et estimation à noyau Reprenons ici le code pour récupérer la base de la taille (et du poids) d’un ensemble de 200 personnes. Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] X=Davis$height
n=length(X)

La première chose que l’on avait vu la semaine passée est le PP-plot,

On va ici ‘tester’ une loi normale (de même moyenne et de même variance que notre échantillon… faisons simple)

PP=data.frame(empirique=(1:n)/n,
theorique=pnorm(sort(X),mean(X),sd(X)))

# Fonction de répartition empirique

Ce matin, en cours, on avait utilisé une série d’observations pour décrire l’estimation nonparamétrique de la fonction de répartition

> Davis=read.table(
"http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")

Comme il y une faute de saisie dans la base, on corrige, à la main,

> Davis[12,c(2,3)]=Davis[12,c(3,2)]

Notre variable d’intérêt sera la taille d’un individu,

> X=Davis\$height
> X[1:10]
 182 161 161 177 157 170 167 186 178 171

Pour tracer la fonction de répartition empirique,

on peut utiliser le code suivant

> n=length(X)
> x=seq(140,200)
> ggplot(data=data.frame(x=x), aes(x = x))+
+ stat_function(fun=Vectorize(function(x) sum(X<=x)/n),
+ geom="step", col = "red", size = 1) +
+ ylab("Fonction de Repartition") # Quiz de probabilité

Ce matin, en fin de cours, j’ai fait passé un petit quiz. Compte tenu de la configuration de la salle, j’avais distribué deux questionnaires, avec chacun 10 questions. Les questions étaient identiques, à l’application numérique près (pour quelques exercices). La première série d’exercice est en ligne, ainsi que des éléments de correction; la seconde série est aussi en ligne, avec là encore des éléments de correction. J’ai rajouté après la correction les statistiques de réponses pour chaque formulaire, mais aussi au niveau agrégé (en incluant l’autre questionnaire). Il va de soi qu’il est possible de continuer à s’entraîner à faire des exercices pendant les semaines à venir, ça ne fera de mal à personne !

# Rappels de statistique

Jeudi et vendredi, on finira les rappels de cours, avec quelques rappels de statistique. Les slides sont en ligne,

Sinon, histoire d’illustrer un peu les différents aspects, j’avais dit que j’aimais bien utiliser des simulations. En particulier, il sera intéressant de générer plusieurs (beaucoup) échantillons de taille donnée. Le code R est le suivant (une ligne correspond à un échantillon, de 20 observations)

> n=20
> m=1000
> M=matrix(rnorm(n*m),m,n)

On utilise ensuite la statistique que l’on souhaite étudier, qui peut être la moyenne empirique,

> Xbarn=apply(M,1,mean)

ou la variance

> Sn=apply(M,1,var)

On peut ensuite visualiser la distribution de la statistique, sur tous les échantillons. Par exemple, le fait que l’estimateur de la moyenne est gaussien, on peut faire un histogramme

> hist(Xbarn,probability=TRUE,
+ breaks=seq(-1,1,by=.1),col="light blue")

et rajouter la densité de la loi normale (car la théorie nous dit que la moyenne suit une loi normale),

> u=seq(-2,2,by=.01)
> lines(u,dnorm(u,sd=1/sqrt(n)),col="red",lwd=2) Cela dit, comme Ewen Gallic fera essentiellement du ggplot2 dans son cours de R, je vais aussi mettre les codes en ggplot2 (et essayer de m’y convertir… ce n’est pas l’envie qui manque !)

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   stat_function(fun=function(x) dnorm(x,sd=1/sqrt(n)),
+ col = "red", size = 1) + ylab("Densite") +
+   ggtitle("Histogramme de Xbarn") Cela dit,  au lieu de représenter la loi normale, on peut aussi visualiser un estimateur non paramétrique de la densité (on en reparlera dans quelques semaines), comme l’estimateur à noyau

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   ylab("Densite") +
+   ggtitle("Histogramme de Xbarn") On peut aussi voir que la variance empirique suit une loi du chi-deux,

> hist((n-1)*Sn,probability=TRUE,
+ breaks=seq(0,56,by=2),
+ col="light blue")
> u=seq(0,56,by=.1)
> lines(u,dchisq(u,df=n-1),col="red",lwd=2) avec là encore une version en ggplot2

> ggplot(data = data.frame(x = (n-1)*Sn), aes(x = x)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = 2, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   xlab(expression((n-1) %*% Sn)) +
+   ylab("Densite") +
+   ggtitle(expression(paste("Histogramme de ", (n-1) %*% Sn))) On peut aussi utiliser ces simulations pour comprendre l’idée d’intervalle de confiance à 90%,

> x1=Xbarn-1.64*1/sqrt(n)
> x2=Xbarn+1.64*1/sqrt(n)
> ind=x1*x2>0
> plot(1:m,Xbarn,ylim=c(-1,1),col="white",
+ ylab="intervalle de confiance",xlab="")
> segments(1:m,x1,1:m,x2,col=c("light blue","red")[ind+1],lwd=1.5)
> abline(h=0,lty=2) ou pour la version ggplot2 (mais on reparlera de tout ça en cours, pas forcément cette semaine cela dit)

> ggplot(data = df) +
+   # Les segments. La couleur doit etre fonction d'un logique
+   geom_segment(aes(x = x, xend = x, y = x1, yend = x2, col = x1*x2>0)) +
+   # La ligne horizontale en y = 0
+   geom_hline(yintercept = 0, linetype = "dashed") +
+   # Changer les couleurs des segments
+   scale_colour_manual(values = c("light blue" ,"red")) +
+   # Modifier le titre de l'axe des y
+   ylab("Intervalle de confiance") +
+   # Modifier le titre de l'axe des x
+   xlab("") +
+   # Cacher la legende
+   theme(legend.position="none") On note que dans 90% des scénarios (notion d’intervalle de confiance à 90%), la vraie valeur (connue ici car on simule nos échantillons) est effectivement dans l’intervalle de confiance.

# Rappels de probabilité

Pour les rappels, j’utiliserais a priori des transparents car il y a beaucoup de matière à couvrir.

Je repasserai au tableau une fois les rappels finis.