Category Archives: Maths

Le retour des abeilles

Suite à mon précédant billet, “Maudites Abeilles“, j’ai eu des commentaires sur le blog, mais aussi sur Twitter pour me dire que malgré tout, il y avait bien une chance sur deux pour avoir un chemin bleu permettant de connecter les deux régions bleues, celle au nord et celle au sud.

En effet, dans le problème de base, les régions sont colorées, et le problème n’est pas l’existence d’un chemin, mais l’existence d’un chemin de la bonne couleur. Il faut juste rajouter une petite ligne pour spécifier la couleur du chemin

> simu2=function(d){
+ C=sample(c(rep(1,d^2/2),rep(2,d^2/2)),size=d^2)
+ B=NULL
+ for(i in 1:(d^2)){
+ x=rep(i,6)
+ y=c(i-d,i-d+1,i-1,i+1,i+d-1,i+d)
+ if(i%%d==1) y=c(i-d,i-d+1,NA,i+1,NA,i+d)
+ if(i%%d==0) y=c(i-d,NA,i-1,NA,i+d-1,i+d)
+ D=data.frame(x,y)
+ D=D[(D$y>=1)&(D$y<=d^2),] + B=rbind(B,D[which((C[D$y]==C[D$x])&(C[D$x]==1)),]) + } + B=B[(B[,2]>=1)&(B[,2]<=d^2),] + G=as.vector(t(B)) + G_iter=connect(make_graph(G), d^2) + connectpt=function(x) as.numeric(adjacent_vertices(G_iter,x)[[1]]) + sum(unlist(Vectorize(connectpt)(1:d))%in%(d^2+1-(1:d)))>0}
> appr=function(d,nsimu) mean(Vectorize(simu2)(rep(d,nsimu)))
> appr(4,10000)
[1] 0.4993

et en effet, il y a une chance sur deux de trouver un chemin bleu permettant de connecter les deux régions bleues. Mais comme je le disais hier, j’ai du mal avec l’argument de symétrie évoqué dans la “correction”.

Maudites Abeilles

Pendant les “vacances”, j’étais avec les enfants (le plus grand et la plus petite) au palais de la découverte, à Paris, à l’exposition temporaire “faîtes vos jeux“. Comme toujours, je suis surpris quand je compare les musées de sciences, et les musées d’art (j’avais fait un billet voilà quelques années en comparant le MoMaths et le Gugenheim, à New York City). Autant ma grande fille est autant intéressée par les tableaux que par les sciences, autant les deux autres préfèrent mille fois les seconds au premier. On a d’ailleurs passé plus de 4 heures au musée (dont une à faire toutes les animations de “faîtes vos jeux“), il a fallu y aller parce que j’avais une réunion de boulot (oui, les parisiens n’étaient pas en vacances, et j’avais une réunion impossible à déplacer). D’où les guillemets à “vacances” (le second jour de notre excursion, c’était une intervention de deux heures en fin d’après midi).

Mais je m’égare… A un des stands de l’animation, il avait le défi des abeilles,

Sur notre “damier” 8×8 (ici), on a 32 cases bleues et 32 rouges (ou oranges), et on cherche la probabilité qu’un chemin permette de passer d’un côté à l’autre. Sur le moment, le problème m’a pas mal pris la tête, et je me disais que la réponse devait être compliquée… en tout cas, que ça n’était pas aussi simple que ça…

En ligne, on peut trouver le guide de l’enseignant qui apporte des compléments aux divers ateliers, et qui propose, en l’occurrence, une solution

La réponse serait-ici 1/2. Mais la réponse ne me satisfait pas, en particulier la phrase “son absence à l’inverse impose un chemin orange” qui suppose qu’il existe forcément un chemin, soit rouge, soit bleu, pour passer d’un côté à l’autre, ce qui ne me semble pas trivial ! Quant

Par exemple sur le damier 4×4 ci-dessous,

il n’existe aucun chemin rouge allant  du bas vers le haut, et aucun chemin bleu non plus ! Et si on regarde, il existe (au moins) un chemin rouge allant de la gauche vers la droite, et un chemin bleu. Bref, la réponse 1/2 me semble douteuse….  Pour ceux qui veulent de vrai hasard (mais qui a autant de chance d’apparaître que celui au dessus), on a le graph ci-dessous

Mais maintenant essayons de la calculer cette probabilité. Par simulations au moins. Regardons un damier,

on peut le visualiser sous forme d’une matrice

ou en décalant un peu d’une matrice carrée 4×4 (ici)

Si on regarde ce que sont les voisins (afin de faire un chemin), sur le schéma initial, on cherchait les configuration suivantes

ce qui, dans notre matrice, correspond aux éléments ci-dessous,

Bref, on va pouvoir connecter des points d’un graph, et regarder si les 4 premiers éléments du graphs sont connectés aux 4 derniers

(le but est désormais de savoir si le nord et le sud connectent, peu importe la couleur). On peut commencer par créer le graph

d=4
set.seed(1)
C=sample(c(rep(1,d^2/2),rep(2,d^2/2)),size=d^2)
M=matrix(C,d,d)
> M
[,1] [,2] [,3] [,4]
[1,] 1 1 2 2
[2,] 1 2 1 1
[3,] 2 2 1 1
[4,] 2 2 1 2

ce qui correspond bien à notre matrice ci-dessous

On va alors créer le graph, c’est à dire chercher les noeuds qui communiquent

