Category Archives: Maths

Geometry, Maths

R0 and the exponential growth of a pandemic

Leave a comment

For some dissemination work, I want to create a nice graph to explain the exponential growth in pandemics, related to the value of R_0. Recall that R_0 corresponds to the average number of people that a contagious person can infect. Hence, with R_0=1.5, 4 people will contaminate 6 people, and those 6 will contaminate 9, etc. After n iteration, the number of contaminated people is simply R_0{}^n. As explained by Daniel Kahneman

people, certainly including myself, don’t seem to be able to think straight about exponential growth. What we see today are infections that occurred 2 or 3 weeks ago and the deaths today are people who got infected 4 or 5 weeks ago. All of this is I think beyond intuitive human comprehension

For different values of R_0 (on each row), I wanted to visualise the number of contaminated people after 3, 5 or 7 iterations, since graphs are usually the most simple way to give some intuition. The graph I had in mind was the following

(to be honest, I am quite sure I had seen it somewhere, but I cannot find where). The main challenge here is pack optimally k identical circles intro a unit circle: we need here the location of the points (center of the disks) and the radius. It seems to be a rather complicated mathematical problem. Nicely, on http://hydra.nat.uni-magdeburg.de/packing, it is possible to get the “best known packings of equal circles in a circle” (up to 5000, but many k‘s are missing). For instance, for k=37, we have

And interestingly, on the same website, we can get the coordinates of the centers, for example with 37 disks, so it is possible to recreate the R graph.

k = 37
base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE)

The problem, as discussed earlier, is that some cases are not solved, yes, for instance k=2^{12}=4096: the next feasable case is 4105. To avoid that issue, one can use

T = "Error"
while(T == "Error"){
    T = substr(try(base = read.table(paste("http://hydra.nat.uni-magdeburg.de/packing/cci/txt/cci",k,".txt",sep=""), header=FALSE),silent = TRUE),1,5)
k=k+1
} 
k=k-1

Now we can almost plot it. The problem is that the radius of the circles is missing, here. But we can compute it

D=as.matrix(dist(x = base[,2:3]))
diag(D)=1e5
i=which(D == min(D), arr.ind = TRUE)
r = D[i[1,1],i[1,2]]

Here the radius is

r
[1] 0.2959118

To plot it, use

plot(base$V2,base$V3,xlim=c(-1,1),ylim=c(-1,1))
n=100
theta=seq(0,pi,length=n+1)
circ= function(x,y,r,h=1){
  vu=x+r*cos(theta)
  vv=r*sin(theta)
  cbind(c(vu,rev(vu))*h,c(y+vv,y-rev(vv))*h)
}
for(i in 1:k) polygon(circ(base[i,2],base[i,3],r/2*.95),col=colr,border=NA)

We can now use that code to create the graph above, with k=R_0{}^n for various values of n

And we can also use it to visualize more subtile differences, like R_0=1.1, R_0=1.3, R_0=1.5 and R_0=1.7

Geometry, Regression

Régression sur des variables corrélées, un peu de géométrie

Leave a comment

Depuis quelques jours, je mets en ligne des billets sur la régression sur des variables corrélés, en commençant par du deuxième effet kiss-cool suivi d’un court billet sur des régressions en cascade. Oui, depuis le premier, mon collègue Olivier ne cesse de demander des compléments. Alors ce soir, on va faire des dessins…

On a ici 3 vecteurs, \color{red}{\vec{x_1}}, \color{green}{\vec{x_2}} et \color{blue}{\vec{y}}. Pour rappel, la longueur de mes vecteurs, c’est la variance de la variable associée, \|{\color{red}{\vec{x_1}}}\|=\text{Var}({\color{red}{x_1}}). Et commençons par supposer que les variables \color{red}{x_1} et \color{green}{x_2} sont indépendantes (ou plus simplement non-corrélées) soit {\color{red}{\vec{x_1}}}\perp\color{green}{\vec{x_2}}. C’est le dessin ci-dessous,

La régression linéaire multiple est ici équivalente aux régressions simples, comme on en parlait dans les billets précédant. Géométriquement, on voit que les coordonnées de la projection de \color{blue}{\vec{y}} sur l’espace engendré par \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} (c’est à dire le plan (x,y)) correpond au point noir. C’est la régression multiple. Les deux régressions simples sont obtenus en projetant respectivement sur sur l’espace engendré par \color{red}{\vec{x_1}} (c’est à dire la droite (x)) et sur \color{green}{\vec{x_2}} (c’est à dire la droite (y)).

Supposons maintenant que \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} ne sont plus orthogonaux. Plus précisément, on va garder \color{red}{\vec{x_1}} et \color{blue}{\vec{y}} inchangés, et on va changer , \color{green}{\vec{x_2}}, tout en gardant la norme \|\color{green}{\vec{x_2}}\| inchangée. C’est ce qu’on a sur le schéma ci-dessous

On avait vu qu’on ne pouvait plus faire deux régressions indépendantes, mais qu’on pouvait garder la première, et en faire ensuite une seconde, mais en faisant un peu attention. Pour la première régression, on va projeter \color{blue}{\vec{y}} sur \color{red}{\vec{x_1}} et ca c’est facile, rien n’a changé,

