Tag Archives: projection

Probabilistic Foundations of Econometrics, part 2

This post is the second one of our series on the history and foundations of econometric and machine learning models. Part 1 is online here.

Geometric Properties of this Linear Model

Let’s define the scalar product in \mathbb{R}^n, ⟨\mathbf{a},\mathbf{b}⟩=\mathbf{a}^T\mathbf{b}, and let’s note \|\cdot\| the associated Euclidean standard, \|\mathbf{a}\|=\sqrt{\mathbf{a}^T\mathbf{a}} (denoted \|\cdot\|_{\ell_2} in the next post). Note \mathcal{E}_X the space generated by all linear combinations of the \mathbf{X} components (adding the constant). If the explanatory variables are linearly independent, \mathbf{X} is a full (column) rank matrix and \mathcal{E}_X is a space of dimension p+1. Let’s assume from now on that the variables \mathbf{x}  and y are centered here. Note that no law hypothesis is made in this section, the geometric properties are derived from the properties of expectation and variance in the set of finite variance variables.

With this notation, it should be noted that the linear model is written m(\mathbf{x})=⟨\mathbf{x},\beta⟩. The space H_z=\{\mathbf{x}\in\mathbb{R}^{p+1}:m(\mathbf{x})=z\} is a hyperplane (affine) that separates the space in two. Let’s define the orthogonal projection operator on \mathcal{E}_X, \Pi_X =\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T. Thus, the forecast that can be made for it is: \widehat{\mathbf{y}}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T\mathbf{y}=\Pi_X\mathbf{y}. As, \widehat{\varepsilon}=\mathbf{y}-\widehat{\mathbf{y}}=(\mathbb{I}-\Pi_X)\mathbf{y}=\Pi_{X^\perp}\mathbf{y}, we note that \widehat{\varepsilon}\perp\mathbf{x}, which will be interpreted as meaning that residuals are a term of innovation, unpredictable in the sense that \Pi_{X }\widehat{\varepsilon}=\mathbf{0}. The Pythagorean theorem is written here: \Vert \mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y} \Vert^2+\Vert \Pi_{ {X}^\perp}\mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y}\Vert^2+\Vert \mathbf{y}-\Pi_{ {X}}\mathbf{y}\Vert^2=\Vert\widehat{\mathbf{y}}\Vert^2+\Vert\widehat{\mathbf{\varepsilon}}\Vert^2which is classically translated in terms of the sum of squares: \underbrace{\sum_{i=1}^n y_i^2}_{n\times\text{total variance}}=\underbrace{\sum_{i=1}^n \widehat{y}_i^2}_{n\times\text{explained variance}}+\underbrace{\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{n\times\text{residual variance}} The coefficient of determination, R^2, is then interpreted as the square of the cosine of the angle \theta between \mathbf{y} and \Pi_X \mathbf{y} : R^2=\frac{\Vert \Pi_{{X}} \mathbf{y}\Vert^2}{\Vert \mathbf{y}\Vert^2}=1-\frac{\Vert \Pi_{ {X}^\perp} \mathbf{y}\Vert^2}{\Vert \mathbf {y}\Vert^2}=\cos^2(\theta)An important application was obtained by Frish & Waugh (1933), when the explanatory variables are divided into two groups, \mathbf{X}=[\mathbf{X}_1 |\mathbf{X}_2], so that the regression becomes y=\beta_0+\mathbf{X}_1 β_1+\mathbf{X}_2 β_2+\varepsilon. Frish & Waugh (1933) showed that two successive projections could be considered. Indeed, if \mathbf{y}_2^\star=\Pi_{X_1^\perp} \mathbf{y} and X_2^\star=\Pi_{X_1^\perp}\mathbf{X}_2, we can show that \widehat{\beta} _2=[{\mathbf{X}_2^\star}^T \mathbf{X}_2^\star]^{-1}{\mathbf{X}_2^\star}^T \mathbf{y}_2^\star In other words, the overall estimate is equivalent to the combination of independent estimates of the two models if \mathbf{X}_2^\star=\mathbf{X}_2, i.e. \mathbf{X}_2\in \mathcal{E}_{X_1}^\perp, which can be noted \mathbf{x}_1\perp\mathbf{x}_2 We obtain here the Frisch-Waugh theorem which guarantees that if the explanatory variables between the two groups are orthogonal, then the overall estimate is equivalent to two independent regressions, on each of the sets of explanatory variables. This is a theorem of double projection, on orthogonal spaces. Many results and interpretations are obtained through geometric interpretations (fundamentally related to the links between conditional expectation and the orthogonal projection in space of variables of finite variance).

This geometric interpretation might help to get a better understanding of the problem of under-identification, i.e. the case where the real model would be y_i=\beta_0+ \mathbf{x}_1^T \beta_1+\mathbf{x}_2^T \beta_2+\varepsilon_i, but the estimated model is y_i=b_0+\mathbf{x}_1^T \mathbf{b}_1+\eta_i. The maximum likelihood estimator of \mathbf{b}_1 is \widehat{\mathbf{b}}_1=\mathbf {\beta}_1 + \underbrace{ (\mathbf {X}_1^T\mathbf {X}_1)^{-1} \mathbf {X}_1^T \mathbf {X}_{2} \mathbf{\beta}_2}_{\mathbf{\beta}_{12}}+\underbrace{(\mathbf{X}_1^{T}\mathbf{X}_1)^{-1} \mathbf{X}_1^T\varepsilon}_{\nu}so that \mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1+\beta_{12}, the bias ( \beta_{12}) being null only in the case where \mathbf{X}_1^T \mathbf{X}_2=\mathbf{0} (i. e. \mathbf{X}_1\perp \mathbf{X}_2 ): we find here a consequence of the Frisch-Waugh theorem.

On the other hand, over-identification corresponds to the case where the real model would be y_i=\beta_0+\mathbf{x}_1^T \beta_1+\varepsilon_i, but the estimated model is y_i=b_0+ \mathbf{x}_1^T \mathbf{b} _1+\mathbf{x}_2^T \mathbf{b}_2+\eta_i. In this case, the estimate is unbiased, in the sense that \mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1 but the estimator is not efficient. Later on, we will discuss an effective method for selecting variables (and avoid over-identification).

From parametric to non-parametric

We can rewrite equation (4) in the form \widehat{\mathbf{y}}=\Pi_X\mathbf{y} which helps us to see the forecast directly as a linear transformation of the observations. More generally, a linear predictor can be obtained by considering m(\mathbf{x})=\mathbf{s}_{\mathbf{x}}^T \mathbf{y}, where \mathbf{s}_{\mathbf{x}} is a weight vector, which depends on \mathbf{x}, interpreted as a smoothing vector. Using the vectors \mathbf{s}_{\mathbf{x}_i}, calculated from the observations \mathbf{x}_i, we obtain a matrix \mathbf{S} of size n\times n, and \widehat{\mathbf{y}}=\mathbf{S}\mathbf{y}. In the case of the linear regression described above, \mathbf{s}_{\mathbf{x}}=\mathbf{X}[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{x}, and in that case \text{trace}(\mathbf{S}) is the number of columns in the \mathbf{X} matrix (the number of explanatory variables). In this context of more general linear predictors, \text{trace}(\mathbf{S}) is often seen as equivalent to the number of parameters (or complexity, or dimension, of the model), and \nu=n-\text{trace}(\mathbf{S}) is then the number of degrees of freedom (see Ruppert et al., 2003; Simonoff, 1996). The principle of parsimony says that we should minimize this dimension (the trace of the matrix \mathbf{S}) as much as possible. But in the general case, this dimension is more to obtain, explicitely.

The estimator introduced by Nadaraya (1964) and Watson (1964), in the case of a simple non-parametric regression, is also written in this form since\widehat{m}_h(x)=\mathbf{s}_{x}^T\mathbf{y}=\sum_{i=1}^n \mathbf{s}_{x,i}y_iwhere\mathbf{s}_{x,i}=\frac{K_h(x-x_i)}{K_h(x-x_1)+\cdots+K_h(x-x_n)} where K(\cdot) is a kernel function, which assigns a value that is lower the closer x_i is to x, and h>0 is the bandwidth. The introduction of this metaparameter h is an important issue, as it should be chosen wisely. Using asymptotic developments, we can show that if X has density f, \text{biais}[\widehat{m}_h(x)]=\mathbb{E}[\widehat{m}_h(x)]-m(x)\sim {h^2}\left(\frac{C_1 }{2}m''(x)+C_2 m'(x)\frac{f'(x)}{f(x)}\right)and \displaystyle{{\text{Var}[\widehat{m}_h(x)]\sim\frac{C_3}{{nh}}\frac{\sigma(x)}{f(x)}}}for some constants that can be estimated (see Simonoff (1996) for a discussion). These two functions evolve inversely with h, as shown in Figure 1 (where the metaparameter on the x-axis is here, actually, h^{-1}). Keep in ming that we will see a similar graph in the context of machine learning models.

Figure 1. Choice of meta-parameter and the Goldilocks problem: it must not be too large (otherwise there is too much variance), nor too small (otherwise there is too much bias).

The natural idea is then to try to minimize the mean square error, the MSE, defined as bias[\widehat{m}_h (x)]^2+Var[\widehat{m}_h (x)], and them integrate over x, which gives an optimal value for h of the form h^\star=O(n^{-1/5}) , and reminds us of Silverman’s rule – see Silverman (1986). In larger dimensions, for continuous \mathbf{x} variables, a multivariate kernel with matrix bandwidth \mathbf{H} can be used, and \mathbb{E}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim m(\mathbf{x})+\frac{C_1}{2}\text{trace}\big(\mathbf{H}^Tm''(\mathbf{x})\mathbf{H}\big)+C_2\frac{m'(\boldsymbol{x})^T\mathbf{H}\mathbf{H}^T \nabla f(\mathbf{x})}{f(\mathbf{x})}while\text{Var}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim\frac{C_3}{n~\text{det}(\mathbf{H})}\frac{\sigma(\mathbf{x})}{f(\mathbf{x})}
If \mathbf{H} is a diagonal matrix, with the same term h  on the diagonal, then h^\star=O(n^{-1/(4+dim(\mathbf{x}))}. However, in practice, there will be more interest in the integrated version of the quadratic error, MISE(\widehat{m}_{h})=\mathbb{E}[MSE(\widehat{m}_{h}(X))]=\int MSE(\widehat{m}_{h}(x))dF(x)and we can prove that MISE[\widehat{m}_h]\sim \overbrace{\frac{h^4}{4}\left(\int x^2k(x)dx\right)^2\int\big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]^2dx}^{\text{bias}^2} +\overbrace{\frac{\sigma^2}{nh}\int k^2(x)dx \cdot\int\frac{dx}{f(x)}}^{\text{variance}}as n→∞ and nh→∞. Here we find an asymptotic relationship that again recalls Silverman’s (1986) order of magnitude, h^\star =n^{-\frac{1}{5}}\left(\frac{C_1\int \frac{dx}{f(x)}}{C_2\int \big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]dx}\right)^{\frac{1}{5}}The main problem here, in practice, is that many of the terms in the expression above are unknown. Automatic learning offers computational techniques, when the econometrician used to searching for asymptotic (mathematical) properties.

To be continued (references mentioned above are online here)…

Manipulation matricielle de données

Pour travailler un peu la manipulation de matrices, dans le cadre du cours d’analyse de données, quelques lignes de code.

La matrice de données est ici

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",header=TRUE)
> M=as.matrix(BD[,2:10])
> colnames(M)=names(BD)[2:10]
> rownames(M)=BD$dep
> head(M[,1:5])
    infra  vols  eco crim  vma
D1  44.11 29.51 4.14 2.79 0.08
D2  45.97 26.29 3.97 4.54 0.07
D3  38.83 21.90 4.61 3.33 0.03
D4  49.68 30.11 5.28 4.54 0.07
D5  47.67 27.19 4.94 3.91 0.05
D6 109.21 77.71 6.87 5.25 0.24

La décomposition spectrale de

> round(t(M)%*%M)
       infra   vols   eco  crim vma  vvi  camb  
infra 308599 191433 29179 21290 642 7042 30951 
vols  191433 119961 17859 12967 404 4427 19335 
eco    29179  17859  2946  2042  59  648  2942  
crim   21290  12967  2042  1564  42  448  2129  
vma      642    404    59    42   2   19    62   
vvi     7042   4427   648   448  19  237   665  
camb   30951  19335  2942  2129  62  665  3200  
roul   55363  34745  5089  3792 113 1217  5612 
auto   24841  15747  2214  1676  53  545  2520

(j’ai viré les deux dernières colonnes pour des raisons de présentation) donne ici

> eigen(t(M)%*%M)
[1] 4.471766e+05 1.239465e+03 3.462936e+02 1.764038e+02 1.303316e+02 6.623592e+01 2.624739e+01 1.761801e+01 1.564292e-01

(etc). Pour rappel, les liens entre la diagonalisation et l’analyse spectrale se voit avec l’écriture

> P=eigen(t(M)%*%M)$vectors
> P%*%diag(eigen(t(M)%*%M)$values)%*%t(P)

Maintenant, deux remarques. Tout d’abord le lien entre la décomposition spectrale de cette matrice, et la décomposition en valeurs singulières

> sqrt(eigen(t(M)%*%M)$values)
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
> svd(M)
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
             [,1]          [,2]          [,3]          [,4]         [,5]          [,6]         [,7]         [,8]          [,9]
 [1,] -0.08156446 -0.0632514738 -2.596265e-02  0.0960587227 -0.064759395 -0.0430947134 -0.024719699 -0.074347281  0.0550356487
 [2,] -0.08087892  0.0541859480 -5.777492e-02 -0.0628249751 -0.099084701 -0.1311683858  0.083484434  0.011872601 -0.0195034904

(etc) et la décomspotion spectrale de l’autre produit matriciel

> sqrt(eigen(M%*%t(M))$values)
 [1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 1.141629e+01 8.138545e+00 5.123221e+00 4.197381e+00 3.955113e-01 7.546505e-06 6.656212e-06 5.285131e-06 4.270035e-06


Maintenant, pour comprendre un peu mieux la recherche de composantes pricinpales, considérons les deux variables

> sM=M[,c(1,3)]
> plot(sM)


On a tout intérêt à centrer et réduire les variables (ou changer la métrique, comme vu en cours)

> sMcr=sM
> for(j in 1:2) sMcr[,j]=(sMcr[,j]-mean(sMcr[,j]))/sd(sMcr[,j])
> plot(sMcr)

Avant de faire des projections sur un axe, introduisons deux fonctions

> prod_scal=function(x,u){
+   ps=rep(NA,nrow(x))
+   for(i in 1:nrow(x)) ps[i]=sum(x[i,]*u)
+   return(ps)
+ }

> proj=function(x,u){
+   px=x
+   for(j in 1:length(u)){
+     px[,j]=prod_scal(x,u)/sqrt(sum(u^2))*u[j]  
+   }
+   return(px)
+ }

Par exemple, si on projette sur l’axe des abscisses,

> points(proj(sMcr,c(1,0)),col="blue")

On peut ensuite chercher la direction de l’axe qui nous donnera le nuage de points qui aura la plus grande inertie

> inertie=function(x) sum(x^2)
> Theta=seq(0,3.141592657,length=101)
> V=unlist(lapply(Theta,function(theta) inertie(proj(sMcr,c(cos(theta),sin(theta))))))
> plot(Theta,V,type='l')

> (angle=optim(0,function(theta) -inertie(proj(sMcr,c(cos(theta),
[1] 0.7853516

Visuellement, on obtient

> plot(sMcr)
> psMcr=proj(sMcr,c(cos(angle),sin(angle)))
> points(psMcr,col="blue")


On notera que l’axe qui donne l’inertie maximale est lié aux vecteurs propres de la décomposition spectrale (celui associé à la plus grande valeur propre)

> c(cos(angle),sin(angle))
[1] 0.7071397 0.7070738
> eigen(t(sMcr)%*%sMcr)
[1] 159.64663  26.35337
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

On continuera demain, à manipuler des matrices de données, de faire des projections, avant de commencer l’analyse en composantes principales, à proprement parler, la semaine prochaine

Visualising a Classification in High Dimension

So far, when discussing classification, we’ve been playing on my toy-dataset (actually, I should no claim it’s mine, it is inspired by the one used in the introduction of Boosting, by Robert Schapire and Yoav Freund). But in ral life, there are more observations, and more explanatory variables.With more than two explanatory variables, it starts to be more complicated to visualise. For instance, consider


where we have observations from people in E.R., for infarctus, and we want to understand who did survive, to get a predictive model. But before running some classifier, let us visualise our data. Since we have seven explanatory variables and our class (survival or death), we can go for a PCA.

library(FactoMineR) # ACP (sur les var continues)

To add the death/survival variable, treat it as numerical 0/1 variable (at least to get a direction)


The nice thing is that we see here where variables are colinear with that one. It is also possible to visualise individuals, and classes, too

plot(acp, habillage = 8,col.hab=c("red","blue"))

Continue reading Visualising a Classification in High Dimension

We are winter

Non, pas de billet pour annoncer la reprise de Game of Thrones dans les jours à venir. Mais un billet rapide pour rappeler qu’au Canada, l’hiver, c’est un sujet important, et sérieux.

C’est ce thème qu’avait repris le Comité Olympique Canadien avant et pendant les Jeux.

Et lors de mon excursion au Mexique la semaine passée (où nous avons retrouvé des collègues chiliens d’ailleurs), je me suis rendu compte que les Jeux Olympiques d’hiver n’intéressaient que peu de monde, en dehors des Canadiens, et des nations qui subissent l’hiver, comme les Russes, ou les Norvégiens (entre autres). Un autre point comment de toutes ces nations est d’avoir un intérêt très fort pour un endroit important d’un point de vue géopolitique, le pôle nord. Et on m’a demandé d’intervenir dans une conférence de Thierry Garcin, et Philippe Reka. J’ai toujours un plaisir fou à discuter avec Philippe, donc c’est avec plaisir que je viendrais, mais j’avoue (une nouvelle fois) avoir du mal à trouver ma place (oui, je suis un peu le sparadrap du capitaine Haddock)

Il y a tout juste un an, Philippe publiait dans le Monde Diplo un article inspirant sur la Cartographie Radicale,

Depuis des décennies, la cartographie traditionnelle revendique le statut de science exacte s’appuyant sur des données fiables. Elle se targue de fournir une image neutre et fidèle de la réalité. Mais une telle approche fait l’impasse sur l’utilisation politique et sociale de la carte, et sur son rôle tant de propagande que de contestation. Depuis le début des années 2000, et dans le désordre, émerge une pratique cartographique qui se dit « radicale » (on parle aussi de « cartographie critique » ou de « contre-cartographie »), riche combinaison revendiquée d’art, de sciences, de géographie, de politique et de militantisme social. Les cartographes orthodoxes ont vu naître avec beaucoup de méfiance les premiers projets exploratoires qui bousculaient, voire cassaient, les codes. « On peut “crobarder” des schémas réducteurs donnant l’illusion de la synthèse, et en faire des figures de style prétentieuses, écrivait un cartographe professionnel en écho à l’une des toutes premières expositions de cette mouvance. La carte géographique, ce fantastique outil de savoir et d’appréhension de l’espace et du monde, se voit ravalée au rang de simple outil de prestidigitation. »

C’est peu dire que le symposium « Art and Cartography — Cartography and Art », organisé par l’Académie des beaux-arts et l’Université technique de Vienne (Autriche) en février 2008, fut le théâtre d’intenses discussions, où s’exprimèrent des positions inconciliables. La cartographie engagée a trouvé la voie d’une forme de contestation au service de la justice sociale, dénonçant des pratiques politiques et économiques douteuses. Elle est un exercice libre de déconstruction de l’espace et des phénomènes sociaux, pour lequel les protagonistes se permettent de pervertir les conventions les plus classiques. Les premiers à s’en emparer furent des artistes, mais aussi des architectes, des urbanistes et des militants, rejoints plus tard par des géographes. Les projets, individuels ou collectifs, se sont surtout multipliés en Argentine, aux Etats-Unis, en Suisse, en Allemagne et en France, couvrant un champ thématique assez large (finance, surveillance et sécurité, marketing, environnement, organisation des territoires, pour ne donner que quelques exemples), et profitant largement de la puissance des outils de la cartographie participative et des réseaux sociaux.

Déchiffrer pour dénoncer, mettre en images des processus peu visibles qui concourent à confisquer l’espace public (voire des biens publics), à compromettre les libertés individuelles, à détourner des lois, tels sont les principaux objectifs des promoteurs de cette cartographie. Dans cet esprit, informer n’est qu’un début. L’étape suivante de la démarche reste l’action pour le changement sur le terrain. Largement informelle, cette initiative se reconnaît quelques affinités avec les mouvements de résistance tels que celui des « indignés » ou Occupy. Ainsi, le collectif d’artistes argentins Grupo de Arte Callejero ne s’est pas contenté de dénoncer les bourreaux de l’ancienne dictature en cartographiant les lieux de leur paisible retraite : il a aussi organisé des itinéraires « touristiques » qui passaient par leurs domiciles et mené campagne sur des panneaux publicitaires… Le projet « Duty Free Shop », en ce qu’il se propose de déconstruire, rendre lisible et donc fragiliser un processus d’aliénation des citoyens dans un espace public, s’inscrit de plain-pied dans ce mouvement pétillant.

En tant qu’amateur (plus qu’expert) de visualisation de données, je rejoins complètement Philippe dans sa lecture critique (je ne vois les cartes que comme un support graphique parmi beaucoup d’autres). Philippe a ainsi proposé une très jolie visualisation de l’arctique.

Et comme souvent, Philippe a raison. Pour comprendre l’arctique, il faut avoir une visualisation claire. Car comme beaucoup de monde, j’ai été éduqué avec des cartes qui ressemblaient à ça

Autrement dit, pour moi, l’arctique et le pôle nord ont longtemps été des concepts flous, un peu comme l’horizon… c’est au nord, tout là haut. Par contre, j’ai bien conscience que l’antarctique existe, par exemple, mais les cartes usuelles ne nous permettent de prendre conscience de sa superficie (un peu moins de deux fois la taille de l’Australie). D’ailleurs, si on va sur https://www.google.ca/maps/place/North+Pole/ on peut voir (en butant contre la partie supérieure de la carte) qu’on ne voit rien, si on cherche à visualiser ce pôle nord (même s’il est sous la glace)

Comme le notait ma fille, le pôle nord est une enjeu stratégique important car c’est là que se trouve le village du père Noël (c’est manifestement la croyance répandue en Amérique). Et si le père Noël n’existe pas (c’est ce que pense ma fille), peut être que le pôle nord n’existe pas non plus ? Le premier défi est donc d’avoir une visualisation aussi pertinente que possible du pôle. Je vais donc revenir un peu sur les problématique de projection. Car les Soviétiques ne se trompaient pas lorsqu’ils représentèrent l’arctique avec un point de vue qui donnait l’impression que l’arctique était russe,

Pour appréhender les projecteurs, deux ouvrages fabuleux expliquent cela très clairement. Flattening the Earth – Two Thousand Years of Map Projections de John Parr Snyder, mais aussi (et surtout) How to Lie with Maps de Mark Monmonier. Mark Monmonier a également publié un livre passionnant sur le sujet (en particulier le chapitre 6), intitulé Air Apparent – How Meteorologists Learned to Map, Predict, and Dramatize Weather

Pour visualiser le pôle nord, on peut utiliser quelques fonctions simples (basées sur un travail en cours avec @3wen), pour visualiser différentes projections. On va ici mettre l’accent sur trois pays importantes quand on parle du pôle nord, le Canada, la Russie, et le Groenland (qui est techniquement un pays consultatif), et sur trois deux régions, une située à moins de 1500km du pôle nord, et l’autre à 3000km du pôle (soit un peu en dessous du cercle arctique)

Si on veut changer de point de vue, on peut visualiser l’animation suivante

D’ailleurs, si on compare la proportion du disque occupé par chacun des pays (j’ai pris ici chacune des huit pays de l’arctique), on obtient, avec en abscisse la distance au pôle, et en ordonnée, la proportion du disque occupée par chacun des pays, avec le Canada (rouge), la Russie (bleu) et le Groenland (vert)

Techniquement, le pôle pose de réels soucis car la plupart des données sont représentées par leur latitude, et leur longitude. Les pôles sont les seuls points – dans cette représentation – à avoir une infinité de coordonnées !

File:Latitude and Longitude of the Earth.svg

La latitude est de 90°, mais toutes les longitudes sont valides. On a ainsi des soucis y compris sur les cartes construites à partir de données satellites, qui pourtant pourraient balayer uniformément la sphère.

Le pôle nord peut ainsi compter beaucoup de discontinuités, uniquement à cause de la représentation des données par la latitude et la longitude : en quelques kilomètres, on passe d’une zone rouge à une zone blanche, ce qui pose des soucis de lissage de la carte (les données n’étant manifestement pas considérées ici comme circulaire pour le lissage).

Maintenant, plus qu’un point sur la carte, on peut voir le pôle nord comme une région, et non pas comme un point sur la sphère, défini très arbitrairement. Bon, j’exagère un peu ici, car le pôle nord est défini par l’axe de rotation de la terre,

qui va induire une espère de symétrie, par exemple sur le climat. Mais le second problème, selon moi, qui va alors apparaître est que le pôle nord vu comme un territoire (un territoire de glace) est alors un territoire qui va changer avec les saisons.

Par exemple, en septembre dernier, la glace occupait la région suivante, ce qui montre que peu de pays peuvent effectivement rejoindre le pôle nord ‘à pieds’ quel que soit le moment de l’année.

On continue de travailler avec @3wen sur le sujet, mais le sujet et vaste, ambitieux… et complexe ! à suivre donc…

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


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


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


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,


qui va impliquer l’inégalité triangulaire


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,


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


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


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,


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


c’est à dire que


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


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


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,


Techniquement, je pense que


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


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


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


ce qui se traduit par la relation bien connue


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


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


et pour les séries temporelles stationnaires,


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


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


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


et la corrélation comme


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


Ce qui se traduit par


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




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


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


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


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.


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)

Equidistant points on a map

This morning, I had a comment on a recent post, regarding a graph I did upload on the blog, which was extracted from a paper now online (see http://hal.archives-ouvertes.fr/hal-00871883). Jo (from KUL, I guess I can share that piece of information) asked me

I was wondering whether you would want to share the R code for plotting figures 1 and 14? W.r.t. the former, the figure-in-figure is a nice touch; as to the latter, I am curious to know how you translated distance in km to the size parameters of the graph (par(“usr”)) for plotting the corresponding concentric circles (and the arrow indicating the radius) on top of your map.

At first, I thought I made a mistake in my plot. I mean, each time I have a question, I start to be suspicious, and I start to wonder if what I did was valid, or not. Here was the graph

Let’s make it clear: I do not draw circles here. So yes, I believe that what I did is valid. What I did is simple. First, I get the background map,

map("world",xlim=c(130,150),ylim=c(25,45),fill=TRUE,col="light green")

Then, I need some function to compute distance from coordinates. The functions I use are

deg2rad = function(deg) return(deg*pi/180)
DISTANCEDEG = function(long1, lat1, long2, lat2) {
R=6371; d=acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R

The center here will be Tokyo (東京),


The idea now is simple: I generate a grid (here 501×501)


I compute the distance from all those points to Tokyo, and check is the distance is larger or smaller than a given value,


If the distance is smaller than 500km, then I put a blue dot on the graph,

points(VtX[D1],VtY[D1],pch=19,cex=.2,col="light blue")

Then I use the same procedure for 250km (obviously, it is more convenient to start from larger and to go to smaller distances)

points(VtX[D],VtY[D],pch=19,cex=.2,col="light yellow")

Then, I did draw an arrow to ilustrate the largest distance

text(((VtX[D1])[k]+X)/2,Y+.35,"500 km")

And now, I have the graph.

Now, the point is that it should depend on the kind of projection we use, right? So here is a function that can be used for different kind of projections (some slight changes are necessary, since the map is now centered on some point, and we cannot use standard coordinates)

mapjapan = function(pr="conic",pm=45){
map("world","japan",fill=TRUE,col="light green",projection=pr, par=pm)
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
map("world","japan",projection=pr, par=pm,add=TRUE)

The default function here produces a map based on a conic projection,


But we can also use a Bonne projection (a pseudo-conic one, named after Rigobert Bonne)


or Lagrange projection,


and (as a last one), Albers projections,


Of course, much more projections are possible !

We do not see much here, right ? So let us play with a larger country to visualize something. Like Canada. And the distance to, say, Winnipeg.

The first thing to do is to define the coordinates of Winnipeg,


Then, we slightly change our function

mapcanada = function(pr="conic",pm=45){
map("world","canada",fill=TRUE,col="light green",projection=pr, par=pm)
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
map("world","canada",projection=pr, par=pm,add=TRUE)

Now, we can have some fun





Fun, isn’t it ? Changing the projection will change the shape of equidistant curves.

Régression, variables explicatives et géométrie

La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections

Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif, deux ensembles (a priori disjoints) de variables explicatives,


Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont


de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par


La première partie correspond à la régression de http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW13.gif, mais il reste un second terme dès lors que http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif ne sont pas orthogonales. Notons http://freakonometrics.hypotheses.org/files/2016/05/FW04.gifla matrice de projection (orthogonale) sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif


Alors on peut écrire


En posant

http://freakonometrics.hypotheses.org/files/2016/05/FW07.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW08.gif,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif , i.e. en transformant les variables),


Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents

Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples


Le graphique ci-dessous au cas non orthogonal,


Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,

> chicago=read.table(
+ "http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
> tail(base)
Y    X1 X2     X3
42  4.8 0.152 19 13.323
43 10.4 0.408 25 12.960
44 15.6 0.578 28 11.260
45  7.0 0.114  3 10.080
46  7.1 0.492 23 11.428
47  4.9 0.466 27 13.731
> regression=lm(Y~X1+X2+X3)
> summary(regression)

lm(formula = Y ~ X1 + X2 + X3)

Min     1Q Median     3Q    Max
-9.737 -4.565 -1.479  3.751 16.079

Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X1          -0.62764    5.28130  -0.119 0.905953
X2           0.22378    0.06161   3.632 0.000744 ***
X3          -1.55059    0.38195  -4.060 0.000204 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.527 on 43 degrees of freedom
Multiple R-squared: 0.4417,	Adjusted R-squared: 0.4027
F-statistic: 11.34 on 3 and 43 DF,  p-value: 1.314e-05

> n=length(Y)
> X=matrix(c(rep(1,n),X1,X2,X3),n,4)
> X[1:5,]
[,1]  [,2] [,3]   [,4]
[1,]    1 0.604   29 11.744
[2,]    1 0.765   44  9.323
[3,]    1 0.735   36  9.948
[4,]    1 0.669   37 10.656
[5,]    1 0.814   53  9.730

Ensuite, on va projeter seulement sur les deux premières variables (et la constante)

> FWX1=X[,1:3]
> regression12=lm(Y~X1+X2)
> solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y
[1,]  0.08069764
[2,] 11.56913900
[3,]  0.15108490
> summary(regression12)$coefficients
Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  0.08069764 3.49182154 0.02311047 0.98166664
X1          11.56913900 5.05014782 2.29085156 0.02681947
X2           0.15108490 0.06854408 2.20420043 0.03278652
> FWX2=X[,4]
> H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1)
> M1=diag(rep(1,n))-H1
> FWX2s=M1%*%FWX2
> FWYs =M1%*%Y
> (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs)
[1,] -1.550594
> summary(regression)$coefficients[4]
[1] -1.550594
> (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y-
+        solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2)
[1,] 22.0752495
[2,] -0.6276442
[3,]  0.2237765
> summary(regression)$coefficients[1:3]
[1] 22.0752495 -0.6276442  0.2237765

On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.