library(igraph)
B=NULL
for(i in 1:(d^2)){
x=rep(i,6)
y=c(i-d,i-d+1,i-1,i+1,i+d-1,i+d)
if(i%%d==1) y=c(i-d,i-d+1,NA,i+1,NA,i+d)
if(i%%d==0) y=c(i-d,NA,i-1,NA,i+d-1,i+d)
D=data.frame(x,y)
D=D[(D$y>=1)&(D$y<=d^2),] B=rbind(B,D[which(C[D$y]==C[D$x]),]) } B=B[(B[,2]>=1)&(B[,2]<=d^2),]
G=as.vector(t(B))
m=make_graph(G)

On peut ainsi voir ici que 1 communique avec 2, et 5

> B[B$x==1,]
x y
4 1 2
6 1 5

On peut d’ailleurs visualiser le graph, et voir qu’il est constitué de 4 graphs disjoints,

> plot(m)

(le coin 16 étant isolé, sans connection avec le reste – étrangement – il n’apparait pas sur le graph)

Pour trouver tous les points qui communiquent entre eux, on peut itérer un peu (voire beaucoup, peu importe)

G_iter=connect(make_graph(G), d^2)
connectpt=function(x) as.numeric(adjacent_vertices(G_iter,x)[[1]])

Ici, le nord et le sud communiquent s’il y a au moins un chemin (peu importe la longueur)

sum(unlist(Vectorize(connectpt)(1:d))%in%(d^2+1-(1:d)))>0

On peut alors une petite fonction, qui va générer une répartition des couleurs sur un damier carré (de longueur paire)

simu=function(d){
C=sample(c(rep(1,d^2/2),rep(2,d^2/2)),size=d^2)
B=NULL
for(i in 1:(d^2)){
x=rep(i,6)
y=c(i-d,i-d+1,i-1,i+1,i+d-1,i+d)
if(i%%d==1) y=c(i-d,i-d+1,NA,i+1,NA,i+d)
if(i%%d==0) y=c(i-d,NA,i-1,NA,i+d-1,i+d)
D=data.frame(x,y)
D=D[(D$y>=1)&(D$y<=d^2),] B=rbind(B,D[which(C[D$y]==C[D$x]),]) } B=B[(B[,2]>=1)&(B[,2]<=d^2),] G=as.vector(t(B)) G_iter=connect(make_graph(G), d^2) connectpt=function(x) as.numeric(adjacent_vertices(G_iter,x)[[1]]) sum(unlist(Vectorize(connectpt)(1:d))%in%(d^2+1-(1:d)))>0}

On peut alors construire la fonction suivante, pour faire un peu de monte carlo

appr=function(d,nsimu) mean(Vectorize(simu)(rep(d,nsimu)))

Par exemple sur notre damier 4×4, on a environ 3 chances sur 4 d’avoir un chemin (peu importe sa couleur) qui permette de rejoindre le sud depuis le nord

> appr(4,10000)
[1] 0.7446

alors que sur le 8×8 du palais de la découverte, on a

> appr(8,10000)
[1] 0.7805

78% de chances (un peu plus). Ce qui semble confirmer qu’il n’y a pas forcément un chemin permettant d’aller du nord au sud.

On peut d’ailleurs faire un petit dessin pour noter que quand la dimension du damier augmente (je n’ai pas tenté ici le cas impair – avec une case de plus d’une certaine couleur), la probabilité de trouver au moins un chemin semblerait augmenter (même si ça reste à prouver, et en particulier, je serais curieux de voire vers quoi ça converge…)

 

Être – en même temps – moyen et extraordinaire

Juste avant 1850, Adolphe Quételet publiait Sur l’appréciation des documents statistiques, et en particulier sur l’application des moyennes, et introdusait alors le concept d’homme moyen (l’ouvrage est cité dans La théorie des moyennes et son emploi dans les sciences d’observation, de Joseph Lottin, paru en 1909 – je n’ai pas trouvé en ligne cet ouvrage de Quételet). Ce concept sera détaillé encore davantage dans Anthropométrie, ou mesure des différentes facultés de l’homme, concept qui sera abondament illustré

Mais si la moyenne est un concept fondamental en statistique, il convient de garder en mémoire que la moyenne n’existe peut être pas, au sein de la population. C’est ce que rappelle la fameuse ‘blague’ “un mec qui à la tête dans le four et les pieds dans le frigo, en moyenne, il est bien“. La moyenne peut être quelquechose d’extraordinaire (au sens premier du terme).

Continue reading Être – en même temps – moyen et extraordinaire

Trouver le millionième chiffre

Dimanche, avec mon fils, on a tenté de résoudre un petit problème posé par Ian Stewart (qu’est-ce qu’il faut pas inventer pour éviter de voir mon fils sur sa tablette toute la matinée). Le problème est très simple, dans sa formulation : considérons la suite suivante

> paste(1:26,collapse="")
[1] "1234567891011121314151617181920212223242526"

Quel serait le millionième chiffre ? Regardons la réponse,

> L=paste(1:999999,collapse="")
> f=function(n) substr(L,n,n)
> f(1000000)
[1] "1"

J’espérais trouver une logique (genre “c’est toujours un 1“), en vain. Le point de départ était “si on va jusqu’à écrire 99“, combien de chiffres faut-il ? Facile

> 9+90*2
[1] 189

et maintenant “si on va jusqu’à 999” ?

> 9+90*2+900*3
[1] 2889

Etc. On peut utiliser la petite fonction

> g=function(k) sum(9*10^(1:k-1)*(1:k))

Pour le millionième chiffre, quel était la taille du nombre ?

> g(5)
[1] 488889
> g(6)
[1] 5888889

