# R0 and the exponential growth of a pandemic

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]]

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$

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

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?

# Enveloppe convexe de points tirés au hasard

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, &gt; 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 ? # 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 Anthropométrie, ou mesure des différentes faculté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).

# 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 $n$ points uniformly distributed on a sphere. What is the probability that the $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 $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 $n$. But it was difficult. I mean, given $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 $n$ points on the unit sphere, which simply means generate $n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$
2. Try to find $\theta\in[0,2\pi]$ such that, after a rotation (with angle $\theta$) all the points lie in the upper part (the North hemisphere, for instance)

So the question here is simply : is

$\max_{\theta\in[0,2\pi]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\geq 0)\right\}$

equal to $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 $\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 $n$ points on the unit sphere, which simply means generate $2n$ random variables $U_1,\ldots ,U_n\sim\mathcal{U}([0,2\pi])$ and $V_1,\ldots,V_n\sim\mathcal{U}([-\pi/2,\pi/2])$
2. Try to find $\theta$ and $\varphi$ such that, after a rotation (with angles $\theta$ and $\varphi$) all the points lie in some given part (say  $\{x\geq0\}$)

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

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

So the question here is simply : is

$\max_{(\theta,\varphi)\in[0,2\pi]\times[-\pi/2,\pi/2]}\left\{\sum_{i=1}^n \boldsymbol{1}(\sin(U_i+\theta)\cos(V_i+\varphi)\geq 0) \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 $d$, the probability that the $n$ points (uniformly distributed on a sphere) lie on a same hemisphere is exactly

$p(d,n)= \frac{1}{2^{n-1}}\sum_{i=0}^{d-1} \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

$p(d,d+1)=1 -\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

$\mathbb{P}(Y=y)=\sum_x \mathbb{P}(Y=y,X=x)$

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

$\mathbb{P}(Y=y)=\sum_x \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)$

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

$\mathbb{E}(Y)=\sum_x \mathbb{E}(Y\vert X=x)\cdot\mathbb{P}(X=x)$

que l’on écrit souvent sous la forme raccourcie $\mathbb{E}(Y) =\mathbb{E}(\mathbb{E}(Y\vert X))$ (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

$\mathbb{E}(Y)=\sum_y y\cdot \mathbb{P}(Y=y)=\sum_y y \cdot \sum_x \mathbb{P}(Y=y, X=x)$

soit

$\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)=\sum_x\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)\cdot\mathbb{P}(X=x)$

ou enfin

$\sum_x \mathbb{P}(X=x)\sum_y y\cdot \mathbb{P}(Y=y\vert X=x)=\sum_x \mathbb{P}(X=x)\cdot \mathbb{E}(Y\vert X=x)$

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 $Y$, et $X$ le sexe, alors

