Category Archives: Probability

Bonferroni, les anniversaires et les terroristes (ou pas)

Le paradoxe des anniversaires est un exercice de probabilité que l’on découvre assez vite. Mais assez souvent, on l’oublie, alors qu’il est en réalité assez profond. Le problème est le suivant: quel est le nombre de personnes qu’il faut réunir pour avoir (au moins) une chance sur 2 pour que deux personnes aient leur anniversaire le même jour. Si on calcule la probabilité que deux personnes dans un groupe de n fêtent leur anniversaire le même jour, on peut calculer la probabilité que personne ne soit né le même jour, parmi n personnes. Un calcul rapide montre qu’avec 23 personnes, la probabilité qu’au moins deux personnes fête son anniversaire le même jour est 50.73%. Autrement dit, il suffit d’une groupe de 23 personnes (c’était la taille d’une classe d’élèves, du moins à mon époque…).

L’intuition, pour comprendre pourquoi nous avons l’impression que ce nombre est petit est que, naturellement, si on a un événement qui survient avec probabilité p, et qu’on répète l’expérience n fois, la probabilité que rien ne se produise est (1-p)n. Aussi, en première approximation (et j’insiste sur ce point parce que c’est avec ce genre d’argument fallacieux qu’on conclue que la probabilité d’avoir une catastrophe nucléaire dépasse 100%, conduisant à du “certitude statistique“), on en déduit que la probabilité qu’un événement survienne est de l’ordre de np. Sauf qu’il, on compte des paires. Et avec n personnes, on constitue environ n2/2 paires. Et donc la probabilité devient proche de n2p/2… ce qui est bien plus grand que np. C’est grossièrement ce qui explique que notre intuition se sent trompée : avec seulement 23 personnes, on peut constituer plus de 250 paires. La probabilité que toutes les paires correspondent à des dates d’anniversaire différente n’est pas si petit qu’il y parait.

Jeff Ullman parle ainsi de principe de Bonferroni,

“if you look in more places for interesting patterns than your amount of data will support, you are bound to find crap”

Et cette intuition faussée quand on travaille sur des paires se retrouve dans les techniques utilisées par les services de renseigngement, tel que le montre Anand Rajaraman dans Mining Massive DatasetsSupposons qu’il existe un milliard d’individus qui pourraient potentiellement être des terroristes. Que les gens, quels qu’il soient vont à l’hôtel 1 jour sur 100. Un hôtel donné peut accueillir 100 personnes. On suppose qu’il existe 100,000 personnes, de sorte qu’en moyenne ils soient plein tous les jours. Et supposons qu’on regarde sur un historique de 3 ans, i.e. 1000 jours pour avoir des calculs ronds. Supposons que pour traquer les terroristes, on cherche les gens qui ont séjourné, ensemble, dans un hôtel, le même jour, au moins deux fois… On imagine que c’est effectivement rare, et donc suspicieux…

La probabilité que deux personnes données dorment à l’hôtel le même jour est faible, 1/100×1/100, soit 0,0001. Et comme il y a 100,000 hôtel, la probabilité que les deux personnes descendent, en plus, dans le même hôtel, il faut diviser 0,0001 par 100,000. Ce qui donne 10-9 . Mais ça, c’est la probabilité pour qu’ils se croisent une fois ! Si on regarde la probabilité qu’ils se croisent deux fois, en trois ans, on arrive à une probabilité de 10-18 Pour ceux qui ont du mal, ça doit faire environ

0,00000000000000001 %

Et encore. Ici il peut en plus se croiser la seconde fois dans un autre hôtel. La probabilité qu’il se croisent dans le même hôtel que la première fois devient encore plus petite !

Supposons maintenant qu’un organisme de surveillance avisé décider de lister les paires de personnes qui, sur des paires de jours, se sont retrouvé dans le même hôtel. Ça doit être suffisamment rare pour être suspicieux. Rappelons que quand on considère n individus, le nombre de paires possibles est de l’ordre de n2/2. Autrement dit, avec un million de personnes, le nombre de paires est de l’ordre de 5×1017. Et pareil, le nombre de paires de jours est de l’ordre de 5×105. Aussi, le nombre d’événements suspicieux que l’on observe, sur ces trois années est de l’ordre de

(5×1017 ) x (5×10) x10-18 = 250000

Autrement dit, 1/4 de millions de paires de personnes sont des terroristes potentiels ! Soit 1/2 millions d’individus.

Maintenant, supposons qu’il y ait 10 terroristes. S’ils sont assez stupides pour se retrouver dans les mêmes hôtels, les mêmes jours, ça veut dire que sur les 500000 personnes qui vont être surveillées, 10 sont effectivement des terroristes, et on a 499990 faux positifs… Enfin, “faux positifs”… ils ont forcément quelque chose d’autre à se reprocher, non ?

Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that  has a similar distribution… so there is no reason to test for a Gamma model for  before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if  has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution  is Gaussian but also the non-conditional distribution of .

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of  is also Gaussian

> library(nortest)
> ad.test(Y)

	Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

	Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here  is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value}
 
> plot(2:100,Vectorize(pv)(2:100),type="l")
> abline(h=.05,col="red")

So here, it could be possible to have also a Gaussian distribution, for . As least to accept that assumption, statistically.

In the context of a Poisson regression, it is well know that it’s not possible to have at the same time  that is Poisson distributed (that’s a Poisson regression) and also  that is Poisson distributed. That simply comes from the fact that

while

and because of the conditional Poisson distribution, then

Thus,

So  cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for .

More generally, it is very difficult to have a distribution family for   that is also the distribution of the non-conditional variable . In the context of a finite mixture ( takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference).

Hence, as a keep saying, before running a regression model on with some given family, it is never a good idea to check if the non-conditional distribution  has the same distribution. Because there is no reason, usually, to remain in the same family.

Collectionner les Figurines Lego avec ses Ami(e)s

Dans mon précédant billet, je regardais combien de figurines Lego je devais acheter pour avoir toute la série, avec ma fille. Pour la série de 16, il fallait, en moyenne, acheter 54 sachets, à 3 euros pièces. Dans les commentaires, @Guillaume et @Gaelle me faisaient noter que les échanges, ça pouvait aider.

Bon, tout d’abord j’aime pas qu’on donne des conseils à mes filles pour avoir raison contre moi ! Parce que oui, c’est sûr qu’en s’y mettant à plusieurs, le coût individuel sera plus faible… Ensuite, parce que les maths me semblent plus compliqués…. Cela dit, en simulant des parties, ça tombe assez vite…

Pour commencer, je fais une petite fonction (histoire d’alléger mes codes) pour incrémenter mes comptages de figurines, à chaque fois que j’en achète une

> f=function(mat,i){
+   mat[i]=mat[i]+1
+   return(mat)
+ }

Ensuite, on peut avoir une petite fonction, correspondant à une simulation, qui me renvoie le nombre de figurine de chaque type que j’ai du acheter pour avoir les 16 une fois. Ou plus. En fait, je veux avoir en tout k fois toutes les figurines (ce qui permet de partager parmi k ami(e)s)

> tirage=function(k=1){
+ test=FALSE
+ cpte=0
+ M=rep(0,16)
+ while(!test){
+   no=sample(1:16,size=1)
+   M=f(M,no)
+   cpte=cpte+1
+   if(min(M)==k) test=TRUE
+ }
+ return(list(compte=M,
+             nbre=cpte))
+ }

Maintenant on peut faire des simulations. Par exemple, pour retrouver le fait qu’il faut, en moyenne, acheter 54 sachet pour que ma fille (seule) ait toutes les figurines

> N=unlist(lapply(1:10000,function(x) tirage(1)$nbre))
> mean(N)
[1] 54.1854

Si maintenant on le fait avec mes deux filles,

> N=unlist(lapply(1:10000,function(x) tirage(2)$nbre))
> mean(N)
[1] 82.7961