Pour la seconde étape, rappelons qu’on doit projeter à la fois \color{blue}{\vec{y}} et \color{green}{\vec{x_2}} sur l’orthogonal de \color{red}{\vec{x_1}}, c’est à dire ici le plan (y,z), en rouge ci-dessous

(je renvois au précédant billet pour les formules algébriques). On obtient alors les projections suivantes, qu’on pourrait noter \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} et \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}}

On le voit, si on regarde sur le mur du fond, en rouge, la projection de \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} sur \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}} est un peu plus grande, autrement dit, le coefficient de la régression simple sera un peu plus grand. En fait, si on veut faire des calculs, on voit que le sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle va intervenir ici… Plus l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle sera petit, plus les variables seront corrélées, et plus la correction sera importante. Graphiquement, la correction dont je parle, c’est le petit segment, à gauche sur le dessin

Mais simulons quelques données pour visualiser tout ça…

library(mnormt)
n=10000
r=0
set.seed(1)
Z=cbind(rnorm(n)*sqrt(2),rnorm(n)*sqrt(3))
df=data.frame(y=Z[,1]+Z[,2]+rnorm(n)*sqrt(3),x1=Z[,1],x2=Z[,2])

On simule ici des variables explicatives indépendamment. Et on peut vérifier que les distances sont bien celles de notre dessin

var(df$x1)
[1] 2.049731
var(df$x2)
[1] 2.810419
var(df$y)
[1] 8.066665

et en plus, les variables sont quasiment orthogonales

cor(df$x1,df$x2)
[1] 0.004845078

Si on fait la régression multiple

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  0.99904    0.01219   81.94   2e-16 ***
x2  1.00090    0.01017   98.40   2e-16 ***

On a des coefficients qui sont identiques (à très peu de choses près, mais le soucis vient du fait que j’ai enlevé la constante du modèle pour coller au mieux avec mon dessin)

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

et pour la seconde régression aussi

reg2=lm(y~0+x2,data=df)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x2  1.00496    0.01315   76.42   2e-16 ***

On a déjà parlé de ce résultat longuement (je renvoie au précédant billet, des régressions en cascade). Regardons maintenant ce qui se passe si on garde nos observations \color{blue}{y_i}, \color{red}{x_{1,i}} mais que maintenant, la variable \color{red}{x_{2,i}} se retrouve corrélée avec \color{red}{x_{1,i}} , disons avec une corrélation de l’ordre de 0.4

r=.4
df$x2=r*df$x1+sqrt(1-r^2)*df$x2
cor(df$x1,df$x2)
[1] 0.3461496

Faire la régression multiple  donne ici

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1   0.5622     0.0130   43.26   2e-16 ***
x2   1.0921     0.0111   98.40   2e-16 ***

ce qui est différent des régressions simples, comme la première

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

Ce que nous dit le théorème de Frisch-Waugh, c’est qu’on peut faire des régressions en cascades: on peut garder cette première régression, mais pour la seconde, on va régresser les projections sur l’orthogonal à \color{red}{\vec{x_1}} (le mur en rouge au fond dans notre dessin). Ici, on obtient

reg2bis=lm(residuals(lm(y~x1,data=df))~residuals(lm(x2~x1,data=df)))
summary(reg2bis)
 
Coefficients:
           Estimate Std. Error t value Pr(|t|)      
residuals(2) 1.0921     0.0111    98.4   2e-16 ***

qui donne effectivement un coefficient de régression un peu supérieur à celui que nous avions auparavant, comme nous l’observions sur le dessin. En fait on peut remarquer que la différence relative s’obtient très facilement avec la corrélation

1/sqrt(1-r^2)
[1] 1.091089

et pour rappel, d’un point de vue géométrique, la corrélation est un cosinus d’un angle, r=\cos(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle). Autrement dit \sqrt{1-r^2}=\sin(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle), c’est à dire que la correction est directement liée au sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle, comme mentionné auparavant… amusant, non?

Maths, Publications

(Petite) histoire du hasard et de la simulation

Leave a comment

Entendre « il y a 10% de chances de pluie aujourd’hui » ou « le test médical a une valeur prédictive positive de 75% » montre que les probabilités sont aujourd’hui partout[i]. Une probabilité est une grandeur difficile à appréhender, mais incontournable quand on cherche à théoriser et mesurer le hasard. Et si la théorie mathématique est finalement arrivée très tard, comme le rappelle Hacking (2006), cela n’a pas empêché l’assurance de se développer assez tôt, et d’avoir les premières tables (actuarielles) de mortalité avant même que la « probabilité de décès » ou « l’espérance de vie » n’ait de fondement mathématique. Et de la même manière, de nombreuses techniques ont été inventées pour « générer du hasard », avant l’explosion des méthodes dites de Monte Carlo, en parallèle avec le développement de l’informatique (et du fait qu’une machine pouvait générer du hasard). Continue reading (Petite) histoire du hasard et de la simulation →

Geometry

Enveloppe convexe de points tirés au hasard

Leave a comment