C’était manifestement un nombre à 6 chiffres: quand on aura écrit le dernier nombre à 5 chiffres – 99999 – il y aura eu, en tout, 488889 chiffres. Pour savoir quel chiffre on aura atteint au millionième chiffre, notons que

> (1000000-g(5))%/%6
[1] 85185

autrement dit, on cherche le 85185ème nombre à 6 chiffres, soit

> 100000+(1000000-g(5))%/%6
[1] 185185

Et pour trouver quel chiffre se trouve précisément à la position 1000000, rien de plus simple

> (1000000-g(5))%%6
[1] 1

autrement dit, c’est le premier chiffre du nombre 185185, qui est un 1. C’est le calcul qu’on avait commencé avec mon fils… mais je me demandais s’il n’y avait pas un truc… Si on regarde pour le millième chiffre, ou le dix millième chiffre, ou le cent millième chiffre, il ne semble pas y avoir de logique

> f(1000)
[1] "3"
> f(10000)
[1] "7"
> f(100000)
[1] "2"
> f(1000000)
[1] "1"

Pourtant,  il y a une forte régularité, quand on regarde la répartition des 10 chiffres (de 0 à 9), dans la séquence

> library(RColorBrewer)
> COL= brewer.pal(10, "Set3")
> k=100
> V=rep(NA,k*k)
> plot(0:k,0:k,col="white",xlab="",ylab="")
> for(i in 1:k){
+   for(j in 1:k){
+     n=(i-1)*k+(j-1)
+     x=as.numeric(f(n))
+     rect(i-1,j-1,i,j,col=COL[x],border=NA)
+     V[n]=x
+   }
+ }

Mais rien de probant… Plus étonnant, si on regarde la distribution des chiffres, sur les 100000 premiers,

> Z=Vectorize(function(x) as.numeric(f(x)))(1:1e5)
> table(Z)[1:8]
    0     1     2     3     4     5     6     7 
 8642 19753 11111  8642  8642  8642  8642  8642

Le 1 ressort très majoritairement… le 2 se distingue aussi, mais ensuite, c’est étonnamment constant

> barplot(table(Z))

 

On peut d’ailleurs noter que la différence entre le nombre de 1 et le nombre de 2 est très exactement le nombre de 0, de 3, de 4, de 5, …

> table(Z)[3]-table(Z)[2]
    2 
-8642

 

étonnant, non ?

What happens if we forget a trivial assumption ?

Last week, @dmonniaux published an interesting post entitled l’erreur n’a rien d’original  on  his blog. He was asking the following question : let https://latex.codecogs.com/gif.latex?a, https://latex.codecogs.com/gif.latex?b and https://latex.codecogs.com/gif.latex?c denote three real-valued coefficients, under which assumption on those three coefficients does https://latex.codecogs.com/gif.latex?ax^2+bx+c has a real-valued root ?

Everyone aswered https://latex.codecogs.com/gif.latex?b^2-4ac\geq%200, but no one mentioned that it is necessary to have a proper quadratic equation, first. For instance, if both https://latex.codecogs.com/gif.latex?a and https://latex.codecogs.com/gif.latex?b are null, there are no roots.

Continue reading What happens if we forget a trivial assumption ?

La formation française, élitiste ?

En discutant avec @tomroud la semaine passée, j’ai découvert l’article Notre élite est toujours excellente, c’est le reste de la classe qui ne suit pas sur le site du Monde. C’est effectivement une idée reçue : si la formation est globalement aussi faible en France, c’est parce que la France serait un pays élitiste, et la formation de ses élites, en revanche, est une réussite. C’est – en gros – ce que raconte cet article. On pourra relire à ce sujet l’échange

Comment expliquer les résultats faibles en mathématiques, malgré la tradition d’excellence de la France dans ce domaine ?

Notre élite est toujours excellente, mais c’est le reste de la classe qui ne suit pas. Nos mauvaises notes à PISA ne nous empêcheront pas d’avoir des médailles Fields, mais elles nous empêchent d’avoir 100 % de citoyens efficaces dans leur travail et capables de se débrouiller dans un monde très mathématiques.

(dans la discussion, avec également, @Lazarrean@rdoume avait signalé un autre article dans Le Monde, la France championne des inégalités scolaires). Pourtant, ce n’est pas tout à fait ce que raconte le graphique, mentionné dans Surprising Test Results For Some Of The World’s Richest Students

où le “Ninth decile” (quantile à 90%) n’est pas si élevé que ça, comparé aux autres pays. Maintenant, pour raconter l’histoire complète, la discussion avec @tomroud est née de la lecture de La crise ne rend pas les riches plus riches et les pauvres plus pauvres, sauf en France,

En effet, @tomroud me disait qu’il avait vu le même genre de conclusion dans la comparaison des études PISA : en France, les bons sont meilleurs avec le temps, même si la moyenne, tout comme les niveaux des déciles les plus faibles, baisse. Dans une vision élitiste de l’éducation, ça ne serait pas gênant.

N’ayant rien trouvé me permettant de mieux comprendre , j’ai voulu aller regarder un peu plus en détails ce que disent vraiment les enquêtes PISA. Car elles sont disponibles en ligne (20122009200620032000). Bon, pour les lire plus rapidement, je vais les récupérer sur d’autres sites, pour avoir d’autres formats, plus simples à lire, en R,

> load(url("https://github.com/pbiecek/PISA2000lite/blob/master/data/math2000.rda"))
> load(url("http://beta.icm.edu.pl/PISAcontest/data/student2012.rda"))