$\mathbb{E}(Y)= \mathbb{E}(Y\vert X=\text{H})\cdot\mathbb{P}(X=\text{H})+\mathbb{E}(Y\vert X=\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

$\text{Var}(Y)=\text{Var}(\mathbb{E}(Y\vert X))+\mathbb{E}(\text{Var}(Y\vert X))$

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 $\mathcal E$, c’est défini par

• une notion de symétrie $<\vec x,\vec y>=<\vec y,\vec x>$
• de bilinéarité $<\alpha \vec x+\beta \vec y,\vec z>=\alpha<\vec x,\vec z>+\beta<\vec y,\vec z>$
• de positivité, $<\vec x,\vec x> \geq 0$
• et le produit scalaire est dit défini $<\vec x,\vec x>=0$ implique $\vec x=\vec 0$

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

$\|\vec x\|=\sqrt{<\vec x,\vec x>}$

En effet, on retrouve la propriété d’homogénéité, $\|\lambda\vec x\|=\vert\lambda\vert \|\vec x\|$ (c’est pour ça que l’on prend la racine carrée), et $\|\vec x\|=0$ si et seulement si $\vec x=\vec 0$. On a aussi l’inégalité de Cauchy-Schwarz,

$\|<\vec x,\vec y>\| \leq \|\vec x\|\|\vec y\|$

qui va impliquer l’inégalité triangulaire

$\|\vec x+y\|\leq \|\vec x\|+\|\vec y\|$

On dira que deux vecteurs sont orthogonaux, noté $\vec x\perp\vec y$ si $<\vec x,\vec y>=0$. Dans ce cas, on a la relation de Pythagore,

$\|\vec x+y\|^2= \|\vec x\|^2+\|\vec y\|^2$

Autre objet intéressant, la projection orthogonale. Rappelons tout d’abord que si $\mathcal F\subset \mathcal E$, on dira que $\vec x\perp\mathcal F$ si $\vec x\perp\vec y$y pour tout $\vec y\in\mathcal F$. Et on notera

$\mathcal F^\perp=\{\vec x\in\mathcal E;\vec x\perp\mathcal F\}$

On a alors des résultats intéressants. En particulier, pour tout $\vec y\in\mathcal E$, il existe un unique $\vec y_\star\in\mathcal F\subset \mathcal E$ tel que

$\|\vec y-\vec y_\star\|= \inf\{\|\vec y-\vec z\| ,\vec z\in\mathcal F\}$

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

$\Pi_{\mathcal{F}}(\vec x)=\Pi_{\mathcal{F}}\big(\Pi_{\mathcal{G}}(\vec x)\big)$

Cette relation est vraie dans $\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 – $L_2$ – peut être muni d’une telle opération, de produit scalaire,

$=\mathbb{E}(XY)$

La norme est alors $\|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 $=0$ si et seulement si $X=0$, et pour avoir une norme il faudrait que $\|X\|=0$ si et seulement si $X=0$. Le soucis technique ici est que $\mathbb{E}(X^2)=0$ signifie que $\mathbb{P}(X=0)=1$ et pas $X(\omega)=0$ pour tout $\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, $L_2$ est précisément l’espace $\mathcal L_2$ (classique en probabilité) quotienté par cette relation d’équivalence $\overset{p.s.}{=}$.

Un sous-espace de $L_2$ est l’espace des constante, $\mathbb{R}^n$, que je noterais aussi $s\{\boldsymbol{1}\}$. On notera que, comme rappelé dans tous les cours que je donne

$\mathbb{E}(Y) =\underset{c\in\mathbb{R}}{\text{argmin}}\{\mathbb{E}([Y-c]^2)\}$

c’est à dire que

$\mathbb{E}(Y)=\Pi_{s\{\boldsymbol{1}\}}(Y)$

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

$\text{Var}(Y) =\underset{c\in\mathbb{R}}{\text{min}}\{\mathbb{E}([Y-c]^2)\}$

Un autre sous-espace de $L_2$ qui sera intéressant est le suivant. Si $X\in L_2$, on notera

$s\{X\}=\{\tilde X=\psi(X)\in L_2,\psi:\mathbb{R}\rightarrow\mathbb{R}\}$

Par exemple, on peut considérer une variable de Bernoulli, $X=\boldsymbol{1}_A$$A\subset \mathbb{R}$. Dans ce cas, le sous-espace est équivalent au sous-espace des combinaisons linéaires,

$sl\{X\}=\{\tilde X=\beta_0+\beta_1 X;(\beta_0,\beta_1)\in \mathbb{R}^2\}\subset L_2$

Techniquement, je pense que

$\overline{\sigma\{X\}} \subset s\{X\}$

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

Si $Y\in L_2$, on posera

$\mathbb{E}(Y\vert X)=\Pi_{s\{\boldsymbol{X}\}}(Y)$

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

$\mathbb{E}(Y\vert X) =\underset{X_\star\in s\{X\}}{\text{argmin}}\{\mathbb{E}([Y-X_\star]^2)\}$

On notera que $\mathbb{E}(Y\vert X) \in s\{X\}}$, et on retrouve ici l’écriture standard d’un modèle de régression, $\mathbb{E}(Y\vert X) =\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éoreme de double projection tout a l’heure. Et comme $s\{\boldsymbol{1}\}\subset s\{X\}$, alors

$\Pi_{s\{\boldsymbol{1}\}}(Y)=\Pi_{s\{\boldsymbol{1}\}}\big(\Pi_{s\{X\}}(Y)\big)$