ce qui est (beaucoup) moins que le double. Si au lieux de faire ça entre mes deux filles, mais que ma fille collectionne les figurines avec une amie à elle (et que les papas partagent les frais), on a en moyenne 41 sachets à acheter, chacun,

> mean(N)/2
[1] 41.39805

et pour trois, on passe à 36 sachets

> N=unlist(lapply(1:10000,function(x) tirage(3)$nbre))
> mean(N)/3
[1] 36.04363

etc

> V=unlist(lapply(1:10,function(k) mean(unlist(lapply(1:10000,function(x) tirage(k)$nbre)))))
> plot(1:10,V/1:10,type="b",pch=19)

On peut aller plus loin, en regardant les quantiles. Avec 4 ami(e)s, si on achète tous 45 sachets, on a 90% de chances que tout le monde ait les 16 figurines

D’ailleurs, si on regarde l’évolution du quantile (à 90%) sur la moyenne, on voit qu’il baisse vite avec la taille…. Avec 2 ami(e)s, il faut 35% de sachets en plus pour s’assurer d’avoir toutes les figurines dans 90% de cas (par rapport au nombre moyen) alors qu’il en faut juste 20% de plus avec 8 ami(e)s.

Maintenant, si quelqu’un a le courage (et le temps) de faire les maths, je suis preneur !

Collectionner les Figurines Lego

 

Tous ceux qui ont des enfants – ou qui jouent encore au Lego – le savent : régulièrement, Lego sort des séries de figurines, vendus en sachets individuels hermétiques, un peu comme les anciennes cartes à jouer. On est à la saison 15, et il y a 16 figurines par séries,

D’ordinaire, les enfants réclament un paquet, quand on va dans un magasin de jeux (genre pour acheter un cadeau pour un anniversaire) mais là, ma fille les a réclamé tous.

Bon, après le moment où je ris dans ma barbe, j’ai quand même fini par me demander combien ça pourrait coûter si on achetait des paquets, de manière à récupérer les 16 figurines.

La bonne stratégie pour faire le calcul me semble d’ordonner les personnage par ordre d’apparition,  et de modéliser le nombre de paquets achetés entre deux premières apparitions, comme indiqué sur le dessin ci-dessous

Si  désigne le nombre total de paquets à acheter, on note  le nombre de paquets à acheter, entre l’apparition de la ()ème nouvelle figurine, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première figurine), on a

On notera que, si  désigne la loi géométrique,
et à partir de là, les calculs sont simples, puisque

soit (par linéarité de l’espérance)

i.e.

ou encore

Le nombre moyen de paquets à acheter est alors
> sum(16/(1:16)) 
[1] 54.09166

Pour ceux qui ont oublié leur calculatrice, on peut se souvenir que

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)
> log(16)*16 
[1] 44.36142

qui diffère du résultat numérique car il manquait la constante d’Euler

https://latex.codecogs.com/gif.latex?gamma%20=%20lim_{n%20<br /> ightarrow%20infty%20}%20left(%201+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20-%20ln(n)%20<br /> ight)

i.e.

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)+gamma

Numériquement, on obtient, en prenant pour https://latex.codecogs.com/gif.latex?gamma 0.57721 (cf Google)

> (log(16)+ 0.57721)*16 
[1] 53.59678

Mais je m’écarte du sujet… à 3€, il faut débourser, en moyenne, 162€ pour avoir les 16 personnages, soit plus de 10€ (en moyenne) par personnage !

> sum(16/(1:16))*3 
[1] 162.275

On est ici au même prix que l’hotel Lego friends (oui, c’est une demande récurrente)

Par simulation, on peut confirmer cet ordre de grandeur (avec ici un million de simulation de la variable),

> rpaquet=function(i) sum(1+rgeom(16,prob=(1:16)/16))
> S=unlist(lapply(1:1e6,rpaquet))
> mean(3*S)
[1] 162.1709

mais surtout, on peut avoir des quantiles, par exemple

> median(S)*3
[1] 150

après avoir déboursé 150 euros, il y a 50% de chances que l’on ait déjà toutes les figurines.

> quantile(S,.9)*3
90% 
237

 

Et il y a 1 chance sur 10 pour qu’on dépasse les 237€. Mine de rien, pour 16 petites figurines, on dépasse là le prix du Falcon Millenium,

Je pense que les enfants ne sont pas prêt que je leur rachète des petites figurines !

Petit exercice de proba

Mardi dernier, on avait fait une série d’exercices de proba, et je voulais reprendre un exercice pour lequel j’avais proposé une vérification sur ordinateur.

Pour résoudre l’exercice, j’avais suggéré la méthode suivante. est tiré parmi , alors que est tiré parmi .  On note assez facilement que le plus petit nombre est , et comme le ppcm de  et  est ,

(je laisse faire les calculs pour montrer qu’effectivement, au delà, on n’est plus dans les ensembles de départ). Bref, seuls  nombres figurent dans les deux ensembles. La probabilité que l’on cherche est alors donnée par

où on somme sur nos  nombres. On a alors

où chacune des probabilités vaut  (on a une chance sur  de tomber sur n’importe quel nombre), et on a une somme de  termes. Donc la probabilité que l’on cherche est

Pour le vérifier, on peut utiliser le petit code (R) suivant

> list_X=seq(2,by=3,length=100)
> list_Y=seq(3,by=4,length=100)
> n=1e8
> x=sample(list_X,size=n,replace=TRUE)
> y=sample(list_Y,size=n,replace=TRUE)
> sum(x==y)/n
[1] 0.00250221

qui confirme le petit calcul que l’on vient de faire.

The odds of a cluster of airplane accidents

Recently, there have been a lot of airplane accidents.

  • July, 17th 2014, Hrabove, Ukraine, Malaysia Airlines, Boeing 777, fatalities 298 (/298)
  • July, 23rd 2014, Magong, Taiwan, TransAsia Airways, ATR 72-500, fatalities 47 (/58)
  • July, 24th 2014, Aguelhok, Mali, Air Algerie, Mc Donnell Douglas MD-83, fatalities 116 (/116)

It is simple to find a lot of datasets about airplane crashes. For instance on http://ntsb.gov/aviationquery. The dataset is nice, with a lot of information,

> planes=read.table(
+ "cbad3ca6-6b8f-4c98-9ee0-601faAviationData.txt",
+ sep="|",header=TRUE)

for instance the exact location of the crashes,

> library(maps)
> map("world", interior = FALSE)
> points(planes$Longitude,planes$Latitude,
+ pch=19,cex=planes$Total.Fatal.Injuries/50,
+ col="red")

Continue reading The odds of a cluster of airplane accidents

Les anniversaires de vos amis sur Facebook

J’ai découvert avec un peu de retard le joli billet Les anniversaires de vos amis sur Facebook, qui tentait de répondre à la question

Si je possède n amis, quelle est la probabilité qu’il y ait au moins un jour dans l’année où je n’ai pas d’anniversaire à souhaiter ?

Ce problème, on peut aussi l’analyses en posant plutôt la question suivante,

Il y a en tout 365 jours dans une année. Combien d’amis doit-on espérer avoir pour qu’il n’y ait pas un jour sans fêter un anniversaire d’un de ses amis ?

Le problème est différent. Mais la résolution est assez jolie, si on regarde attentivement.

Si je note  le nombre d’amis qu’il me faut pour couvrir tous les anniversaires, et que je commence à regarder, ami après ami le nombre de jours couverts. Supposons que j’ai  jours de l’année couverts, et notons  le nombre de nouveaux amis qu’il va me falloir avant de couvrir  jours de l’année. On peut voir que

Le point intéressant avec cette notation, c’est qu’on connaît la loi de . En fait,  on reconnaît une loi géométrique

Aussi, on peut assez facilement avoir le nombre moyen d’amis qu’il va falloir pour couvrir tous les jours de l’année,

