Category Archives: Geometry

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

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

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

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

Random points on the Earth

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

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

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

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

So the question here is simply : is

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

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

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

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

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

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

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

while, for another sample of points, it was possible

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

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

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

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

But there is a simple technique.

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

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

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

So the question here is simply : is

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

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

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

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

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

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

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

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

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

Probabilité et géométrie

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

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

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

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

La démonstration est simple

soit

ou enfin

tout simplement.

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

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

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

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

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

  • Rappels de géométrie

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

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

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

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

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

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

qui va impliquer l’inégalité triangulaire

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

c’est à dire que

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

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

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

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

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

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

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

Techniquement, je pense que

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

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

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

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

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

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

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

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

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

ce qui se traduit par la relation bien connue

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

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

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

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

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

et pour les séries temporelles stationnaires,

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

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

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

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

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

  • et la constante ?

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

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

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

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

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

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

et la corrélation comme

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

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

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

Ce qui se traduit par

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

si les variables sont orthogonales.

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

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

et

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

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

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

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

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

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

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

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

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

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

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

Newton-Raphson avec un dessin

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

https://latex.codecogs.com/gif.latex?T=\frac{1}{8}\left(\begin{matrix}1&%200%20&%201%20\\%200%20&%204%20&%200%20\\%201%20&%200&1\end{matrix}\right)
Here, in all areas with mass, we spread it uniformly (say), i.e. the support of https://latex.codecogs.com/gif.latex?T(C^\perp) is the one below, i.e. https://latex.codecogs.com/gif.latex?1/8th of the mass is located in each corner, and https://latex.codecogs.com/gif.latex?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 https://latex.codecogs.com/gif.latex?[0,1/4]https://latex.codecogs.com/gif.latex?[1/4,3/4] and https://latex.codecogs.com/gif.latex?[3/4,1],

Then the idea, then, is to consider https://latex.codecogs.com/gif.latex?T^2=\otimes^2T, where  https://latex.codecogs.com/gif.latex?\otimes^2T is the tensor product (also called Kronecker product) of https://latex.codecogs.com/gif.latex?T with itself. Here, the support of https://latex.codecogs.com/gif.latex?T^2(C^\perp) is

Then, consider https://latex.codecogs.com/gif.latex?T^3=\otimes^3T, where https://latex.codecogs.com/gif.latex?\otimes^3T is the tensor product of https://latex.codecogs.com/gif.latex?T with itself, three times. And the support of https://latex.codecogs.com/gif.latex?T^3(C^\perp) is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}=(a_{i,j}), then

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}\otimes\mathbf{B}%20=%20\begin{pmatrix}%20a_{11}%20\mathbf{B}%20&%20\cdots%20&%20a_{1n}\mathbf{B}%20\\%20\vdots%20&%20\ddots%20&%20\vdots%20\\%20a_{m1}%20\mathbf{B}%20&%20\cdots%20&%20a_{mn}%20\mathbf{B}%20\end{pmatrix}

So, we need a transformation matrix: consider the following https://latex.codecogs.com/gif.latex?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 https://latex.codecogs.com/gif.latex?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 https://latex.codecogs.com/gif.latex?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

http://freakonometrics.blog.free.fr/public/latex/thales.gif

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…