Commençons par regarder un peu la distribution des notes, en mathématiques (je vais me restreindre aux tests de maths, qui sont souvent utilises pour dire que la France est effectivement un pays libelliste : les maths ont longtemps servi de critère, me semble-t-il).

> py="France"
> xf=sb[(sb$CNT==py)&(sb$ST04Q01=="Female")==TRUE,"PV1MATH"]
> xm=sb[(sb$CNT==py)&(sb$ST04Q01=="Male")==TRUE,"PV1MATH"]

Sur le graphique ci-dessous, on a la distribution des notes des filles (en haut à gauche) et des garçons (en bas à gauche). A droite, on a la fonction de survie : à une note donnée, on a en ordonnées la proportions d’élèves qui dépassent cette note. Par exemple 18,72% des garçons dépassent 600, et 12,42% des filles

> mean(xf>600)
[1] 0.1242105
> mean(xm>600)
[1] 0.1872207

Enfin, en bas à gauche, j’ai rajouté la fonction de moyenne en excédant (ou pour être plus précis, la Tail-Value-at-Risk), qui est la note moyenne que l’on peut espérer avoir si dans un certain décile. Par exemple, si on s’intéresse aux élèves dans le décile supérieur (dépassant le quantile à 90%), la note moyenne des garçons est 672 et celle des filles 649,

> mean(xm[xm>quantile(xm,.9)])
[1] 672.7608
> mean(xf[xf>quantile(xf,.9)])
[1] 649.281

Cette mesure est assez classique en tant que downside risk measure et permet de focaliser justement sur les queues de distribution (ce qui nous intéresse ici car on veut étudier les élites).

Si on zoom dans la partie à droite (les meilleures notes), on obtient

Mais il est difficile de faire des comparaisons avec d’autres pays sur ces graphiques. Par exemple, pour l’Australie, on obtient

 

Pour comparer les pays entre eux, on peut regarder la différence entre la note moyenne, et le niveau du quantile à 99% (on va chercher la vraie élite, les 1% des meilleurs… car c’est manifestement ce qui été évoque dans l’article, et qui intéresse du monde sur les inégalités de richesse)

L’élite en France est pas incroyablement meilleure que le cœur de la distribution, on en conviendra !

Si on compare un pays avec la France, on peut regarder la différence entre le quantile de niveau u (e.g. 50%, ou 90% ou 99%) entre ce pays et la France. En niveau absolu, ou en pourcentage.

> qpays("France",c(.5,.99))
      Female     Male
50% 496.8141 507.5635
99% 691.0734 717.2885

Sur les graphique ci-dessous, on a le différence absolue à gauche et relative à droite, avec l’Allemagne en haut, et l’Australie en bas,

On peut aussi regarder la Thaïlande et le Mexique (on est toujours en 2012)

Regardons un peu la Thaïlande par exemple. Les élèves (scolarises) médians sont 50 points en dessous des élèves Français. Mais pour le top 1%, les notes sont comparables… Donc manifestement, la formation des élites semble fonctionner dans d’autres pays que la France. En Australie, les inégalités scolaires semble encore plus fortes : les élèves médians sont globalement moins bon en Australie qu’en France (disons 15 points en dessous), mais l’élite est au dessus (de presque 50 points quand on regarde le top 1%).

Cela dit, comme j’avais dit qu’il pouvait aussi être intéressant de regarder les fonctions de moyennes en excédant (plus que les quantiles), on peut comparer ces fonctions,

Ici, les graphiques sont présentés en variation relative : au Mexique, la note moyenne d’une fille dans le top 20% (au dessus du quantile à 80%) est 17.5% plus faible qu’en France. On retrouve ici qu’en Australie par exemple, l’élite a – en moyenne – de meilleures notes qu’en France. Et en Thaïlande, plus on regarde les meilleurs élèves, plus les résultats sont comparables avec ceux des élèves français.

Qu’en est-il maintenant de l’évolution dans le temps ? Si on compare dans le temps, entre 2000 et 2012.

> sb2012=student2012[,c("PV1MATH","CNT","ST04Q01")]
> sb2000=math2000[,c("PV1MATH","CNT","ST03Q01")]
> sb2000=sb2000[!is.na(sb2000$PV1MATH),]
> sb2000$cnt=tolower(sb2000$CNT)
> sb2012$cnt=tolower(sb2012$CNT)
> ext2000=function(py){
+ x=sb2000[sb2000$cnt==py,"PV1MATH"]
+ v=c(mean(x),quantile(x,c(.1,.9)))
+ return(v)
+ }
> ext2012=function(py){
+ x=sb2012[sb2012$cnt==py,"PV1MATH"]
+ v=c(mean(x),quantile(x,c(.1,.9)))
+ return(v)
+ }
> listepays2000=unique(sb2000$cnt)
> listepays2012=unique(sb2012$cnt)
> R2000=sapply(listepays2000,ext2000)
> R2012=sapply(listepays2012,ext2012)
> D2000=data.frame(pays=listepays2000,
+ m2000=R2000[1,],
+ q12000=R2000[2,],
+ q92000=R2000[3,])
> D2012=data.frame(pays=listepays2012,
+ m2012=R2012[1,],
+ q12012=R2012[2,],
+ q92012=R2012[3,])
> D=merge(D2000,D2012)

Si on fait un dessin proche de celui sur les revenus (avec les évolutions des moyennes en bâtons, des 10% les plus bas en rond, et des 10% les plus hauts en rond), on obtient