Je laisse les amateurs de calculs essayer de simplifier, mais ça me semble compliqué. Par contre, numériquement on obtient facilement

> sum(365/(1:365))
[1] 2364.646

En fait, la somme, on peut quand même la simplifier un peu, car on reconnaît

Et on sait que

où  est la constante d’Euler. Et effectivement, avec cette approximation, on obtient

> (log(365)+ 0.57721)*365
[1] 2364.144

Ce nombre d’amis est à rapprocher du calcul fait sur https://blogdemaths.wordpress.com/, qui obtenait

on a une chance sur deux pour qu’il y ait au moins un jour sans anniversaire à souhaiter, il faut avoir environ 2285 amis

Calculer un quantile est un peu plus compliqué pour moi (je ne peux calculer que des moyennes), mais on voit que les ordres de grandeurs sont assez proches.

On Hoeffding’s identity

In 1940, Wassily Hoeffding published Masstabinvariante Korrelationstheorie, which was an impressive paper. For those (like me) who unfortunately barely speak German, an English translation could be found in The Collected Works of Wassily Hoeffding, published a few years ago. As I keep saying in my courses about copulas, almost everything was in that paper, by Wassily Hoeffding. For instance, we can see the following graph, of a cumulative distribution function,

What is the difference with a copula? A copula (in dimension 2) is the cumulative distribution function of a random pair with uniform on , as defined by Abe Sklar

But Wassily Hoeffding considered a random pair with uniform on . But everything else is the same. He can even derive the level curves of the density of the Gaussian copula,

> library(mnormt)
> r=.6
> dc=function(u,v) return(
+ as.numeric(dmnorm(cbind(qnorm(u),qnorm(v)),varcov=
+ matrix(c(1,r,r,1),2,2))/dnorm(qnorm(u))/dnorm(qnorm(v))))
> n=500
> vectu=seq(1/n,1-1/n,length=n-1)
> matdc=outer(vectu,vectu,dc)
> contour(vectu,vectu,matdc,levels=
+ c(.325,.944,1.212,1.250,1.290,1.656,3.85),lwd=2)

 

But another interesting point is that there is the so-called Hoeffding’s equality

which is interesting, and quite important, actually, to understand that the covariance (or the correlation) can be seen as some ‘distance‘ to the independence. More precisely, observe that

where  would be the joint cumulative distribution function of some independent variables, with the same marginal distributions.

Of course, it is not exactly a distance, since it can be negative. But still. Now, the thing is that the proof is not trivial. But it is using interesting identities. For instance, in 1885, Franklin wrote a nice paper, Proof of a Theorem of Tchebycheff’s on Definite Integrals, in the American Journal of Mathematics. To get some heuristics about the identity, consider some (finite) sequences  and , then one can prove that

And there is a continuous version of that identity. Consider two bounded functions  and , on some interval,  then

is equal to

In 1979, in Monotone Regression and Covariance Structure, Gerald Shea gave a more probabilistic interpretation of that results, using the fact that

and using a different measure. More precisely, assume now that  functions  and  are integrable, with respect to some measure , on some set . Then

is equal to

In the case where  is a probability measure of , i.e. , this equality is the one used by Wassily Hoeffding, in 1940. The interpretation in terms of random variable is simple that

(with standard assuptions of existence of those quantitites) where  and  are two independent vectors, with identical distribution, . Actually, this relationship can also be found in Some Concepts of Dependence, by E. L. Lehmann, published in 1966. Oh, and by the way, the connection with Chebyshev inequality (claimed in the title of seminal paper by Franklin) come from the fact that if  and  are monotonic, then the left part of the identity is positive, and thus,

But let’s get back to Hoeffiding’s result. How do we get it from that lemma. The idea is to write

as

i.e.

We can then intervert the integral and the expectation, use the fact that

and then, and some integral calculus can be used to rewrite that expression as

So we get here Hoeffding’s identity. Actually, as mentioned by Ben Derrett about the equality above, it can be observed (see http://math.stackexchange.com/105713) that

can also be written

where again,  and  are two independent vectors, with identical distribution, . The later can be writen

Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are  distributions.

  • The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins , and with correlation . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

That was the simple part.

  • The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair , we would like to consider the pair , so that the marginal distributions are . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

  • The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where  and  with independent random variables, given , and that both variables are exponentially distributed, with parameter . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value))
Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value
f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y))
vx=seq(-3,3,length=151)
vy=seq(-3,3,length=151)
z=matrix(NA,length(vx),length(vy))
for(i in 1:length(vx)){
for(j in 1:length(vy)){
z[i,j]=f(vx[i],vy[j])}}
set.seed(1)
Theta=rgamma(1500,a,b)
Z=cbind(rexp(1500,Theta),rexp(1500,Theta))
X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2]))

There is a small technical problem, but no big deal.

Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian.

37% chance

I don’t know if you ever realized, before, but it is quite common to have 37% chance that something happened (or actually “not happened” if we want to be more rigorous). For instance, consider a  grid, and draw  points randomly (and uniformely). Then, around 37% cells are empty. Or if you consider a cell, on that grid, there is 37% chance, that the cell is empty. You can look, on the animation below,

Actually, it is quite simple to prove this result. This come from the fact that

And  is a common probability. To be more specific, we can write it

which is  when  is a Poisson distribution with parameter  (or with mean ). In our case,  is the number of points in a given cell, and since there are   points for   cells (it is a  grid), then . And the Poisson distribution arise here as a limit of a binomial distribution. More precisely, here . We have here the law of small number. Nice isn’t it?

On Wigner’s law (and the semi-circle)

There is something that I love about mathematics: sometimes, you discover – by chance – a law. It has always been there, it might have been well known by some people (specialized in some given field), but you did not know it. And then, you discover it, and you start wondering how comes you never heard about it before… I experienced that feeling this evening, while working on the syallbus for my course on copulas and extreme values. I discovered the so-called Wigner’s Semicircle Law (see e.g. Fan Zhang’s notes, or Fraydoun Rezakhanlou’s notes on that topic). Consider some  random matrice, with  large (say 100) where elements are centered, such as a collection of random variable taking value  with equal probability. Then, eigenvalues can be visualized below

n=100
M=matrix(sample(c(-1,1),size=n*n,replace=TRUE),n,n)
E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))

Consider the symmetric matrix obtained from that matrix,

and more precisely, let us look at its eigenvalues,

E=eigen(.5*(M+t(M)))$values

Then the distribution of those eigenvalues is the so-called semi-circle distribution

hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

Now, if we consider some  distribution, instead of our binomial one, we got exactly the same

M=matrix(rnorm(n*n),n,n)
E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))
E=eigen(.5*(M+t(M)))$values
hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

Actually, it is a very general result, see the second chapter of an Introduction to Random Matrices by Greg Anderson, Alice Guionnet and Ofer Zeitouni, for instance. If entries of the random matrix are independent centred random variables, symmetric, such that higher moments exist, then this property is valid. That’s awesome, isn’t it? Because if the distribution has too heavy tails, then this property is no longer valid. For instance, if we consider a random matrix where entries have a Student distribution, we get something different…

M=matrix(rt(n*n,df=2.1),n,n)
M=M/sd(M)
E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))
E=eigen(.5*(M+t(M)))$values
hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="",
main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65))
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)*2/pi
lines(u,v,col=CL[6],lwd=2)

(here, I do normalize by the standard deviation to get something comparable with the previous graph, where variables were centered, with unit variance)

and if we consider a distribution with infinite variance, we get

M=matrix(rt(n*n,df=1.75),n,n)

we get

I guess I will get back on that property in my course!

Probabilité et géométrie

Une des formules les plus importantes en probabilité (je trouve) est la “formule des probabilités totales” qui dit tout simplement que

que l’ont peut aussi écrire, à l’aide de la formule de Bayes

Une des conséquences de ce résultat est la “law of total expectation“, souvent appelé théorème de double projection,