ce qui se traduit par la relation bien connue

$\mathbb{E}(Y)= \mathbb{E}\big(\mathbb{E}(Y\vert X)\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 $s\{X\}$. On introduit souvent l’opérateur d’espérance linéaire, avec

$\text{EL}(Y\vert X)=\Pi_{sl\{{X}\}}(Y)$

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

$sl\{X\}=\{\tilde X=\beta_0+\beta_1 X;(\beta_0,\beta_1)\in \mathbb{R}^2\}\subset L_2$

et pour les séries temporelles stationnaires,

$\text{EL}(X_t\vert X_{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 $L_2$ est un espace de Hilbert, ce n’est pas le cas de $L_1$. Pourtant, nul besoin d’être dans $L_2$ pour définir une espérance conditionnelle… Je pense que l’astuce peut être de noter que si $X$ appartient à $L_1$ mais pas $L_2$, on peut quand même s’en sortir car $L_1$ est un espace de Banach et on a malgré tout des notions de convergences. Et $L_2$ est dense dans $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 $X\perp Y$ se traduit ici par $\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

$\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 $\vec x=X-\mathbb{E}(X)$, et $\vec y=Y-\mathbb{E}(Y)$, c’est à dire que l’on va travailler sur les variables centrées. On peut noter que

$\text{Var}(X)=\mathbb{E}([X-\mathbb{E}(X)]^2)=\| \vec x \|^2$

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

$\text{Cov}(X,Y)=<\vec x,\vec y>$

et la corrélation comme

$\text{corr}(X,Y)=\frac{<\vec x,\vec y>}{\| \vec x \|\cdot \| \vec y \|} =\cos(\theta)$

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

$\| \vec x +\vec y \|^2=\| \vec x \|^2+\| \vec y \|^2$

$\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

$\vec x=\mathbb{E}(Y\vert X)-\mathbb{E}(Y)$

et

$\vec y=Y-\mathbb{E}(Y\vert X)$

On note facilement que $\vec x \perp \vec y$, car $\vec x\in s\{X\}$ alors que $\vec y\in s\{X\}^\perp$, par construction de la projection orthogonale. Aussi le théorème de Pythagore nous dit que

$\| Y-\mathbb{E}(Y) \|^2=\| \mathbb{E}(Y\vert X)-\mathbb{E}(Y) \|^2+\| Y-\mathbb{E}(Y\vert X) \|^2$

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

$\| Y-\mathbb{E}(Y) \|^2= \text{Var}(Y)$

Pour le premier terme de droite, la aussi, c’est facile, car $\mathbb{E}[\mathbb{E}(Y\vert X)]=\mathbb{E}(Y)$, par le théorème de double projection, et donc

$\| \mathbb{E}(Y\vert X)-\mathbb{E}(Y) \|^2=\text{Var}(\mathbb{E}(Y\vert X))$

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.

$\| \vec Y-\mathbb{E}(Y\vert X) \|^2 = \mathbb{E}(\text{Var}(Y\vert X))$

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 $sl\{\boldsymbol{X}\}$, on aura donc une partie de la variance qui sera expliquée par nos variables explicatives $\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)

# 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,

$a_{k+1} = a_k - \frac{f(a_k)}{f'(a_k)}$

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

$a_{k+1} = a_k - \frac{f'(a_k)}{f''(a_k)}$

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

$\boldsymbol{a}_{k+1} = \boldsymbol{a}_k - [H f(\boldsymbol{a}_k)]^{-1} \nabla f(\boldsymbol{a}_k)$

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

# Pavage et magie, rotations et symétries

(oui, mon fils voulait une rime dans le titre… promis on va tenter de faire une comédie musicale à partir de mes billets). La semaine passée, un article intéressant était publié (http://www.bbc.co.uk/news/…) sur sciences et magie. Cet article m’a fait penser à une réflexion que je m’étais faite lors de notre visite au MoMaths, cet hiver. Pour revenir sur l’histoire que j’avais raconté il y a quelques mois, sur la fin de la visite, les enfants avaient joué avec les pavages, à la Escher (les pièces étaient magnétiques, et s’assemblaient sur un grand tableau blanc)

Le soucis c’est qu’à aucun moment je n’ai pu voir d’explication du  pourquoi ça marche ? ou pour être plus explicite comment a été construit le singe pour que tout s’emboîte aussi bien. Car il y a un truc… La magie ça n’existe pas ! ce sont juste des maths (c’est un peu ce que racontent Persi Diaconis et Ron Graham dans Magical Mathematics, même si j’avoue ne pas l’avoir fini, faute de temps, et d’avoir laissé des chapitres en attente). Mon fils travaillant les rotations et les symétries, à l’école, j’ai voulu qu’on prenne le temps de revenir sur les pavages. Mais avant de me lancer avec lui, j’ai essayé de m’entraîner un peu tout seul (dimanche après midi, pendant le temps calme).

La difficulté (ou la magie) est qu’il y a des symétries dans tous les sens. Regardons un exemple simple (ou disons pas trop compliqué) avec des chinois (identifiés par les chapeaux dits chinois). En regardant rapidement, on note que le triangle ci-dessous est reproduit à l’infini, par un jeu de rotations. Il suffit donc d’avoir un tel triangle. Le soucis est que le découpage (parfaitement symétrique) n’est pas trivial (il est délicat de faire un personnage en faisant un découpage dans un triangle, même si c’est la méthode retenue dans l’activité mentionnée sur https://irem.univ-lille1.fr/…)

Le secret part de l’observation suivante: sur le graphique de gauche, en dessous, les trois points (qui sont des intersections des trois personnages, ou disons des trois couleurs) sont les sommets d’un triangle équilatéral. C’est là que la magie va opérer, car c’est cette figure (avec une autre, juxtaposée) que l’on va tenter de reproduire. En fait, un personnage (disons un petit chinois bleu) c’est deux triangles collés (on doit pouvoir parler de parallélogramme isocèle).

On ne le voit pas très bien, mais pourtant, au départ, on avait un parallélogramme isocèle (dont les contours sont tracés en noir) qui était  tout bleu. Bien entendu, on a découpé des parties et en les a envoyé ailleurs (par rotation autour d’un sommet). Le plus simple est de construire la figure, à la main. On part donc du parallélogramme suivant, tout bleu :

Oui, j’ai repris les fonctions de base utilisées l’autre jour pour expliquer R aux enfants,

source("http://freakonometrics.free.fr/RforKIDZ.R")

En fait, j’ai souhaité définir une fonction créant un polynôme à partir d’une liste de points (ou de sommets).

polyg=function(L,couleur="blanc",tour=couleur,ep=1,st=1,...){
b=unlist(L)
n=length(b)
polygon(b[seq(1,n,by=2)],b[seq(2,n,by=2)],col=couleurc(couleur),border=couleurc(tour),lwd=ep,lty=st,...)}

La fonction pour faire une rotation d’une liste de points est ici

rotation=function(L,angle=90,centre=L[[1]]){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]-centre[1],b[seq(2,n,by=2)]-centre[2])
Y=matrix(c(cos(angle/180*pi),sin(angle/180*pi),-sin(angle/180*pi),cos(angle/180*pi)),2,2)%*%t(X)
Z=Y+centre
L2=L
for(i in 1:(n/2)){L2[[i]]=Z[,i]}
return(L2)}

On a alors ici les quatre sommets

A=c(4,0)
B=c(4,10)
C=c(4+sqrt(3)*5,5)
D=c(4+sqrt(3)*5,15)

et on peut tracer le parallélogramme avec

L1=list(A,B,D,C)
polyg(L1,"bleu",tour="noir",ep=3)

On peut alors jouer à découper des parties, puis à les pivoter (on va faire des rotations, quoi). Par exemple, le prend le triangle en haut à droite, et je pivote de 240 degrés autour d’un des sommets du parallélogramme. C’est la figure ci-dessus à droite.  Le triangle, il est défini par le code suivant

E=c(10,10)
F=c(7,10+sqrt(3))
L1=list(B,E,F)
polyg(L1,tour="rouge",ep=3)

ce qui enlève la portion – en traçant un triangle blanc – et pour faire la rotation, j’utilise le code suivant

L2=rotation(L1,240)
polyg(L2,"bleu","rouge",ep=3)

(la seconde ligne permet de tracer le triangle bleu). Il va de soi qu’on peut découper des objets plus complexes que des triangles, voire non connexes. Par exemple, je peux découper la forme verte ci-dessous, puis la pivoter par rapport au sommet en haut à gauche, de 240 degrés.

G=c(10,2*sqrt(3)+10)
H=c(4+sqrt(3)*5,2*sqrt(3)+10)
F1=c(8,10+2/sqrt(3))
F2=c(8,10+4/sqrt(3))
L1=list(B,E,F1,F2,B,G,H,D,G)
polyg(L1,"vert",tour=NA)
L2=rotation(L1,240)
polygL(L2,"bleu",tour=NA)

(il n’est pas nécessaire de revenir au point B, mais peu importe). Ou la figure jaune, qui va pivoter autour du sommet en bas à droite, de 120 degrés.

F3=c(8,10-2/sqrt(3))
F4=c(8,10-4/sqrt(3))
G2=c(10,-2*sqrt(3)+10)
H2=c(4+sqrt(3)*5,-2*sqrt(3)+10)
L1=list(H,G,F2,F1,E,F3,F4,G2,H2,H)
polyg(L1,"jaune",tour=NA)
L2=rotation(L1,120,centre=C)
polygL(L2,"bleu",tour=NA)

Avec des ciseaux, ça sera plus simple, mais le but est aussi de voir que programmer, c’est facile. En gros, on a deux fonctions de base: une qui trace un polynôme et l’autre qui fait des rotations.

Si on fait les deux en même temps (à peu près), on enlève toute la partie en haut à droite, et on la répartie: une partie à gauche, et une partie en bas à droite.  On a alors une forme qui ressemble à un chinois. Non ? On peut aller creuser un peu sous les bras, pour améliorer les mains, mais on y est presque, me semble-t-il…

On peut continuer à améliorer, par exemple en enlevant des parties, puis en faisant des rotations par rapport à un des sommets initiaux. Comme le petit triangle jaune, qui va pivoter de 120 degrés par rapport au sommet en bas à droite. Puis le triangle vert par rapport au sommet en bas à gauche (on note qu’on pivote par rapport aux sommets du triangle équilatéral qu’on avait repéré, initialement)

S=c(8-.4,-2*sqrt(3)+10-.4)
L1=list(F4,G2,S,F4)
polyg(L1,"jaune")
L2=rotation(L1,120,centre=C)
polyg(L2,"bleu")

Bon, là on fait une pause, parce qu’on a, par morceaux, notre chinois ! On un paquet de lettre (voire un paquet de sommets qu’il va falloir nommer, suite à nos rotations). On les nomme avec un code qui ressemble à (oui, on aurait pu pu les nommer dès le départ, lors des rotations)

rE=rotation(list(E),240,centre=B)[[1]]
rF1=rotation(list(F1),240,centre=B)[[1]]
rH=rotation(list(H),240,centre=B)[[1]]

Une fois qu’on a identifié nos sommets, on a notre chinois,

chinois=list(B,rE,rF1,LG,S2,LF2,rH,A,r2H,r2G,r2S2,
r2F2,r2F1,r2E,rF3,rF4,r2S,rG2,rH2,C,H2,G2,S,F4,F3,E,B)

Bon, ensuite, reste à paver, ce qui est est le plus simple. Par exemple, on peut faire des rotations,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")

voire définir une translation,

translation=function(L,h){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]+h[1],b[seq(2,n,by=2)]++h[2])
L2=L
for(i in 1:(n/2)){L2[[i]]=X[i,]}
return(L2)}

puis translater notre chinois,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")
polyg(translation(r1chinois,15/sqrt(3),15),couleur="jaune")

Il manque l’autre, le vert. On peut vérifier qu’il se superpose parfaitement :

Bon, ben ça y est, on a constitué notre pavage. Reste à travailler un peu le code, essentiellement sur la forme, et on pourra jouer avec les enfants, à faire les lézards, des singes, et toutes sortes de pavages possibles, à la Escher. En notant qu’on peut utiliser d’autres formes de bases plus complexes, comme expliqué sur http://xavier.hubaut.info/…

# Fractals and Kronecker product

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

$T=\frac{1}{8}\left(\begin{matrix}1& 0 & 1 \\ 0 & 4 & 0 \\ 1 & 0&1\end{matrix}\right)$
Here, in all areas with mass, we spread it uniformly (say), i.e. the support of $T(C^\perp)$ is the one below, i.e. $1/8$th of the mass is located in each corner, and $1/2$ is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals $[0,1/4]$$[1/4,3/4]$ and $[3/4,1]$,

Then the idea, then, is to consider $T^2=\otimes^2T$, where  $\otimes^2T$ is the tensor product (also called Kronecker product) of $T$ with itself. Here, the support of $T^2(C^\perp)$ is

Then, consider $T^3=\otimes^3T$, where $\otimes^3T$ is the tensor product of $T$ with itself, three times. And the support of $T^3(C^\perp)$ is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if $\mathbf{A}=(a_{i,j})$, then

$\mathbf{A}\otimes\mathbf{B} = \begin{pmatrix} a_{11} \mathbf{B} & \cdots & a_{1n}\mathbf{B} \\ \vdots & \ddots & \vdots \\ a_{m1} \mathbf{B} & \cdots & a_{mn} \mathbf{B} \end{pmatrix}$

So, we need a transformation matrix: consider the following $4\times4$ matrix,

> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1

Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a $4^2\times4^2$ matrix,

> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1

And then, we continue,

> for(s in 1:3){N=kronecker(N,M)}

After only a couple of loops, we have a $4^5\times4^5$ matrix. And we can plot it simply to visualize the support,

> image(N,col=c("white","blue"))

As we zoom in, we can visualize this fractal property,

# 11/11/11 11:11

Aujourd’hui, c’est la fête du parallélisme ! Déjà que 11:11 mettait Pierre Desproges dans un état proche de l’orgasme, mais qu’il soit 11:11 le 11/11/11… c’est la fête ! Tous les amateurs de parallélisme vont apprécier.

On va aussi dire que c’est la fête du théorème de Thalès.

Monsieur Thalès, alias Thalès de Millet (Θαλῆς ὁ Μιλήσιος) était un brave monsieur qui faisait des mathématiques il y a un peu plus de 2500 ans (et contrairement à ce que croit ma fille, je n’ai pas connu cette époque là, ni celle des – vrais – dinosaures). Selon Plutarque, dans Le Banquet des Sept Sages, un pharaon (en l’occurrence Ahmès II) aurait dit que personne ne serait capable de calculer la hauteur de la grande pyramide de Khéops. “Challenge accepted” se serait alors écrié Thalès (πρόκληση αποδεκτή à l’époque) et il aurait relevé le défi en calculant le rapport entre son ombre et celle d’un bout de bois…

vous avez mesuré la pyramide sans le moindre embarras et sans avoir eu besoin d’aucun instrument. Après avoir dressé votre bâton à l’extrémité de l’ombre que projetait la pyramide, vous construisîtes deux triangles par la tangence d’un rayon, et vous démontrâtes qu’il y avait la même proportion entre la hauteur du bâton et la hauteur de la pyramide qu’entre la longueur des deux ombres

(Denis Guedj en parle dans le Théorème du Perroquet). Formellement, si on utilise le dessin à droite, ce théorème dit que si les droites rouges sont parallèles, alors

Ce théorème est fondamental non seulement en géométrie, mais plus généralement, car c’est la base de la règle de trois ! En effet, le parallélisme induit de la proportionnalité. Et la proportionnalité est la base de tout. J’avais fait un billet sur la règle de trois et les techniques de provisionnement en assurance non-vie, et j’en ai parlé dans le cours de tarification, lorsque l’on faisait des régressions avec un modèle multiplicatif. Mais je parle trop, cet instant magique ne dure qu’une seconde, et ça serait dommage de le gâcher avec mes petites histoires…