Si effectivement le niveau moyen des étudiants baisse en France (comme dans plusieurs autres pays), on notera que l’idée que l’élite française s’en sort bien est un cliche sans réel fondement.

L’échelle logarithmique expliquée aux juristes (et aux travailleurs)

Comme le rappelle https://sciencenews.org/blog/, il y a 400 ans, John Napier (que les français ont décidé d’appeler Neper) a ‘inventé’ le logarithme. Il publia son travail dans Mirifici logarithmorum canonis descriptio, et pour lui rendre hommage, le logarithme que l’on utilise le plus souvent est le logarithme népérien (qui est l’inverse de la fonction exponentielle).

Comme souvent, quand je parle de maths sur mon blog, on me reproche de ne m’adresser qu’à une élite (le commun des mortels dira des nerds) en employant un jargon incompréhensible. Pourtant, les logarithme sont présents partout !

Continue reading L’échelle logarithmique expliquée aux juristes (et aux travailleurs)

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

Le Père Noël et le pari de Pascal (spoiler alert)

Les enfants ont beau grandir, Noël reste un moment magique. Bien sûr, le plus grand nous dit depuis des années qu’il n’y croit plus, et sa sœur le suit. Mais quand même… on a toujours l’impression que juste au cas où, ça pourrait valoir le coup d’y croire. Comme le Dieu de Pascal. Ouais, c’est bien de faire le grand en disant qu’on n’y crois plus, mais s’il existait, même avec une très faible probabilité, ça ne vaudrait pas la peine d’y croire encore un peu ? Pour ceux qui ne connaîtra pas le pari de Pascal, il existe une version courte, et illustrée, avec Calvin et Hobbes,

Pour la version originale, signée Blaise Pascal, on pourra relire le paragraphe 233 des Pensées,

Dans ce jeu, la stratégie minimax pousserait à croire en Dieu, pour éviter d’aller en enfer (correspondant à une perte infinie, que rien ne peut compenser). Il y a plusieurs points commun entre les deux situations. Par exemple, je trouve que les parents reproduisent bien le comportement de l’église, en faisant planer le doute et en utilisant des rituels soi-disant millénaires. Bon, la grande différence, c’est qu’on n’a qu’une vie (ou qu’une mort), mais il me semble que Noël revient tous les ans. En plus, si Dieu n’existe pas, il n’y a rien à la place… alors que si le père Noël n’existe pas…

(pour d’autres réflexions sur Noël, http://buzzfeed.com/ariannarebolini/…)

Random points on some hemisphere

In my previous post, I tried to answer the following question Consider https://latex.codecogs.com/gif.latex?n points uniformly distributed on a sphere. What is the probability that the https://latex.codecogs.com/gif.latex?n points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?

If I have been able to use Monte Carlo simulations in dimension 2 (on a circle, not on a sphere), I could not get it in dimension 3. Hopefully, my colleague Simon gave me a nice solution (much more efficient than my previous code). Unfortunately, it was in Maple,

clear all
taillesim=5000;
for N=3:10
nb=0; 
for k=1:taillesim
X=randn(N,1);
Y=randn(N,1);
Z=randn(N,1);
pts=[X Y Z];
for i=1:N  
  pts(i,:)=pts(i,:)/norm(pts(i,:),2);
end
tol=1e-07;
A=-pts;
f=[0;0;0];
b=zeros(N,1)-tol;
%options=optimset('Display','off');
[x,val,exitflag]=linprog(f,A,b);
if(exitflag==1)
  nb=nb+1;    
end
end
probapprox(N-2)=nb/taillesim;
end
ns=[3:10];
for i=3:10
probhs(i-2)=probh(3,i);
end
scatter(ns,probhs)
hold on
plot(ns,probapprox);

The idea is very clever (I did not know how to code it, but as we will see, it is actually simple – or say not too difficult – to code actually). The idea is based on the idea that all the points lie on the same hemisphere if there is some  such that  for all .

 

That is not a big improvement, compared with the previous post, because we have to find such a vector. The idea suggested by Simon is to use some constraint optimization algorithm. We can try to optimize (maximize, or minimize, actually, we do not care) something like

If the set of constraint is empty, then there is no solution, while we can find one (maybe more, but we do not care) if there is such a vector . I wanted to add some constraint on , like  (for the Euclidean norm), but that would be a non-linear constraint. So here, I chose to assume that  was such that . With R, the optimization routine would be like

library(lpSolve)
A=rbind(pts,c(1,1,1))
C=c(1,2,3)
B=c(rep(1e-10,N),1)
slp=lp("min",C,A,c(rep(">=",N),"="),B)

So that I can solve here

and then, we simply have to check for the value of

slp$status

0 means that we found a vector  . But when you run the code, it does not work. Because actually, what is solved here is not exactly the program above, but the one below,

(yes, it is mentioned in the help of the optimization function that only positive values are considered). So, a simple strategy, since we focus on an orthant here, is to consider all possible orthants. For instance, using

test=FALSE
x1=0:(2^3-1)%/%(2^(3-1))
x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2))
x3=((0:(2^3-1))-(2^2)*x1-2*x2)
for(u in 1:(2^3)){
pts2=cbind(pts[,1]*(-1)^x1[u],
           pts[,2]*(-1)^x2[u],
           pts[,3]*(-1)^x3[u])
A=rbind(pts2,c(1,1,1))
C=c(1,2,3)
B=c(rep(1e-10,N),1)
slp=lp("min",C,A,c(rep(">=",N),"="),B)
if(slp$status==0) test=TRUE}
if(test==TRUE) nb0=nb0+1