que l’on écrit souvent sous la forme raccourcie  (dans la formule de droite, le premier symbole est un espérance, c’est à dire une intégrale, et donc un nombre réel, alors que le second indique que l’on travaille sur une espérance conditionnelle, c’est à dire une variable aléatoire). Mais comme toujours dans les relations simplifiée, il faut savoir de quoi on parle.

La démonstration est simple

soit

ou enfin

tout simplement.

L’interprétation est encore plus simple (pour utiliser un exemple que j’aime de plus en plus) : si est le poids d’une personne https://latex.codecogs.com/gif.latex%20?Y, et https://latex.codecogs.com/gif.latex%20?X le sexe, alors

https://latex.codecogs.com/gif.latex%20?\mathbb{E}(Y)=%20\mathbb{E}(Y\vert%20X=\text{H})\cdot\mathbb{P}(X=\text{H})+\mathbb{E}(Y\vert%20X=\text{F})\cdot\mathbb{P}(X=\text{F})

c’est à dire que le poids moyen d’une personne prise au hasard est un barycentre (damned, je fais déjà de la géométrie) entre le poids moyen des hommes et le poids moyen des femmes, les poids pour le calcul du barycentre étant liés aux proportions d’hommes et de femmes. Tout simplement.

Il existe une autre relation classique en probabilité (connue en statistiques sous le nom de “formule de décomposition de la variance”) qui dit que

Là aussi, on a une écriture un peu simplificatrice, que l’on va essayer de décortiquer un peu. Pour ça, le plus simple est de faire un peu de géométrie. Car oui, dans les espaces de variables aléatoire, on peut faire de la géométrique. En particulier des projections (orthogonales). Mais avant, de faire des projections, il faut des distances, des angles, une notion d’orthogonalité, etc.

  • Rappels de géométrie

Pour parler d’orthogonalité, il faut une notion de produit scalaire. Pour rappel, un produit scalaire, sur un espace https://latex.codecogs.com/gif.latex?\mathcal%20E, c’est défini par

  • une notion de symétrie https://latex.codecogs.com/gif.latex?%3C\vec%20x,\vec%20y%3E=%3C\vec%20y,\vec%20x%3E
  • de bilinéarité https://latex.codecogs.com/gif.latex?%3C\alpha%20\vec%20x+\beta%20\vec%20y,\vec%20z%3E=\alpha%3C\vec%20x,\vec%20z%3E+\beta%3C\vec%20y,\vec%20z%3E
  • de positivité, https://latex.codecogs.com/gif.latex?%3C\vec%20x,\vec%20x%3E%20\geq%200
  • et le produit scalaire est dit défini https://latex.codecogs.com/gif.latex?%3C\vec%20x,\vec%20x%3E=0 implique https://latex.codecogs.com/gif.latex?\vec%20x=\vec%200

De ce produit scalaire, on peut déduire une norme

https://latex.codecogs.com/gif.latex?\|\vec%20x\|=\sqrt{%3C\vec%20x,\vec%20x%3E}

En effet, on retrouve la propriété d’homogénéité, https://latex.codecogs.com/gif.latex?\|\lambda\vec%20x\|=\vert\lambda\vert%20\|\vec%20x\| (c’est pour ça que l’on prend la racine carrée), et https://latex.codecogs.com/gif.latex?\|\vec%20x\|=0 si et seulement si https://latex.codecogs.com/gif.latex?\vec%20x=\vec%200. On a aussi l’inégalité de Cauchy-Schwarz,

https://latex.codecogs.com/gif.latex?\|%3C\vec%20x,\vec%20y%3E\|%20\leq%20\|\vec%20x\|\|\vec%20y\|

qui va impliquer l’inégalité triangulaire

https://latex.codecogs.com/gif.latex?\|\vec%20x+y\|\leq%20\|\vec%20x\|+\|\vec%20y\|

On dira que deux vecteurs sont orthogonaux, noté https://latex.codecogs.com/gif.latex?\vec%20x\perp\vec%20y si https://latex.codecogs.com/gif.latex?%3C\vec%20x,\vec%20y%3E=0. Dans ce cas, on a la relation de Pythagore,

https://latex.codecogs.com/gif.latex?\|\vec%20x+y\|^2=%20\|\vec%20x\|^2+\|\vec%20y\|^2

Autre objet intéressant, la projection orthogonale. Rappelons tout d’abord que si https://latex.codecogs.com/gif.latex?\mathcal%20F\subset%20\mathcal%20E, on dira que https://latex.codecogs.com/gif.latex?\vec%20x\perp\mathcal%20F si https://latex.codecogs.com/gif.latex?\vec%20x\perp\vec%20yy pour tout https://latex.codecogs.com/gif.latex?\vec%20y\in\mathcal%20F. Et on notera

https://latex.codecogs.com/gif.latex?\mathcal%20F^\perp=\{\vec%20x\in\mathcal%20E;\vec%20x\perp\mathcal%20F\}

On a alors des résultats intéressants. En particulier, pour tout https://latex.codecogs.com/gif.latex?\vec%20y\in\mathcal%20E, il existe un unique https://latex.codecogs.com/gif.latex?\vec%20y_\star\in\mathcal%20F\subset%20\mathcal%20E tel que

https://latex.codecogs.com/gif.latex?\|\vec%20y-\vec%20y_\star\|=%20\inf\{\|\vec%20y-\vec%20z\|%20,\vec%20z\in\mathcal%20F\}

On peut aussi montrer que https://latex.codecogs.com/gif.latex?\vec%20y-\vec%20y_\star%20\in\mathcal%20F^\perp. On parlera alors de projection orthogonale de https://latex.codecogs.com/gif.latex?\vec%20y sur https://latex.codecogs.com/gif.latex?\mathcal%20F, et on pourra noter https://latex.codecogs.com/gif.latex?\vec%20y_\star=\Pi_{\mathcal%20F}(\vec%20y). Tout cela est assez standard dans les espaces de dimension finie (on pensera à https://latex.codecogs.com/gif.latex?\mathbb{R}^n pour avoir un peu d’intuition). On connait plein de choses sur les projections orthogonales, il suffit de penser dans https://latex.codecogs.com/gif.latex%20?\mathbb{R}^n pour avoir un peu d’intuition. Par exemple, si , alors

Cette relation est vraie dans https://latex.codecogs.com/gif.latex%20?\mathbb{R}^n, mais aussi dans des espaces plus généraux. C’est ce qu’on appelle le théoreme de double projection (on qu’on devrait voir réapparaître sur les variables aléatoires).

  • Géométrie dans les espaces de variables aléatoires

L’espace des variables aléatoires de variance finie – https://latex.codecogs.com/gif.latex?L_2 – peut être muni d’une telle opération, de produit scalaire,

https://latex.codecogs.com/gif.latex?%3CX,Y%3E=\mathbb{E}(XY)

La norme est alors https://latex.codecogs.com/gif.latex?\|X\|=\sqrt{\mathbb{E}(X^2)}. En fait, on voit qu’on s’avance sur un terrain glissant ici, car pour avoir un produit scalaire, il faudrait que https://latex.codecogs.com/gif.latex?%3CX,X%3E=0 si et seulement si https://latex.codecogs.com/gif.latex?X=0, et pour avoir une norme il faudrait que https://latex.codecogs.com/gif.latex?\|X\|=0 si et seulement si https://latex.codecogs.com/gif.latex?X=0. Le soucis technique ici est que https://latex.codecogs.com/gif.latex?\mathbb{E}(X^2)=0 signifie que https://latex.codecogs.com/gif.latex?\mathbb{P}(X=0)=1 et pas https://latex.codecogs.com/gif.latex?X(\omega)=0 pour tout https://latex.codecogs.com/gif.latex?\omega. Bref, l’égalité est a comprendre au sens presque partout. Mais c’est un point de détail (ici). Techniquement, comme l’explique justement Williams (1991) dans le chapitre 6, https://latex.codecogs.com/gif.latex?L_2 est précisément l’espace  (classique en probabilité) quotienté par cette relation d’équivalence .