Le week-end dernier, Jean-Baptiste qui était de passage à la maison, me présentait un problème amusant de géométrie, lié à un papier mis en ligne l’an dernier, monotonicity of facet numbers of random convex hulls. Dans cet article, ils montrent que quand on tire n points au hasard (dans un espace de dimension d) alors P_n, le nombre moyen de face de l’enveloppe convexe est strictement croissant avec n. Si on n’a pas regardé la démonstration (on avait mieux à faire), Jean-Baptiste me disait que ce problème très simple était en fait très complexe. Et bien entendu, comme ça m’a interpelé, j’ai voulu regarder plus en détails. En dimension d=2, et en tirant en plus uniformément sur le carré unité (oui, j’ai fait très très simple).

Pour tirer des points au hasard, et récupérer l’enveloppe convexe, c’est assez simple, par exemple avec 15 points

library(sp)
library(geosphere)
n=15
UV=matrix(runif(2*n),n,2)
CH=chull(UV)
PLCH=UV[c(CH,CH[1]),]
plot(c(0,0,1,1),c(1,0,0,1))
polygon(PLCH,border="blue",col=rgb(0,0,1,.3))
points(UV,pch=19,cex=2,col="red")

On le voit, ici l’enveloppe convexe est un polygône de 8 côtés (ou 8 sommets). On peut alors faire un petit code pour tirer des points au hasard, construire l’enveloppe convexe, et sortir des infos (nombre de points extrémaux, surface de l’enveloppe convexe, présence ou non de certains points – sur la diagonale, etc)

simu=function(n,isplot=FALSE){
UV=matrix(runif(2*n),n,2)
CH=chull(UV)
PLCH=UV[c(CH,CH[1]),]
nb_ex=length(CH)
p_in=function(u) point.in.polygon(u,u,PLCH[,1],PLCH[,2])
pts_in=Vectorize(p_in)(seq(.5,.95,by=.05))
if(isplot==TRUE) lines(PLCH,col=rgb(0,0,1,.25))
return(list(nb=nb_ex,area=areaPolygon(PLCH),pts=pts_in))}

par exemple

plot(c(0,0,1,1),c(1,0,0,1),col="white")
for(s in 1:1000){S=simu(5,isplot=TRUE)}

(on pourrait bien entendu stocker tout plein de choses)

ou encore avec n= 20 points au lieu de 5

plot(c(0,0,1,1),c(1,0,0,1),col="white")
for(s in 1:1000){S=simu(20,isplot=TRUE)}

Essayons de boucler un peu sur n maintenant

Np=c(3,4,5,6,7,8,10,15,20,30,40,50,75,100,200)
VN=VA=rep(NA,15)
NN=matrix(NA,20000,15)
VPT=matrix(NA,15,10)
for(i in 1:15){
N=A=rep(NA,20000)
PT=matrix(NA,20000,10)
np=Np[i]
for(s in 1:20000){
S=simu(np,isplot=FALSE)
N[s]=S$nb
PT[s,]=S$pts
A[s]=S$area
}
NN[,i]=N
VN[i]=mean(N,na.rm=TRUE)
VA[i]=mean(A,na.rm=TRUE)
VPT[i,]=apply(PT,2,function(x) mean(x,na.rm=TRUE))
}

Cette fois on stocke tout plein de choses. On peut juste faire un boxplot du nombre des points extrémaux en fonction de la taille de l’échantillon

VV=rep(Np,each=20000)
boxplot(as.vector(NN)~as.factor(VV))

Oui, en moyenne, ça semble croitre. Plus amusant, si on regarde la moyenne en fonction de \log(n)

plot(Np,VN,type="l",log="x",col="blue")

on obtient… une belle droite ! Le nombre moyen de points extrémaux croit en \log(n). On peut même avoir la pente de cette droite,

> lm(VN~log(Np))
Coefficients:
(Intercept) log(Np)
0.05224 2.58717717

Je laisse les plus courageux trouver du sens à ce 2.58717… Si on continue un peu, on peut regarder la probabilité que (u,u) soit à l’intérieur, pour plusieurs valeurs de u.

plot(Np,VPT[,10],type="l",log="x",col="blue")
lines(Np,VPT[,8])
lines(Np,VPT[,5],col="red",lwd=2)

On retrouve des fonctions croisantes, en fonction de n, mais la convexité semble dépendre de l’endroit où se trouve le point u. Amusant, non ?

Maths

Le retour des abeilles

2 Comments

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”.

Geometry, Maths

Maudites Abeilles

7 Comments

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…)

 

Geometry

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

5 Comments

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 →

Maths

Trouver le millionième chiffre

7 Comments

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 ?

Maths

What happens if we forget a trivial assumption ?

Leave a comment

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 ? →

Maths

La formation française, élitiste ?

2 Comments

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 (2012, 2009, 2006, 2003, 2000). 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.

Maths

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

Leave a comment

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) →

Copulas, MAT8595, Maths, Probability

On Hoeffding’s identity

1 Comment

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) that2\text{cov}(X,Y)=2\big(\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y]\big)can also be written

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

Maths

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

Leave a comment

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/…)

Optimization

Random points on some hemisphere

Leave a comment

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!