So here, the code would be something like

taillesim=5000
probaap=rep(NA,10)
probath=rep(NA,10)
p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))
for(N in 3:10){
nb0=0
for(k in 1:taillesim){
MZ=matrix(rnorm(n*3),n,3)
d=apply(MZ,1,function(z) sqrt(sum(z^2)))
X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d; 
pts=cbind(X,Y,Z)
test=FALSE
x1=0:(2^3-1)%/%(2^(3-1))
x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2))
x3=((0:(2^3-1))-(2^2)*x1-2*x2)
for(u in 1:(2^3)){
pts2=cbind(pts[,1]*(-1)^x1[u],
           pts[,2]*(-1)^x2[u],
           pts[,3]*(-1)^x3[u])
A=rbind(pts2,c(1,1,1))
C=c(1,2,3)
B=c(rep(1e-10,N),1)
slp=lp("min",C,A,c(rep(">=",N),"="),B)
if(slp$status==0) test=TRUE
}
if(test==TRUE) nb0=nb0+1  
}
probaap[N]=nb0/taillesim
probath[N]=p(3,N)
}

and this time, the probability obtained using Monte Carlo simulation is extremely close to the theoretical value.

And the algorithm is extremely fast! So using optimization routines to see if sets defined by linear constraints are empty – or not – is truly a great idea!

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!

Random points on the Earth

The problem with puzzles is that you keep it in your head for days, until you find an answer. Or at least some ideas about a possible answer. This is what happened to me a few weeks ago, when a colleague of mine asked me the following question : Consider https://latex.codecogs.com/gif.latex?n points uniformly distributed on a sphere. What is the probability that the https://latex.codecogs.com/gif.latex?n points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?

Analogously, what is the probability to see the https://latex.codecogs.com/gif.latex?n points on the Earth, at the same time, from somewhere in the galaxy ? (even extremely far away, so we can see a complete hemisphere) I wanted to use Monte Carlo simulations to estimate that probability, for some https://latex.codecogs.com/gif.latex?n. But it was difficult. I mean, given https://latex.codecogs.com/gif.latex?n points on the sphere, in  can you easily determine if they lie on a common hemisphere, or not ? I did try with distance, or angle, but I could not find a simple answer. So I tried a technique I did learn a few years ago : if you cannot do something in dimension 3, try first in dimension 2.

Again, I could not find a simple answer. But there is a simple technique.

  1. Draw https://latex.codecogs.com/gif.latex?n points on the unit sphere, which simply means generate https://latex.codecogs.com/gif.latex?n random variables https://latex.codecogs.com/gif.latex?%20U_1,\ldots%20,U_n\sim\mathcal{U}([0,2\pi])
  2. Try to find https://latex.codecogs.com/gif.latex?\theta\in[0,2\pi] such that, after a rotation (with angle https://latex.codecogs.com/gif.latex?\theta) all the points lie in the upper part (the North hemisphere, for instance)

So the question here is simply : is

https://latex.codecogs.com/gif.latex?\max_{\theta\in[0,2\pi]}\left\{\sum_{i=1}^n%20\boldsymbol{1}(\sin(U_i+\theta)\geq%200)\right\}

equal to https://latex.codecogs.com/gif.latex?n ? Of course, from a computational point of view, it is slightly more complex, since this function is not differentiable.

n=5
Theta=runif(n)*2*pi
top=Vectorize(function(theta) sum(sin(Theta+theta)>=0))

So a simple strategy can be to compute those values on a finite grid, and to check if, for some https://latex.codecogs.com/gif.latex?\theta, all the points lie in the upper part

max(top(seq(0,2*pi,length=6001)))==n

Hence, with the following sample, all the points cannot be on a common hemisphere

set.seed(2)
Theta=runif(5)*2*pi

while, for another sample of points, it was possible

set.seed(7)
Theta=runif(5)*2*pi

With this simple code, we get get the probability, but only in dimension 2 (so far)

SIM=Vectorize(function(n) simul(n,1000))
plot(3:10,SIM(3:10))

In dimension 3, it is still possible to use also a polar representation. Things are easier to generate, but also it is simple to consider rotations. And again, a simple algorithm can be derived,

But there is a simple technique.

  1. Draw https://latex.codecogs.com/gif.latex?n points on the unit sphere, which simply means generate https://latex.codecogs.com/gif.latex?2n random variables https://latex.codecogs.com/gif.latex?%20U_1,\ldots%20,U_n\sim\mathcal{U}([0,2\pi]) and https://latex.codecogs.com/gif.latex?V_1,\ldots,V_n\sim\mathcal{U}([-\pi/2,\pi/2])
  2. Try to find https://latex.codecogs.com/gif.latex?\theta and https://latex.codecogs.com/gif.latex?\varphi such that, after a rotation (with angles https://latex.codecogs.com/gif.latex?\theta and https://latex.codecogs.com/gif.latex?\varphi) all the points lie in some given part (say  https://latex.codecogs.com/gif.latex?\{x\geq0\})

As mentioned by Dominique, using this technique, points are not uniformly distributed on the sphere. Instead we can use

  1. Draw https://latex.codecogs.com/gif.latex?n points on points from a trivaraite Gaussian distribution, and normalize it
  2. Get the polar coordinates, and try to find https://latex.codecogs.com/gif.latex?\theta and https://latex.codecogs.com/gif.latex?\varphi such that, after a “rotation” (with angles https://latex.codecogs.com/gif.latex?\theta and https://latex.codecogs.com/gif.latex?\varphi) all the points lie in some given part (say  https://latex.codecogs.com/gif.latex?\{x\geq0\})

