Tag Archives: distance

D’un abus de langage…

Dans mon cours d’analyse des données, je commençais par présenter les diverses distances et normes, en particulier, celles qui dérivent d’un produit scalaire, comme la norme Euclidienne \|\boldsymbol{X}\|_{\mathbb{I}}^2=\boldsymbol{X}^\top\boldsymbol{X}, ou plus intéressante, la norme que l’on peut construire à partir d’une matrice \boldsymbol{M} symétrique définie positive, \|\boldsymbol{X}\|_{\boldsymbol{M}}^2=\boldsymbol{X}^\top\boldsymbol{M}\boldsymbol{X}. Un exemple bien connu étant la norme de Mahalanobis, ou la distance associée, où on prend pour \boldsymbol{M}  la matrice \boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}=\text{Var}[\boldsymbol{X}]. Sauf que, par abus de langage, j’appelle aussi distance de Mahalanobis le distance obtenue en prenant pour matrice \boldsymbol{M} la matrice diagonale \boldsymbol{\Delta}_{\boldsymbol{\Sigma}}^{-1} définie par \boldsymbol{\Delta}_{\boldsymbol{\Sigma},ii}=\boldsymbol{\Sigma}_{ii}. En fait, c’est cette matrice que l’on utilise principale en analyse de données, par exemple en analyse des correspondances principales… Prenons un petit exemple, avec quelques données centrées

library(FactoMineR)
library(clusterGeneration)
library(mnormt)
n=15
S=genPositiveDefMat("unifcorrmat",dim=4)
X=round(rmnorm(n,varcov=S$Sigma)+8,2)
for(j in 1:ncol(X)) X[,j]=X[,j]-mean(X[,j])
rownames(X)=1:n
colnames(X)=LETTERS[1:4]

et seulement 4 variables. Par défaut, quand on fait une analyse en composantes principales, on va réduire les données

acp = PCA(X,scale.unit = TRUE,graph = FALSE)
fviz_pca_var(acp, col.var = "cos2")

Le fait de centrer et réduire les données fait que l’on travaille avec la matrice de corrélation (\boldsymbol{X}^\top\boldsymbol{X}, à un facteur n, ou n-1, près). D’où le cercle de corrélation qui apparaît sur la figure. Et abusivement, quand je présente l’ACP, je dis qu’on a deux manière de voir : (1) utiliser le produit scalaire et la distance usuelle (Euclidienne) sur les données centrées et réduites, (2) utiliser la distance de Mahalanobis sur les données centrées. Mon abus de langage tient du fait qu’on utilise la version avec la matrice diagonale \boldsymbol{\Delta}_{\boldsymbol{\Sigma}}^{-1} , et pas la vraie distance de Mahalanobis… En discutant avec Gabriel après mon cours de ce matin, je me suis rendu compte de la confusion possible (même si je me faisais la réflexion lors de mes enregistrements de capsules, les jours passés). Mais en réalité, faire une analyse en composantes principales avec la vraie distance de Mahalanobis n’aurait aucun sens.

Pour s’en convaincre, utilisons les données centrées et réduites,

S = var(X)
Sd = diag(diag(S))
X1 = X%*%solve(chol(Sd))
colnames(X1)=colnames(X)
var(X1)
A B C D
A 1.0000000 0.10113165 -0.88021475 0.4733021
B 0.1011316 1.00000000 0.07885109 -0.7400999
C -0.8802147 0.07885109 1.00000000 -0.5743660
D 0.4733021 -0.74009993 -0.57436600 1.0000000

x

acp1 = PCA(X1,scale.unit = FALSE,graph = FALSE)
fviz_pca_var(acp1, col.var = "cos2")

On n’a plus ici, le cercle des corrélations car la fonction ne sait pas que mes données ont été centrées et réduite au préalable… mais on a la même sortie qu’auparavant à part ça…

Si maintenant on utilise la norme de Mahalanobis, on peut là aussi tranformer les données, en notant que \boldsymbol{X}^\top\boldsymbol{M}\boldsymbol{X} peut aussi s’écrire \boldsymbol{X}^\top\boldsymbol{M}^{1/2}\boldsymbol{M}^{1/2}\boldsymbol{X} soit encore la distance euclidienne standard \tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}} sur les données transformées \tilde{\boldsymbol{X}}={\boldsymbol{X}}{\boldsymbol{M}^{1/2}},