Un sous-espace de https://latex.codecogs.com/gif.latex?L_2 est l’espace des constante, https://latex.codecogs.com/gif.latex?\mathbb{R}, que je noterais aussi https://latex.codecogs.com/gif.latex?s\{\boldsymbol{1}\}. On notera que, comme rappelé dans tous les cours que je donne

https://latex.codecogs.com/gif.latex?\mathbb{E}(Y)%20=\underset{c\in\mathbb{R}}{\text{argmin}}\{\mathbb{E}([Y-c]^2)\}

c’est à dire que

https://latex.codecogs.com/gif.latex?\mathbb{E}(Y)=\Pi_{s\{\boldsymbol{1}\}}(Y)

L’espérance est la projection orthogonale sur l’espace des constantes. Et

https://latex.codecogs.com/gif.latex?\text{Var}(Y)%20=\underset{c\in\mathbb{R}}{\text{min}}\{\mathbb{E}([Y-c]^2)\}

Un autre sous-espace de https://latex.codecogs.com/gif.latex?L_2 qui sera intéressant est le suivant. Si https://latex.codecogs.com/gif.latex?X\in%20L_2, on notera

https://latex.codecogs.com/gif.latex?s\{X\}=\{\tilde%20X=\psi(X)\in%20L_2,\psi:\mathbb{R}\rightarrow\mathbb{R}\}

Par exemple, on peut considérer une variable de Bernoulli, https://latex.codecogs.com/gif.latex?X=\boldsymbol{1}_Ahttps://latex.codecogs.com/gif.latex?A\subset%20\mathbb{R}. Dans ce cas, le sous-espace est équivalent au sous-espace des combinaisons linéaires,

https://latex.codecogs.com/gif.latex?sl\{X\}=\{\tilde%20X=\beta_0+\beta_1%20X;(\beta_0,\beta_1)\in%20\mathbb{R}^2\}\subset%20L_2

Techniquement, je pense que

https://latex.codecogs.com/gif.latex?\overline{\sigma\{X\}}%20\subset%20s\{X\}

mais il semblerait qu’on ait égalité stricte… Mais comme je veux parler de projection sans parler de l’espace https://latex.codecogs.com/gif.latex?\Omega, on va continuer avec mon interprétation heuristique, et je ne parlerais pas de https://latex.codecogs.com/gif.latex?\sigma\{X\}. Je fais des projections sur des sous-espaces de https://latex.codecogs.com/gif.latex?L_2, pas sur des https://latex.codecogs.com/gif.latex?\sigma-algèbres.

Si https://latex.codecogs.com/gif.latex?Y\in%20L_2, on posera

https://latex.codecogs.com/gif.latex?\mathbb{E}(Y\vert%20X)=\Pi_{s\{\boldsymbol{X}\}}(Y)

qui sera la solution du problème de moindres carrés

https://latex.codecogs.com/gif.latex?\mathbb{E}(Y\vert%20X)%20=\underset{X_\star\in%20s\{X\}}{\text{argmin}}\{\mathbb{E}([Y-X_\star]^2)\}

On notera que https://latex.codecogs.com/gif.latex?\mathbb{E}(Y\vert%20X)%20\in%20s\{X\}}, et on retrouve ici l’écriture standard d’un modèle de régression, https://latex.codecogs.com/gif.latex?\mathbb{E}(Y\vert%20X)%20=\psi(X). On parlera alors d’espérance conditionnelle. Qui est ici une variable aléatoire, par construction. En fait, l’unicité de la variable aléatoire est possible précisément parce que tout a l’heure, on définissait une égalité presque sure.