So the question here is simply : is

https://latex.codecogs.com/gif.latex?\max_{(\theta,\varphi)\in[0,2\pi]\times[-\pi/2,\pi/2]}\left\{\sum_{i=1}^n%20\boldsymbol{1}(\sin(U_i+\theta)\cos(V_i+\varphi)\geq%200)%20\right\}

(I am not sure about the set of angles for the rotations, so I tried a larger one, just in case). Again, it would be complex, or more complex than before, because we need here a joint grid. For instance

MZ=matrix(rnorm(n*3),n,3)
d=apply(MZ,1,function(z) sqrt(sum(z^2)))
X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d; 
Theta=acos(Z)
Phi=acos(X/sqrt(X^2+Y^2))*(Y>=0)+(2*pi-acos(X/sqrt(X^2+Y^2)))*(Y<0)
top=function(theta,phi) sum(sin(Theta+theta)*cos(Phi+phi)>=0)
TOP=mapply(top,rep(seq(0,2*pi,length=1001),1001),rep(seq(-pi,pi,length=1001),each=1001))
max(TOP)==n

As we can see with the red curve, there might be some problems here. Because, (as mention in kmath327), this problem was solved in any dimension in Wendel (1962) with the following simple expression : in dimension https://latex.codecogs.com/gif.latex?d, the probability that the https://latex.codecogs.com/gif.latex?n points (uniformly distributed on a sphere) lie on a same hemisphere is exactly

https://latex.codecogs.com/gif.latex?p(d,n)=%20\frac{1}{2^{n-1}}\sum_{i=0}^{d-1}%20\binom{n-1}{i}

p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))

Note that Leonard Savage proved (a few years before) that

https://latex.codecogs.com/gif.latex?p(d,d+1)=1%20-\frac{1}{2^d}

(which can be obtained easily actually). I do not see what’s wrong with my Monte Carlo simulations… and if anyone has a nice Monte Carlo strategy to get that probability, I’d be glad to hear it !

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)

Mathématiques de l’Assurance Non-Vie (2)

« Dans ce contexte d’incertitude, il est particulièrement réconfortant de revenir aux sources, aux fondamentaux, c’est-à-dire aux mathématiques et de rappeler que le risque naît de l’aléa et s’appréhende grâce aux développements les plus avancés du calcul des probabilités. » (Claude Bébéar, dans la préface du tome 1)

Après avoir été épuisé plusieurs semaines (voire plusieurs mois ?), une nouvelle impression du tome 2 de Mathématiques de l’Assurance Non-Vie, coécrit avec Michel Denuit, est de nouveau disponible. J’ai été surpris quand on nous a dit – en début d’année – que le second tome serait bientôt épuisé, et que plus d’un millier d’exemplaires avaient été vendus (ou donnés, ou volés… peu importe). Surpris qu’autant de monde pouvait être intéressé par un gros pavé (les deux tomes doivent faire plus de mille pages, au total) aussi sérieux. Comme le notait Gilles Bénéplanc dans la revue Risques, « s’ils sont riches et rigoureux (le titre présente l’avantage ne pas être trompeur quant au contenu), ces deux volumes d’actuariat se lisent toujours avec intérêt et même souvent avec plaisir ». Car effectivement, on a passé du temps à choisir le titre, et on est vite tombé d’accord sur le fait qu’il fallait insister sur le fait qu’il s’agissait de livre de Mathématiques (avant tout). Pour la petite anecdote, j’étais allé le chercher chez Gibert (juste après la publication du premier tome… après avoir fait un détour par les rayons de bandes dessinées, bien sûr). En vain. Là, un vendeur m’a expliqué que le livre était considéré comme traitant d’assurance, et que l’assurance était rangé dans les rayons de droit. Bref, notre livre vert était perdu au milieu des livres de Dalloz (alors que les livres de mathématiques financières se trouvaient dans la section mathématiques). Mais c’est un détail, tant qu’en le cherchant, on le trouve…

Je dois aussi dire que je suis très fier du livre. Très fier que des étudiants et des praticiens me disent le lire. Très fier d’apprendre que des étudiants au Brésil avaient photocopié (et traduit) l’introduction aux modèles de ruine (que l’on introduit avec l’approche de de Finetti, beaucoup plus intuitive je trouve…). Très fier quand des chercheurs que j’admire me disent avoir trouvé un petit résultat qu’ils cherchaient pour un papier de recherche (je pense a des discussions avec Alain, à Paris 1, sur la transformée d’Esscher). Je trouve aussi dangereux qu’autant de monde le lisent, car il est truffé d’erreurs. Mais après près de 6 mois de relecture complète, presque tous les jours, on fini par ne plus rien voir. Et malgré le travail incroyable de Christian, éditeur de la collection, il reste des coquilles (des coquillettes j’espère). Mais si on devait attendre que les livres soient parfaits avant de les publier, rien ne sortirait. Après presque dix ans, je pense que je réécrirais des paragraphes d’une autre manière, ou que je mettrais plus d’aspects algorithmiques. Mais ce qui est fait est fait, et j’en suis encore assez fier. OK, très fier.