S = var(X)
X2 = X%*%solve(chol(S))
colnames(X2)=colnames(X)
> var(X2)
              A             B             C             D
A  1.000000e+00 -3.582611e-17  5.307644e-16 -1.004784e-16
B -3.582611e-17  1.000000e+00 -3.021333e-16  4.164692e-17
C  5.307644e-16 -3.021333e-16  1.000000e+00 -3.649967e-16
D -1.004784e-16  4.164692e-17 -3.649967e-16  1.000000e+00

ce qui donne des données… orthogonales (pour ne pas dire indépendantes, on retrouve ici le processus de whitening, utilisé en ICA, analyse des composantes indépendantes. Et comme on a rendu les données indépendantes… l’analyse des composantes principales n’aura plus grand chose à raconter…

acp2 = PCA(X2,scale.unit = FALSE,graph = FALSE)
fviz_pca_var(acp2, col.var = "cos2")

Chaque axe explique ici pour 25% de l’inertie… bref, on ne va jamais faire une ACP avec la distance de Mahalanobis (la vraie). Elle ne va pas normaliser les variables, mais les orthonormaliser. Et pourtant, par abus de langage, je vais continuer à utiliser ce terme… Désolé…

INF7100, “rappels” de maths

En parallèle des trois parties de mon intervention sur la science des données, dans le cadre du cours INF7100, j’ai mis en ligne une série des capsules de “rappels” de maths.

  • 401: Concepts généraux (logarithme, exponentielle, dérivée, intégrale, approximation) pdf video (38:57)
  • 411: Optimisation (condition du premier ordre, descente de gradient, Lagrangien) pdf video (31:33)
  • 412: Ensembles (intersection, unions, fusion de bases, convexité, centroïdes) pdf video (10:13)
  • 421: Vecteurs et Matrices (orthogonalité, distance, norme, tenseurs) pdf video (43:20)
  • 422: Algèbre Linéaire (déterminant, valeurs propres, vecteurs propres, formes quadratiques, approximations) pdf
  • 423: Coordonnées (données circulaires, coordonnées polaires, projection, spatial) pdf
  • 431: Loi Gaussienne (loi normale, cas multivariée, limite et approximation, loi des grands nombres) pdf video (17:37)
  • 432: Lois Puissance (changement d’échelle, loi puissance, Zipf, Benford) pdf
  • 433: Loi de Durée (loi géométrique, loi exponentielle, loi de Poisson, loi des petits nombres) pdf
  • 434: Modèle de Markov (données séquentielles) pdf

Where hiding if you don’t want to get wet ?

Following the previous post, two additional remarks. Following a comment by@cosi, I have investigated quickly a binomial fit to the distribution of the number of people not getting wet, with a fixed number of players on the field. It looks like it should be a binomial distribution with a fixed probability (2/3) and with size parameter affine in the number of players. @guigui suggested some connexion with with “birds on a wire” problem (see e.g. http://www.cut-the-knot.org/)

n=p=rep(NA,20)
for(i in 1:20){
NSim=10000
N=Vectorize(NOTWET)(n=rep(3+2*i,NSim))
n[i]=mean(N)/(1-var(N)/mean(N))
p[i]=1-var(N)/mean(N)
}
plot(seq(5,43,by=2),n,col="red",type="b")

for the implied size parameter above, and below the implied probability parameter.

plot(seq(5,43,by=2),p,col="blue",type="b")

(as functions of the number of players). I’d be glad to get more details on that 2/3 probability.

Now, let us investigate another question sent by email: “Where should you hide if you don’t want to get wet ?” A first idea could be the following: given that some players are already on the field, where should I go if I do not want to get wet ? Below are some simulations for 7 or 25 players (already on the field). The red area is the area so that I will become someone’s target (perhaps even the target of two players…). The green area is the safe zone.

(with 7 players above, and 25 below)

It looks like, on the border, it might be safer than in the middle of the field. But we have to confirm that intuition… or at least see if that intuition is valid.

Based on what was done the other day, it is possible to look where people that got wet were located (instead of counting dry players as done in the previous function). So here, we simply look where non wet players were standing

NOTWET=function(n,p){
x=runif(n)
y=runif(n)
(d=as.matrix(dist(cbind(x,y), method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
whonotwet=( (1:n) %notin% names(table(dmin)) )
#plot(x[-whonotwet],y[-whonotwet],pch=19,col="blue",type="p")
#points(x[whonotwet],y[whonotwet],pch=19,col="red")
M=matrix(NA,p,p);u=seq(0,1,by=1/p)
for(i in 1:p){
for(j in 1:p){
M[i,j]=sum((x[whonotwet]>=u[i])&(x[whonotwet]<u[i+1])&
(y[whonotwet]>=u[j])&(y[whonotwet]<u[j+1]))
}}
return(M)}

based on function

"%notin%" <- function(x, y) x[!x %in% y]

On a given grid, we count people playing the game that ended dry (with might avoid boundary bias on nonparametric smooth estimator of distribution, as we’ll see later on). For instance with 11 players,

M11=matrix(0,25,25);
for(s in 1:100000){
M11=M11+NOTWET(11,25)
}

Then we can plot the distribution, on the field,

COL=rev(heat.colors(101)); p=25
u=seq(0,1,by=1/p)
plot(0:1,0:1,col="white",xlab="",ylab="")
for(i in 1:p){
for(j in 1:p){
polygon(c(u[i],u[i],u[i+1],u[i+1]),
 c(u[j],u[j+1],u[j+1],u[j]),border=NA,
col=COL[trunc(M11[i,j])/max(M11)*100+1])
}}

Red means a lot of non-wet people (i.e safer zones). Graphs below are with 7 and 11 players respectively (from the left to the right)

with the following distribution on the diagonal: corners are almost 4 times safer than the middle of the field, with 7 players,

Below are plotted distributions of locations of non-wet players when the total number of players was either 25 (on the left) and 101 (on the right)

with again on the diagonal

Hence, the border is rather safe, but next to the border, it is no safe any longer: is someone is standing right on the border, he will probably shoot at you: there is no one behind him ! This explains the stange behavior on the borders (and corners, thanks JP for the intuitive explanation).
But would it be completely different with a field shaped as a disk ?

using the previous technique of working on a fixed grid (or correcting for boundary bias, since the disk might cover only a fraction of the grid-square), or keeping coordinates of non-wet players, and using standard kernel-based estimator of the distribution (the light yellow circle outside the disk is simply due to bias of the kernel estimator on the border)

NOTWET=function(n){
x=(runif(n*20)*2-1)*1
y=(runif(n*20)*2-1)*1
I=which((x^2+y^2<1))
x=x[I];y=y[I]
x=x[1:n];y=y[1:n]
(d=as.matrix(dist(cbind(x,y),
method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
whonotwet=( (1:n) %notin% names(table(dmin)) )
return(cbind(x[whonotwet],y[whonotwet]))
}

M=t(c(0,0))
for(s in 1:10000){
M=rbind(M11,NOTWET(25))
}
M=M[-1,]

library(ks)
HP=matrix(c(.001,0,0,.001),2,2)
K=kde(x=M11, H=HP)
image(K$eval.points[[1]],K$eval.points[[2]],K$estimate2,
col=rev(heat.colors(101)),xlim=c(-1,1),ylim=c(-1,1))

 

And note that the distribution of the number of players ending the game dry is the same, for a square field, or a disk,

NOTWET2=function(n){
x=(runif(n*20)*2-1)*1
y=(runif(n*20)*2-1)*1
I=which((x^2+y^2<1))
x=x[I];y=y[I]
x=x[1:n];y=y[1:n]
(d=as.matrix(dist(cbind(x,y), 
method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
notwet=n-length(table(dmin))
return(notwet)}

NOTWET=function(n){
x=runif(n)
y=runif(n)
(d=as.matrix(dist(cbind(x,y), 
method = "euclidean",upper=TRUE)))
diag(d)=999999
dmin=apply(d,2,which.min)
notwet=n-length(table(dmin))
return(notwet)}

NSim=100000
Nsquare=Vectorize(NOTWET)(n=rep(25,NSim))
Ndisk=Vectorize(NOTWET2)(n=rep(25,NSim))
Tsq=table(Nsquare)
Tdk=table(Ndisk)
plot(as.numeric(names(Tsq)),Tsq/NSim,
type="b",col="red")
lines(as.numeric(names(Tdk)),Tdk/NSim,
type="b",pch=4,col="blue")


But so far, it was still simple… I wonder what it might become if we consider a non-convex place, with walls, where player might hide…. Next time, a post on indoor paint-ball !

Vitesse et comportement dangereux au volant (partie 3)

Après plusieurs billets sur les accidents de la route (ici et ), parlons un peu de vitesse des véhicules. J’ai pu récupérer des données via l’inrets, d’observations près de Saclay en région parisienne, sur une 2×2 limitée à 110 km/h.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/.n118-saclay2_m.jpg

J’ai ainsi le passage de plusieurs millions de véhicules, sur quelques jours, incluant tous les véhicules. Mais si j’ai signé des papiers quant à la confidentialité des données, je dois pouvoir mettre des dessins construits à partir de ces données…

  • A quelle vitesse roulent les conducteurs parisiens ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, qui poussent les gens à ralentir. La courbe en noir est la vitesse moyenne des véhicules qui passent sur ce tronçon, et la courbe en bleu, la vitesse des 5% des véhicules les plus rapides,

Ces pics sont légèrement atténués le week end,

avec ici un trafic très dense tout l’après midi, entre 15 et 20 heures.

  • Quelle proportion de véhicules dépasse la limite autorisée ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, avec en gras la proportion des véhicules qui dépasse 110 km/h, et en trait fin, la proportion qui dépasse 115 km/h (soit 5% de plus),

avec des choses finalement assez proches le week end,

Autrement dit, aux heures de pointe, peu de monde dépasse la vitesse autorisée, mais vers midi, ou surtout la nuit, environ un véhicule sur quatre dépasse la vitesse autorisée (mais seulement un sur six la dépasse de plus de 5%).

  • Vitesse et respect des distances de sécurité

La base est très complète, puisque j’ai tous les véhicules qui sont passés près de Saclay, avec l’heure d’observation, la vitesse, ce qui permet d’obtenir la distance entre deux véhicules qui se suivent (en supposant que le premier véhicule maintienne sa vitesse constante).
Rappelons qu’il existe un lien théorique entre la vitesse d’un véhicule et la distance de freinage (expliqué ici ou ). Pour faire simple, la distance de freinage https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein02.png est la somme d’une distance parcourue avant de réagir https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein04.png, et d’une distance de freinage proprement dite https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein03.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein01.png

Si on suppose que la personne met une seconde à réagir, la distance de réaction est simplement la vitesse https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein05.png(exprimée en m/h), i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein06.png

Quant au freinage à proprement parler, rappelons que les équations qui caractérisent le mouvement d’un objet sont de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein10.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein11.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein12.png

Aussi, en mettant la deuxième au carré, et en combinant avec la troisième, on a

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein13.png

avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein15.png (on cherche à arrêter le véhicule) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein14.png, aussi, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein16.png, on obtient quelque chose qui peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein17.png

à condition d’avoir des vitesses en km/h. Aussi en agglomération, il faut 23 m pour s’arrêter. Et sur notre route à Saclay, pour une personne qui roule à 110 km/h, il faut compter 110 m (ce qui correspond au message véhiculé par la sécurité routière).
Sur nos données, on peut représenter en rouge les individus qui ne respectent pas la distance de sécurité.

La courbe bleue est la distance moyenne des véhicule en fonction de la distance qui les sépare du véhicule qui les précède,
Aussi, en faisant une régression logistique (lissée), on obtient que 60% des conducteurs qui roulent entre 50 et 110 km/h ne respectent pas les distances de sécurité.