On avait rappelé le théor`eme de double projection tout a l’heure. Et comme https://latex.codecogs.com/gif.latex%20?s\{\boldsymbol{1}\}\subset%20s\{X\}, alors

https://latex.codecogs.com/gif.latex%20?\Pi_{s\{\boldsymbol{1}\}}(Y)=\Pi_{s\{\boldsymbol{1}\}}\big(\Pi_{s\{X\}}(Y)\big)

ce qui se traduit par la relation bien connue

https://latex.codecogs.com/gif.latex%20?\mathbb{E}(Y)=%20\mathbb{E}\big(\mathbb{E}(Y\vert%20X)\big)

En économétrie (linéaire) et dans les premiers cours de séries temporelles, on introduit un autre type d’opérateur car on ne projette pas dans des espaces aussi gros que https://latex.codecogs.com/gif.latex%20?s\{X\}. On introduit souvent l’opérateur d’espérance linéaire, avec

https://latex.codecogs.com/gif.latex?\text{EL}(Y\vert%20X)=\Pi_{sl\{{X}\}}(Y)

où (comme on l’avait introduit tout à l’heure)

https://latex.codecogs.com/gif.latex?sl\{X\}=\{\tilde%20X=\beta_0+\beta_1%20X;(\beta_0,\beta_1)\in%20\mathbb{R}^2\}\subset%20L_2

et pour les séries temporelles stationnaires,

https://latex.codecogs.com/gif.latex?\text{EL}(X_t\vert%20X_{t-1},X_{t-2},\cdots)=\Pi_{\overline{sl\{\boldsymbol{X}\}}}(Y)

Oui, si le processus est stationnaire, je peux régresser sur tout le passé, mais dans ce cas, je pense qu’il faut fermer l’espace afin de projeter dessus…

Bon, maintenant, on voit que mon interprétation géométrique est bancale… en effet, je ne peux définir mon espérance conditionnelle (voire mon espérance tout court) à l’aide de projections orthogonales que si je suis dans un espace muni qu’un tel opérateur. Et si https://latex.codecogs.com/gif.latex?L_2 est un espace de Hilbert, ce n’est pas le cas de https://latex.codecogs.com/gif.latex?L_1. Pourtant, nul besoin d’être dans https://latex.codecogs.com/gif.latex?L_2 pour définir une espérance conditionnelle… Je pense que l’astuce peut être de noter que si https://latex.codecogs.com/gif.latex?X appartient à https://latex.codecogs.com/gif.latex?L_1 mais pas https://latex.codecogs.com/gif.latex?L_2, on peut quand même s’en sortir car https://latex.codecogs.com/gif.latex?L_1 est un espace de Banach et on a malgré tout des notions de convergences. Et https://latex.codecogs.com/gif.latex?L_2 est dense dans https://latex.codecogs.com/gif.latex?L_1. Bref, en faisant un peu de limite, pour peut étendre ce qu’on vient de faire pour les variables qui ne sont pas de carré intégrable…

Visuellement (il serait peut-être temps de faire des dessins, non ?) on a

On voit apparaître un paquet de relations sur ce dessin, que ce soit des doubles projections, mais aussi le théorème de Pythagore. Mais prenons notre temps…

  • et la constante ?

Maintenant, si on regarde un peu détails, il y a des choses étranges. Par exemple https://latex.codecogs.com/gif.latex?X\perp%20Y se traduit ici par https://latex.codecogs.com/gif.latex?\mathbb{E}(XY)=0. Mais ce n’est pas ce qu’on utilise classiquement comme notion d’orthogonalité. On a davantage l’habitude de voir

https://latex.codecogs.com/gif.latex?\text{Cov}(X,Y)=\mathbb{E}(XY)-\mathbb{E}(X)\mathbb{E}(Y)=0

Qu’est ce qu’on a raté ? L’idée est de poser https://latex.codecogs.com/gif.latex%20?\vec%20x=X-\mathbb{E}(X), et https://latex.codecogs.com/gif.latex%20?\vec%20y=Y-\mathbb{E}(Y), c’est à dire que l’on va travailler sur les variables centrées. On peut noter que

https://latex.codecogs.com/gif.latex%20?\text{Var}(X)=\mathbb{E}([X-\mathbb{E}(X)]^2)=\|%20\vec%20x%20\|^2

et on va définir la covariance entre X et Y a partir des vecteurs translatés,

https://latex.codecogs.com/gif.latex%20?\text{Cov}(X,Y)=%3C\vec%20x,\vec%20y%3E

et la corrélation comme

https://latex.codecogs.com/gif.latex%20?\text{corr}(X,Y)=\frac{%3C\vec%20x,\vec%20y%3E}{\|%20\vec%20x%20\|\cdot%20\|%20\vec%20y%20\|}%20=\cos(\theta)

où https://latex.codecogs.com/gif.latex%20?\theta est l’angle entre les vecteurs https://latex.codecogs.com/gif.latex%20?\vec%20x et https://latex.codecogs.com/gif.latex%20?\vec%20y. On dira que https://latex.codecogs.com/gif.latex?X\perp%20Y si https://latex.codecogs.com/gif.latex%20?\vec%20x%20\perp%20\vec%20y (au sens géométrique du terme), qui correspond a la relation classique https://latex.codecogs.com/gif.latex%20?\text{corr}(X,Y)=0). Le théroreme de Pythagore nous dit que si https://latex.codecogs.com/gif.latex%20?\vec%20x%20\perp%20\vec%20y, alors

https://latex.codecogs.com/gif.latex%20?\|%20\vec%20x%20+\vec%20y%20\|^2=\|%20\vec%20x%20\|^2+\|%20\vec%20y%20\|^2

Ce qui se traduit par

https://latex.codecogs.com/gif.latex%20?\text{Var}(X+Y)=\text{Var}(X)+\text{Var}(Y)

si les variables sont orthogonales.

Maintenant, on était allé un peu plus loin tout a l’heure, avec des variables correspondant à des espérances conditionnelles. Par exemple, la formule de décomposition de la variance. On peut reprendre l’égalité du théorème de Pythagore, avec

https://latex.codecogs.com/gif.latex%20?\vec%20x=\mathbb{E}(Y\vert%20X)-\mathbb{E}(Y)

et

https://latex.codecogs.com/gif.latex%20?\vec%20y=Y-\mathbb{E}(Y\vert%20X)

On note facilement que https://latex.codecogs.com/gif.latex%20?\vec%20x%20\perp%20\vec%20y, car https://latex.codecogs.com/gif.latex%20?\vec%20x\in%20s\{X\} alors que https://latex.codecogs.com/gif.latex%20?\vec%20y\in%20s\{X\}^\perp, par construction de la projection orthogonale. Aussi le théorème de Pythagore nous dit que

https://latex.codecogs.com/gif.latex%20?\|%20Y-\mathbb{E}(Y)%20\|^2=\|%20\mathbb{E}(Y\vert%20X)-\mathbb{E}(Y)%20\|^2+\|%20%20Y-\mathbb{E}(Y\vert%20X)%20\|^2

Pour le terme de gauche, c’est assez facile,

https://latex.codecogs.com/gif.latex%20?\|%20Y-\mathbb{E}(Y)%20\|^2=%20\text{Var}(Y)

Pour le premier terme de droite, la aussi, c’est facile, car https://latex.codecogs.com/gif.latex%20?\mathbb{E}[\mathbb{E}(Y\vert%20X)]=\mathbb{E}(Y), par le théorème de double projection, et donc

https://latex.codecogs.com/gif.latex%20?\|%20\mathbb{E}(Y\vert%20X)-\mathbb{E}(Y)%20\|^2=\text{Var}(\mathbb{E}(Y\vert%20X))

Pour le dernier terme, a droite, c’est plus vicieux. Disons qu’on va l’identifier a une variance conditionnelle (je n’ai toujours pas défini formellement de variable aléatoire appelée “variance conditionnelle“), i.e.

https://latex.codecogs.com/gif.latex%20?\|%20\vec%20Y-\mathbb{E}(Y\vert%20X)%20\|^2%20=%20\mathbb{E}(\text{Var}(Y\vert%20X))

On retrouve alors la formule de décomposition de la variance, en variance intra et variance inter.

Cette formule est classique en statistique, en régression (les résidus sont orthogonaux a https://latex.codecogs.com/gif.latex%20?sl\{\boldsymbol{X}\}, on aura donc une partie de la variance qui sera expliquée par nos variables explicatives https://latex.codecogs.com/gif.latex%20?\boldsymbol{X} et une partie qui sera dite non expliquée) ou en crédibilité (je peux renvoyer cette fois à l’exemple des jeux de fléchettes de Philbrick (1982)). Promis, on reparlera de géométrie cet hiver, quand je ferais les lois multivariées en cours (avec cette fois des notions d’invariance, de rotations, de symétries, etc)

Generating functions

Today, I wanted to publish a post on generating functions, based on discussions I had with Jean-Francois while having our coffee after lunch a couple of times already. The other reason is that I publish my post while my student just finished their Probability exam (and there were a few questions on generating functions).

  • A short introduction (back on a specific exercise)

In the Probability exam, I included an exercise we’ve seen in class, last week. The question is the following (question 16 in the form – in French). Let  for  and  for  be the cumulative distribution function of some random variable , i.e. . What is the moment generating function of , i.e.  ?

Consider some  (we’ll see later on if some additional constraint are necessary). The tricky part of this exercice appears extremely fast, actually: how could you write  ? I mean, in any probability textbook, the standard answer is

  • if  is discrete,

  • if  is (absolutely) continuous,

where  is the density of . Here,  is clearly not a discrete variable. But is it (absolutely) continuous. My (strong) belief is that you need to plot that distribution function to see how it looks like, , for all 

(following recent discussions with Philippe Reka, I will try to post more hand-made graphs)

Ooops. It looks like we have a discontinuity in 0. So we have to be a bit carefull here :  is neither continuous nor discrete. Let us use the double projection formula,

which can also be writen, if ,

This is simply the idea of saying that the overall average is a barycenter of the average per subgroup. Here, and let  while  (note that ). Thus,

Let us consider the three different components.

 

and

(since it is is a real-valued constant), and here . So finally, we should compute . Observe that  given  is a (absolutely) continuous random variable, with a density. To get it, observe that for all ,

and , i.e.  given  is an exponential distribution.

Hence,  is a mixture between an exponential variable and a Dirac mass in . This was actually the tricky part of the question since it is not obvious when we see (only) the formula above.

From now on, it is just high-school level computations,

if  (for the first time, we see that the function is not defined everywhere). If we put all the expressions together,

  • Monte Carlo computations

If we are lazy (and trust me, I am extremely lazy), it is possible to use Monte Carlo simulations to compute that function,

> F=function(x) ifelse(x<0,0,1-exp(-x)/3)
> Finv=function(u) uniroot(function(x) F(x)-u,c(-1e-9,1e4))$root

or (to avoid the problem of the discontinuity)

> Finv=function(u) ifelse(3*u>1,0,uniroot(function(x)
+ F(x)-u,c(-1e-9,1e4))$root))

Here, the inverse is simple to get, so we can faster the code using

> Finv=function(u) ifelse(3*u>1,0,-log(3*u))

Then, we use

> rF=function(n) Vectorize(Finv)(runif(n))
> M=function(t,n=10000) mean(exp(t*rF(n)))
> Mtheo=function(t) (3-2*t)/(3-3*t)
> u=seq(-2,1 ,by=.1)
> v=Vectorize(M)(u)
> plot(u,v,type="b",col='blue')
> lines(u,Mtheo(u),col="red")

The problem with Monte Carlo simulations is that they should be used only if they are valid. If mean, I can compute

> set.seed(1)
> M(3)
[1] 5748134

Finite sum can always be computed, numerically. Even if here, https://latex.codecogs.com/gif.latex?\mathbb{E}(e^{3X}) does not exist (or to be more precise, is not finite). It is like the average of a Cauhy sample… I can always compute it, even if the expected value does not exists…

> set.seed(1)
> mean(rcauchy(1000000))
[1] 0.006069028

This is related to questions I tried to ask a few years ago in a paper, where I wanted to test if  (or not). Almost all the tests I know are actually based on that assumption… But this is not the point here. My point is that those generating functions are interesting, when then exist. And perhaps working with characteristic function is a better idea.

  • Generating functions

Now, to get back on the begining of last course, generating functions are interesting for a lot of reasons. But first of all, let us define those function properly.

The moment generating function  exists if it is finite on a neighbourhood of https://latex.codecogs.com/gif.latex?0 (there is an https://latex.codecogs.com/gif.latex?a%3E0 such that for all https://latex.codecogs.com/gif.latex?t\in[-a,+a], https://latex.codecogs.com/gif.latex?M_X(t)%3C\infty). In that case, there exists some (open) interval https://latex.codecogs.com/gif.latex?(a,b)\in\overline{R} such that for all https://latex.codecogs.com/gif.latex?t\in(a,b), https://latex.codecogs.com/gif.latex?M_X(t)%3C\infty, called the convergence strip of the moment generating function.

This function is said to be moment generating, since if https://latex.codecogs.com/gif.latex?M_X(\cdot) exists (as defined in the previous paragraph), then all moments exist, for all https://latex.codecogs.com/gif.latex?k\in%20\mathbb{N}\backslash\{0\}, https://latex.codecogs.com/gif.latex?\mathbb{E}\left(\vert%20X\vert^k\right)%3C\infty. This is basically due to the fact that, for all https://latex.codecogs.com/gif.latex?k\in%20\mathbb{N}\backslash\{0\}https://latex.codecogs.com/gif.latex?x^k\exp(-\vert%20t\vert%20x)\rightarrow%200 as https://latex.codecogs.com/gif.latex?x\rightarrow\infty, so, for all https://latex.codecogs.com/gif.latex?x large enough, https://latex.codecogs.com/gif.latex?x^k%20\leq%20\exp(\vert%20t\vert%20x). And before, it is always possible to use a multiplicative constant,

for some https://latex.codecogs.com/gif.latex?K. Thus,

if https://latex.codecogs.com/gif.latex?t is small enough (namely https://latex.codecogs.com/gif.latex?[-t,+t] belongs to the convergence strip).

Now, if we use Taylor’s expansion,

and

If we look at the value of the derivative of that function at point 0, then

As we’ve seen last week in class, it is possible to define a moment generating function in higher dimension, for some random vector https://latex.codecogs.com/gif.latex?\boldsymbol{X}=(X_1,\cdots,X_d),

for some https://latex.codecogs.com/gif.latex?\boldsymbol{t}\in\mathbb{R}^d. It is again a moment generating function since crossed derivatives (taken a point https://latex.codecogs.com/gif.latex?\boldsymbol{0}) are cross-moments. For instance,

 

Some, moment generating functions are interesting if you want to derive moments of a given distribution. Another interesting feature is that this moment generating function (under certain conditions) fully characterize the distribution of the random variable, in the sense that if for some https://latex.codecogs.com/gif.latex?%20h%3E0,
https://latex.codecogs.com/gif.latex?%20M_X(t)=M_Y(t) for all https://latex.codecogs.com/gif.latex?%20t\in(-h,+h), then https://latex.codecogs.com/gif.latex?X\overset{\mathcal{L}}{=}Y.

  • From moment generating functions to characteristic functions

The problem with the moment generating function is that the function is defined (only) on some neighborhood of https://latex.codecogs.com/gif.latex?%200, and we should be careful. The other problem is that it does exist only for distribution in https://latex.codecogs.com/gif.latex?%20L_\infty. Which might be a strong assumption.

Thus, an interesting idea is to consider https://latex.codecogs.com/gif.latex?%20\mathbb{E}\left(%20e^{tX}%20\right) not on the real line, but on the imaginary line.

Thus, let https://latex.codecogs.com/gif.latex?%20\phi_X(t)=\mathbb{E}\left(%20e^{i%20tX}%20\right) for some https://latex.codecogs.com/gif.latex?%20t\in\mathbb{R}. Actually, not some, but all https://latex.codecogs.com/gif.latex?%20t\in\mathbb{R}, since

so the characteristic function always exists. Paul Lévy proved in 1925 that the characteristic function completely characterizes the distribution.

Now, if we look at it quickly, it looks like we did not change a lot of things here, and we should be able to write

https://latex.codecogs.com/gif.latex?%20\phi_X(t)=M_X(i%20t)

If we want to do things properly, let us look at Gut (2005) for instance. Assume that https://latex.codecogs.com/gif.latex?%20M_X(\cdot) is defined on some interval https://latex.codecogs.com/gif.latex?%20(-a,+a). It is then possible to define a function  (this time, it is no longer a real-valued function) as

which is well defined on some strip .
https://latex.codecogs.com/gif.latex?%20\phi_X(\cdot)and https://latex.codecogs.com/gif.latex?%20M_X(\cdot) are then restriction of that function respectively on the imaginary line, and the real line. That function https://latex.codecogs.com/gif.latex?%20\Gamma_X(\cdot) is clearly holomorphic, and thus, the value it takes on such a strip is fully determined by the values it takes on the real interval https://latex.codecogs.com/gif.latex?%20(-a,+a). Thus, the moment generating function will completely characterize the distribution.

But it has to be defined on some neighbourhood of https://latex.codecogs.com/gif.latex?%200. Which is not trivial actually… I mean, I nonlife insurance, we see a lot a Pareto distributions.

  • Fast Fourier Transform

Recall Euler’s formula,

Thus, we should not be surprised to see Fourier’s transform. From this formula, we can write

Using some results in Fourier analysis, we can prove that probability function satisfies (if the random variable has a Dirac mass in x)

which can also be written,

And a similar relationship can be obtained if the distribution is absolutely continuous at point ,

Actually, since we work with real-valued random variables, the complex area was just a detour, and we can prove that actually,

It is then possible to get the cumulative distribution function using Gil-Peleaz’s inversion formula, obtained in 1951,

Nice isn’t it. It means, anyone working on financial markets know those formulas, used to price options (see Carr & Madan (1999) for instance). And the good thing is that any mathematical or statistical software can be used to compute those formulas.

  • Characteristic function and actuarial science

Now, what is the interest of all that in actuarial science ? Characteristic functions are interesting when we deal with sums of independent random variables, since the characteristic function of the sum is simple the product of the characteristic functions. They are also interesting when dealing with compound sums1. Consider the problem of computing the 99.5% quantile of the compound sum of Gamma random variable, i.e.

https://latex.codecogs.com/gif.latex?%20S=\sum_{n=1}^N%20X_i

where https://latex.codecogs.com/gif.latex?%20X_i\sim\mathcal{G}(\alpha,\beta) are i.i.d. and https://latex.codecogs.com/gif.latex?%20N\sim\mathcal{P}(\lambda). The strategy is to discretize the loss amounts,

> n <- 2^20; 
> p <- diff(pgamma(0:n-.5,alpha,beta))

Then, the code to compute https://latex.codecogs.com/gif.latex?%20\tilde%20f(s)=\mathbb{P}(S\in[s\pm1/2]), we use

> f <- Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n

To compute the 99.5% quantile, we just use

> sum(cumsum(f)<.995)

That’s extremely simple, isn’it. Want me to do it for real ? Consider the following losses amounts

> set.seed(1)
> X <- rexp(200,rate=1/100)
> print(X[1:5])
[1] 75.51818 118.16428 14.57067 13.97953 43.60686

Let us fit a gamma distribution. We can use

> fitdistr(X,"gamma")
      shape         rate    
  1.309020256   0.013090411 
 (0.117430137) (0.001419982)

or

> f <- function(x) log(x)-digamma(x)-log(mean(X))+mean(log(X))
> alpha <- uniroot(f,c(1e-8,1e8))$root
> beta <- alpha/mean(X)
> alpha
[1] 1.308995
> beta
[1] 0.01309016

Whatever, we have the parameters of our  Gamma distribution for individual losses. And assume that the mean of the Poisson counting variable is

> lambda <- 100

Again, it is possible to use monte carlo simulations, if we can easily generate a compound sum. We can use the following generic code: first we need functions to generate the two kinds of variables of interest,

> rN.P <- function(n) rpois(n,lambda)
> rX.G <- function(n) rgamma(n,alpha,beta)

then, we can use (see here for a discussion on possible codes)

> rcpd4 <- function(n,rN=rN.P,rX=rX.G){
+ return(sapply(rN(n), function(x) sum(rX(x))))}

If we generate one million variables, we can get an estimator for the quantile,

> set.seed(1)
> quantile(rcpd4(1e6),.995)
   99.5% 
13651.64

Another idea is to remember a proporty of the Gamma distribution: a sum of independent Gamma distributions is still Gamma (with additional assumptions on the parameters, but here we consider identical Gamma distributions). Thus, it is possible to compute the cumulative distribution function of the compound sum,

> F <- function(x,lambda=100,nmax=1000) {n <- 0:nmax
+ sum(pgamma(x,n*alpha,beta)*dpois(n,lambda))}

(or at least a approximation). If we invert that function, we get our quantile

> uniroot(function(x) F(x)-.995,c(1e-8,1e8))$root
[1] 13654.43

Which is consistent with our monte carlo computation. Now, we can also use fast Fourier transform here,

> n <- 2^20; lambda <- 100
> p <- diff(pgamma(0:n-.5,alpha,beta))
> f <- Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n
> sum(cumsum(f)<.995)
[1] 13654

Now, if it is simple, is it efficient ? Let us compare for instance computation time to get those three outputs,

> system.time(quantile(rcpd4(1e5),.995))
       user      system     elapsed 
      2.453       0.106       2.611 
> system.time(uniroot(function(x) F(x)-.995,c(1e-8,1e8))$root)
       user      system     elapsed
      0.041       0.012       0.361 
> system.time(sum(cumsum(Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n)<.995))
       user      system     elapsed
      0.527       0.020       0.560

Computations here are comparable with the (numerical) inversion of the cumulative distribution function. Except that here, we were lucky: if the distribution is not Gamma but log normal, the second algorithm cannot be used.

1. This numerical example is taken from the first chapter of Computational Actuarial Science with R, to appear in a few months.

Halloween and candies (a ballot problem)

This year, for Halloween, a post on candies (I promise, next year I will write another post on zombies). But I don’t want to focus on the kids problems (last year, we tried to minimize their walking distance to collect as much candies as possible, with part 1 and part 2), I want to discuss my own problems. Because usually, the kids wear their costumes, and they go in the streets, they knock on the doors, while I stay at home. So I’m the one, with a bag full of candies, waiting for kids to knock on our door, and then I give them some candies (if they wear a costume). Consider the following problem. Assume that we start with  red candies, and  black ones, with . The thing is that no one like those black candies. What could be the probability that for the   kids that will get candies after knocking at my door (with  for convenience, but we will also consider the more general case where I have to many candies, , later on), the probability to get a red candy is always larger than the probability to have a black candy ? This is somehow related to the popular ballot problem, proposed (and solved) by Whitworth in 1878, but he wrote it only in the fourth edition of Choice and Chance, in 1886 (this is what the legend told us). In 1887, Joseph Bertrand proposed a similar problem, and Désiré André introduced the reflection principle to solve it. The problem is simple : consider an election between two candidates, A (who receives  votes) and B (who receives  votes). A wins the election (). If the ballots are cast one at a time, what is the probability that A will lead throughout the voting? For those who don’t remember the conclusion, the probability is here quite simple,

Observe that some geometry proofs were given, later on, by Aebly or Mirimanoff, both in 1923, as well as Howard Grossman in the 1950’s (see the discussion on http://academiclogbook.blogspot.ca/…).  Actually, http://futilitycloset.com// produced the following geometric proof (with no clear reference),

ballot box lattice diagram

We start at O, where no votes have been cast. Each vote for A moves us one point east and each vote for B moves us one point north until we arrive at E, the final count, (mn). If A is to lead throughout the contest, then our path must steer consistently east of the diagonal line OD, which represents a tie score. Any path that starts by going north, through (0,1), must cut OD on its way to E.

If any path does touch OD, let it be at C. The group of such paths can be paired off as p and q, reflections of each other in the line OD that meet at C and continue on a common track to E.

This means that the total number of paths that touch OD is twice the number of paths p that start their journey to E by going north. Now, the first segment of any path might be up to m units east or up to units north, so the proportion of paths that start by going north is n/(m + n), and twice this number is 2n/(m + n). The complementary probability — the probability of a path not touching OD — is (m –n)/(m + n).

But let’s try to solve our problem. Let  and  denote the number of black and red candies, respectively after the th kid git his (or her) candy. Yes, one at a time. Here,  and . What we want is

Using this formulation, we recognize the ballot problem. Almost. Actually, in the original ballot problem (see Bertrand (1887)), we have to compute the probability that one candidate remains strictly ahead the other one throughout the count. With a strict condition, we get the well-known probability (given previously)

Here, ties are allowed, and we can prove (easily) that

(again, there is some nice geometric interpenetration of that result). It is also possible to get numerically that value using the following function, which will generate a trajectory, and return some indicators (with or without ties)

> red_black=function(sd){
+ set.seed(sd)
+ vectcandy=sample(c(rep("R",r),rep("B",b)))
+ v1=rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B"))
+ v2=cumsum(rev(vectcandy)=="R")<= cumsum(rev(vectcandy)=="B") 
+ return(list(evol=cbind(rev(cumsum(rev(vectcandy)=="R")),
+ rev(cumsum(rev(vectcandy)=="B")),v1),list=vectcandy,test=(sum(v1)==0),
+ ballot=(sum(v2)==0),when=min(which(v1==1))))}

(here I compute the ballot-type index, where ties are not allowed, and the candy-type index). If we generate 100,000 scenarios, starting with 50 red and 25 black candies, we get

> r=50 
> b=25 
> M=sapply(1:100000,red_black) ­­
> mean(unlist(M[3,]))
[1] 0.50967

which can be compared with the theoretical value

> (r+1-b)/(r+1)
[1] 0.5098039

We can also get the distribution of the first time we have more black candies than red ones left (given that this event occur)

> Z=unlist(M[5,])
> Z=Z[Z<Inf]
> hist(Z,breaks=seq(0,80),probability=TRUE,col="light blue", 
+ border=NA,xlab="",main="")

There might be some analytically formula that can be derived, but I have to confess that I am becoming extremely lazy,

Assume now that this year, kids do not show up at my door (for some reason). Assume that  kids show up. We can see how the probability

will change, with ,

> r=50
> b=25
> impact_n = function(n){
+ red_black=function(sd,nb=n){
+ set.seed(sd)
+ vectcandy=sample(c(rep("R",r),rep("B",b)))
+ v=(rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")))[1:nb]
+ return(list(list=vectcandy,test=(sum(v)==0),when=min(which(v==1))))}
+ M=sapply(1:10000,red_black)
+ return(mean(unlist(M[2,])))}

Yes, not only I am too lazy to derive analytic formulas, I am so lazy that I do not try to optimize my code. Here, the evolution of the probability, as a function of  is

> V=Vectorize(impact_n)(25:75)
> plot(25:75,V)

Fun isn’t it? But now, I have to conclude my post, to work a little bit on my make-up : I have learnt so many thinks at the Montreal Zombie Walk a few days ago that kids willing to knock at my door will be scared to death. I guess I will keep all the candies for me this year !