Lorsque l’on a été contacté, on nous a demandé s’il fallait un  nouveau tirage, ou si une seconde édition (revue et corrigée) serait nécessaire. Soyons réalistes. Le livre n’est plus forcément à jour, et comme je l’ai dit, il reste beaucoup de coquilles, voire de fautes (il faudrait que l’on prenne le temps de les recenser). Mais une nouvelle édition ferait probablement doubler la taille du livre, et prendrait beaucoup de temps ! Sur les deux premiers chapitres, Michel a publié un livre complet sur les données de comptage en assurance. Sur les aspects computationnels, je vais éditer un livre à venir. Et s’il y a 10 ans, il y a avait encore peu de livres sur le sujet, plusieurs permettent maintenant de rentrer plus en détails sur certains aspects. Je pense au livre de Mario Wüthrich et Michael Merz sur le provisionnement (plus clair que le chapitre de notre livre, et plus complet). Pareil pour les chapitres de micro économie de l’assurance, ou sur les valeurs extrêmes ! On avait mis ensemble nos notes de cours, et on a pris plaisir à travailler ensemble sur ce projet. Et le livre reste pertinent pour nos cours. Il faudra juste aller voir ailleurs pour aller plus loin sur certains aspects ! Donc toutes nos excuses d’avoir juste demandé un retirage du livre, et de ne pas avoir pris le temps de faire une révision complète du livre…

Newton-Raphson avec un dessin

Dans un précédant billet, publié il a quelques mois maintenant, je trouvais dommage que l’on oublie de citer des personnes qui ont introduit des concepts fondamentaux. En l’occurence, je pensais à Bruno de Finetti. En traînant sur Google Books l’autre jour, j’ai eu l’impression inverse, en trouvant dommage que l’on ne se souvienne que de la personne qui a introduit un concept, et pas celui qui en a compris la portée générale. Mais je devrais peut-être revenir un peu sur le contexte.

En travaillant avec un étudiant cet été, on devait résoudre des problèmes d’optimisation (en l’occurence l’estimation de paramètres dans un mélange de lois par maximum de vraisemblance) et je me souviens encore lui avoir dit, “c’est facile, tu programmes une descente de gradient, la méthode de Newton-Raphson, et c’est bon, on réglera l’histoire des contraintes plus tard” (oui, on avait un ensemble – simple – de contraintes qu’on pouvait intégrer en reparamétrant le programme d’optimisation). La méthode de Newton repose sur l’idée que pour trouver le zéro d’une fonction, on utilise une suite définie par récurence,

La fameuse illustration de cette formule étant un dessin de la forme

De là, si on se souvient que la condition du premier ordre pour cherche un extremum d’une fonction revient à annuler la dérivée, on peut utiliser cette méthode pour chercher un extremum, avec une suite définie par

Cette méthode se généralise en dimension quelconque, avec la méthode dite de la descente de gradient

Bref, ces méthodes sont très géométriques. Et j’avais toujours cru que les vieux mathématiciens étaient tous des géomètres. Aussi, j’espérais trouver plein de dessins dans les textes de Newton et Raphson. Mais non….

La méthode de Newton a été introduite par Isaac Newton dans De analysi per aequationes numero terminorum infinitas, écrit en 1669 et publié en 1711. Mais c’est surtout dans De metodis fluxionum et serierum infinitarum, écrit en 1671, traduit (du latin) et publié sous le titre Methods of Fluxions en 1736 que l’on va trouver l’idée de la méthode dite de Newton pour trouver les zéros d’une fonction à valeurs réelles. C’est ce texte que l’on peut trouver en ligne, numérisé. Mais Newton n’a aucunement proposé une méthode générale, car il ne l`appliquait qu’aux seuls polynômes. A sa décharge, la notion de dérivée n’était pas (clairement) définie à cette époque… donc imaginer une descente en suivant la tangente, en 1710, ne serait peut-être pas raisonnable.

C’est Thomas Simpson qui proposa de généraliser cette méthode pour mettre en place une méthode de calcul itératif des solutions d’une équation non linéaire, en utilisant les dérivées (qu’il appelait fluxions, comme Newton).  En 1690,  Raphson publie une description de sa méthode dans Analysis aequationum universalis. Comme Newton, il va proposer une méthode de calcul récursif des approximations successives d’un zéro d’un polynôme (avec une petite nuance, car Newton essayait de construire une suite de polynômes, comme l’expliquent Ypma (1995) ou Deuflhard (2005))

Maintenant, si on regarde les publications, le dessin avec la tangente, on ne le voit pas ! Loin de là ! On voit des pages de calculs algébriques,

Il m’a fallu un certain temps pour découvrir un traité abordant ce problème avec des dessins. Pour les amateurs de mathématiques marseillaises, on peut découvrir un traité rafraîchissant, datant de 1768, publié d’ailleurs par un homme qui allait devenir par la suite maire de Marseille (à l’époque où les hommes politiques faisaient des maths, comme j’en parlais dans un vieux billet, avec une preuve du théorème de Pythagore proposé par un mathématicien amateur, qui allait devenir ensuite Président des États-Unis d’Amérique), monsieur Jean-Raymond Pierre  Mouraille.

Ce traité est passionnant, avec un style d’une autre époque

Mourailles se positionne clairement dans la lignée de Newton, avec quelques bémols toutefois,

Le plus intéressant, ce sont les pages d’annexes, remplies de dessins, où enfin, on peut découvrir le dessin avec la suite définie par récurrence, et la descente le long de la tangente,

Bref, tout le monde pense à ces dessins quand il entend parler de la méthode de Newton-Raphson, même si aucune représentation géométrique n’est présentée.

Je ne sais pas si Jean-Raymond Pierre  Mouraille a le premier a avoir représenté la recherche de zéros d’une fonction de cette manière, mais il est dommage qu’on n’en parle pas plus souvent…