Heuristics on Correspondance Analysis

This week, in the course on non-supervised techniques for data science, we’ve been using a dataset, with candidate for the presidential elections in 2002 (per row) and newpapers (per column). In order to visualize that dataset, consider three candidates, and three newspapers

> base=read.table(
"http://freakonometrics.free.fr/election2002.txt",header=TRUE)
> sb=base[,c(2,3,4)]
> sb=sb[c(4,12,7),]
> (N=sb)
       LeFigaro Liberation LeMonde
Jospin        7         41      26
Chirac       35          9      18
Mamere        1         10       7

The first part is based on a description of rows. Consider here rows are conditional probabilities, in the set of newspapers,

> (L=N/apply(N,1,sum))
         LeFigaro Liberation   LeMonde
Jospin 0.09459459  0.5540541 0.3513514
Chirac 0.56451613  0.1451613 0.2903226
Mamere 0.05555556  0.5555556 0.3888889

The “average row” is the marginal distribution of newspapers

> (Lbar=apply(N,2,sum)/sum(N))
  LeFigaro Liberation    LeMonde 
 0.2792208  0.3896104  0.3311688

If we visualize those individuals, in the set of newpapers (in the simplexe in the newspapers space), we have

Here it is,

But actually, we will not stay in the simplexe. A PCA is considered, with weights on individuals, that take into account the importance of the different candidates, and weights for the scalar product (in order to have a distance related to the chi-square distance, and not a standard Euclidean distance)

> matL0=t(t(L)-Lbar)
> library(FactoMineR)
> acpL=PCA(matL0,scale.unit=FALSE,
+   row.w=(apply(N,1,sum)),
+   col.w=1/(apply(N,2,sum)))
> plot.PCA(acpL,choix="ind",ylim=c(-.02,.02))

 

The second part is based on a description of columns. Here Columns are conditional probabilities, in the set of candidates,

> (C=t(t(N)/apply(N,2,sum)))
         LeFigaro Liberation   LeMonde
Jospin 0.16279070  0.6833333 0.5098039
Chirac 0.81395349  0.1500000 0.3529412
Mamere 0.02325581  0.1666667 0.1372549

Here again, we can compute the “average column”

> (Cbar=apply(N,1,sum)/sum(N))
   Jospin    Chirac    Mamere 
0.4805195 0.4025974 0.1168831

In the simplex, points are

i.e.

But here again, we won’t use that simplexe. We consider a PCA, with two vectors of weights, some to take into account the weights of the newspapers, and some to get a chi-square distance

> Cbar=apply(N,1,sum)/sum(N)
> matC0=C-Cbar
> acpC=PCA(t(matC0),scale.unit=FALSE,
+          row.w=(apply(N,2,sum)),
+          col.w=1/(apply(N,1,sum)))

 

Now, we can almost overlap the two projections. Almost because we might, sometime switch right and left, top and bottom. Because if  is a (unit) eigenvector, so is . Here, for instance, we should swich them

> CA(N)

Contraintes sur les autocorrélations d’une série temporelle stationnaire

En cours de séries temporelles, ce matin, on a fini avec les contraintes sur les auto-corrélations, compte tenu du fait que la matrice d’autocorrélations doit être positive. On l’a vu, la matrice d’autocorrélations et la matrice de Toeplitz de déterminant

Avec la règle qu’on apprend tout petit pour calculer les déterminants 3×3, on a la somme des produits des termes pourpres, bleus et verts

à laquelle on retranche la somme des produits des termes verts, bleus et oranges,

Bref, on obtient assez facilement que

soit

avec de plus les corrélations dans l’intervalle ,

La région admissible pour la paire de corrélations  est alors la région bleue ci-dessous

 

Niort, Data Day

Lundi et mardi, je serais à Niort où la MAIF organise un Data Day avec plusieurs conférences (tables rondes) autour de la donnée

L’arrivée du digital créée de nouvelles opportunités, de nouveaux besoins, et de nouveaux risques. Dans ce contexte fortement évolutif et très complexe comment les personnes que nous sommes vont-ils évoluer autour de l’information de son usage et de sa gouvernance ? Quelles tendances peut-on observer ? Quelles compétences à quels horizons et pour quelles attentes ?

Analyse des Correspondances et ACP

Considérons les données suivantes, correspondant à des professions au sein d’un couple. On dispose du tableau de contingence suivant

> base=read.table(
"http://freakonometrics.free.fr/epoux.csv",
sep=";",header=TRUE)
> rownames(base)=base$Nom

Classiquement, avec ce genre de données, on est habitué à utiliser la distance du chi-deux, et les contributions au chi-deux pour voir les modalités qui correspondent fortement

> M=base[1:9,2:10]
> CT=chisq.test(M)

On peut le visualiser de la manière suivante,

> mosaicplot(t(M),main="",las=2,shade=TRUE)

avec les maris en ligne, et les épouses en colonnes. Les contribution importantes sont en bleu ou rouge foncé, les couleurs correspondant respectivement à un lien ‘positif’ (probabilité jointe plus importante que le cas indépendant) ou un lien ‘négatif’ (probabilité jointe moins importante que le cas indépendant)

ou dans l’autre sens

> mosaicplot(M,main="",las=2,shade=TRUE)

mais avec la même conclusion que précédemment: il y a de fortes valeurs bleues sur la diagonale.

Autrement dit, les couples sont relativement homogènes en terme de professions.

En analyse des correspondance, on va regarder notre tableau de contingence, en ligne ou en colonne. Par exemple, on peut définir les profils lignes, qui sont des vecteurs de probabilité \boldsymbol{L}_i=\left(\frac{p_{i1}}{p_{i\cdot}},\cdots,\frac{p_{iJ}}{p_{i\cdot}}\right)\in\mathbb{R}^J

> N=M 
> L=N/apply(N,1,sum)

En notant \boldsymbol{D}_I=\text{diag}(p_{1\cdot},\cdots,p_{I\cdot}), on peut écrire \boldsymbol{L}=\boldsymbol{D}_I^{-1}\boldsymbol{N}Le barycentre de nos vecteurs de lignes est ici\overline{\boldsymbol{L}}=\left(p_{\cdot 1},\cdots,p_{\cdot J}\right)\in\mathbb{R}^J

Là encore, en notant \boldsymbol{D}_J=\text{diag}(p_{\cdot 1},\cdots,p_{\cdot J}), on peut écrire ce barycentre de manière matricielle, \overline{\boldsymbol{L}}=\mathbb{I}_{IJ}\boldsymbol{D}_J.\boldsymbol{L}_{0,i}=\boldsymbol{L}_i-\overline{\boldsymbol{L}}

> Lbar=apply(N,2,sum)/sum(N)
> matL0=t(t(L)-Lbar)

A chaque point, on a associé la fréquence (relative) comme poids, p_{i\cdot}, ce qui correspond à utiliser la matrice \boldsymbol{D}_I. Pour mesurer la distance entre deux points, on va pondérer la distance euclidienne par l’inverse des des probabilités, p_{\cdot j}, ce qui revient à utiliser la matrice \boldsymbol{D}_J^{-1}. La distance entre deux lignes estd(\boldsymbol{L}_{i_1},\boldsymbol{L}_{i_2})=\sum_{j=1}^J\frac{1}{p_{\cdot j}}\left(\frac{p_{i_1j}}{p_{i_1\cdot}}-\frac{p_{i_2j}}{p_{i_2\cdot}}\right)^2

On va alors faire une analyse en composantes principales avec ces différents poids. D’un point de vue matriciel, on va alors l’analyse spectrale de\boldsymbol{H}=\boldsymbol{D}_J^{-1/2}\boldsymbol{L}_0^{\text{T}}\boldsymbol{D}_I\boldsymbol{L}_0\boldsymbol{D}_J^{-1/2}

En particulier, on note \boldsymbol{U}_1,\cdots,\boldsymbol{U}_k les vecteurs propres, et on défini les composantes principales\boldsymbol{C}=\boldsymbol{L}_0\boldsymbol{D}_J^{-1/2}\boldsymbol{U}

La projection sur les deux premières composantes des lignes donne ici

> library(FactoMineR)
> acpL=PCA(matL0,scale.unit=FALSE,
+          row.w=(apply(N,1,sum)),
+          col.w=1/(apply(N,2,sum)))

L’idée est de visualiser le nuage des individus correspondants aux lignes. Dans un second temps, on va faire exactement pareil, mais sur les colonnes\boldsymbol{C}_j=\left(\frac{p_{1j}}{p_{\cdot j}},\cdots,\frac{p_{jI}}{p_{\cdot j}}\right)\in\mathbb{R}^I

> C=t(t(N)/apply(N,2,sum))

Le barycentre est ici\overline{\boldsymbol{C}}=\left(p_{1\cdot},\cdots,p_{I\cdot}\right)\in\mathbb{R}^I

> Cbar=apply(N,1,sum)/sum(N)
> matC0=C-Cbar

et on peut alors faire une analyse en composante principale

> acpC=PCA(t(matC0),scale.unit=FALSE,
+          row.w=(apply(N,2,sum)),
+          col.w=1/(apply(N,1,sum)))

ce qui donne, en regardant le nuage des individus,

La magie de l’analyse des correspondances est qu’on ‘peut’ représenter les deux projections des individus sur le même plan,

> matC=acpC$ind$coord
> matL=acpL$ind$coord
> plot(matC[,1:2],col="red",xlim=c(-.002,.015))
> text(matC[,1:2],rownames(matC),col="red",pos=4)
> points(matL[,1:2],col="blue")
> text(matL[,1:2],rownames(matL),col="blue")

C’est exactement ce que fait la fonction

> afc=CA(N)

Note: en réalité, si on regarde les coordonnées des lignes par ACP on obtient

acpL$ind$coord
                    Dim.1     Dim.2     Dim.3     Dim.4     Dim.5
Agriculteur (M) -0.000056 -0.000414  0.002940  0.000150  0.000011
Artisan (M)     -0.000014  0.000152 -0.000057  0.000323  0.001279
Cadres (M)      -0.000262  0.001734 -0.000015  0.000359 -0.000219
Prof_Int (M)    -0.000194  0.000450 -0.000003 -0.000504  0.000041
Employe (M)     -0.000104 -0.000317 -0.000083 -0.000188  0.000011
Ouvrier (M)     -0.000052 -0.000764 -0.000080  0.000133 -0.000141
Retraite (M)     0.007771  0.000502 -0.000001 -0.000089 -0.000065
Inactif (M)     -0.000018 -0.000150 -0.000163  0.000751 -0.000007
Inconnu(M)       0.000261 -0.000492 -0.000554  0.000398  0.000648

alors que

CA(M, graph = FALSE)$row$coord
                 Dim 1  Dim 2  Dim 3  Dim 4  Dim 5
Agriculteur (M) -0.028 -0.209  1.481  0.076  0.005
Artisan (M)     -0.007  0.077 -0.029  0.163  0.644
Cadres (M)      -0.132  0.873 -0.007  0.181 -0.110
Prof_Int (M)    -0.098  0.227 -0.001 -0.254  0.020
Employe (M)     -0.053 -0.160 -0.042 -0.095  0.006
Ouvrier (M)     -0.026 -0.385 -0.040  0.067 -0.071
Retraite (M)     3.914  0.253 -0.001 -0.045 -0.033
Inactif (M)     -0.009 -0.075 -0.082  0.378 -0.004
Inconnu(M)       0.132 -0.248 -0.279  0.201  0.327

ce qui donne un ratio constant,

acpL$ind$coord/CA(M, graph = FALSE)$row$coord
                      Dim.1       Dim.2       Dim.3       Dim.4       Dim.5
Agriculteur (M) 0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Artisan (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Cadres (M)      0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Prof_Int (M)    0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Employe (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Ouvrier (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Retraite (M)    0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Inactif (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Inconnu(M)      0.001985182 0.001985182 0.001985182 0.001985182 0.001985182

qui correspond à 1/\sqrt{n}

1/sqrt(sum(M))
[1] 0.001985182

Les Espérances de Vie (à la naissance)

Ce soir, nous avions une discussion passionnante avec @@ et @ au sujet d’un article paru dans Le Monde

Comme souvent, seul un point très précis est retenu d’une longue étude sur la démographie française (ici l’espérance de vie à la naissance) publié par l’INSEE. Au delà du débat sur les chiffres, j’ai pu comprendre ce que l’INSEE mettait derrière le terme “espérance de vie à la naissance” à savoir “la durée de vie moyenne d’une génération fictive qui connaîtrait tout au long de son existence les conditions de mortalité par âge de l’année considérée“. En fait, je n’avais pas compris jusqu’à ce soir ce que cachait cette idée de “génération fictive soumise aux conditions de mortalité de l’année“. En formalisant un peu, ça veut dire que les calculs sont fait de manière transversale, et pas du tout prospective. Car pour faire simple, des données de mortalité, ce sont des données sous la forme suivante,

A différentes dates, disons qu’on fait un recensement. On a alors le nombre de personnes en vie, et décédée à un âge donné. Par exemple, p_{x,t} désigne ici la probabilité qu’une personne d’âge x survive jusqu’à son prochain anniversaire, l’année t. Un peu de calculs permet de voire que

e= \sum_{x=0}^{120} \prod_{n=0}^x p_{n,t}peut être interprété comme une espérance de vie à la naissance. C’est la version discrète de la formule

\mathbb{E}[T]=\int_0^\infty \mathbb{P}(T>t)dt(quand on travaille sur une variable aléatoire positive) pour la somme, et le produit est juste lié au fait que

\mathbb{P}[T>t]=\mathbb{P}[T>t\vert T>t-1]\cdot \mathbb{P}[T>t-1\vert T>t-2] \cdots

(avec des mots, la probabilité d’atteindre 60 ans, c’est la probabilité de survivre à son 1er anniversaire, puis à son second, …, puis à son 59ème). Bref, on a toutes ces probabilités de survie dans nos données. Des données, on peut en trouver par exemple sur mortality.org.

> library(demography)
> france <- hmd.mx("FRACNP", username, password, "France")

La méthode décrite par l’INSEE ressemble à la suivante:

On considère toutes les probabilités de survie brute, et on cumule,

> kPx_obs=1-france$rate$total[,"2013"]
> sum(cumprod(kPx_obs))
[1] 81.24053

(je ne fais pas ici de distinction Homme/Femme, les plus curieux pourront creuser par eux-mêmes). On note que si en 2013, il y a eu une grosse grippe, qui a touché tout le monde, l’espérance de vie en sera forcement impactée ! On n’a pas ici de correction ici, pour lisser les décès tout au long de l’année, mais c’est possible de le faire.

Une autre stratégie peut être d’ajuster un modèle sur les données de mortalité, par exemple le modèle de Lee & Carter (1992), qui dit que la probabilité qu’une personne d’âge x décède avant son prochain anniversaire, l’année t, est de la forme

1-p_{x,t}=\exp[a_x+b_x\cdot \kappa_t]
On peut ajuster ce modèle et considérer alors

\widehat{p}_{x,t}=1-\exp[\widehat{a}_x+\widehat{b}_x\cdot \widehat{\kappa}_t]

On peut utiliser une version transversale comme auparavant,

> france.LC <- lca(france,adjust="e0",years=1950:2013)
> france.fcast <- forecast(france.LC,jumpchoice="actual",h=150)
> Ax=france.LC$ax
> Bx=france.LC$bx
> Kt=c(france.LC$kt,france.fcast$kt$mean)
> years=c(1950:2013,2014:2163)

L’espérance de vie à la naissance (telle qu’on l’avait définie auparavant) est alors

> kPx_model=1-exp(Ax+Bx*Kt[years==2013])
> sum(cumprod(kPx_model))
[1] 84.27459

J’avoue être ici surpris par l’écart avec le précédant. Mais pourquoi pas…

Enfin, autre méthode, basée sur un modèle prospectif. L’idée est de noter que le vieillissement d’un individu, dans une cohorte donnée, se fait en diagonale, et pas verticalement, dans le diagramme de Lexis

On utilise ici le fait que quelqu’un né en 2015 autant 40 ans en 2055, et qu’il faut prévoir la probabilité de décéder l’année de ses 40 ans dans 40 ans. En intégrant les gains sur l’espérance de vie, on obtient cette fois l’espérance de vie à la naissance

> kPx_pred=1-exp(Ax+Bx*Kt[years%in 13:(2012+length(Ax))])
> sum(cumprod(kPx_pred))
[1] 90.15392

pour quelqu’un né en 2013. Cette fois, on a un gain vraiment significatif, correspondant à des améliorations des conditions de vie avec le temps, et la baisse de la mortalité.

Je reste toutefois surpris que l’INSEE continue à utiliser une version transversale (et pas du tout prospective) dans ses calculs d’espérance de vie, même si, comme le notait @l’espérance de vie à la naissance est une statistique internationale publiée par tous les pays, et donc doit être calculée de façon simple (ou pas trop complexe). Et je renvoie vers un vieux billet pour ceux qui veulent aller plus sur ces calculs prospectifs, et le fameux “tous les ans, on gagne un trimestre d’espérance de vie“. Ou encore un article sur les tables de mortalité, paru dans Risques en 2007.

Décomposition(s) de Séries Temporelles

La semaine dernière, pour le TP du cours de série temporelle, on avait parlé de décomposition de séries temporelles. On part d’une série temporelle , et on va considérer une décomposition de la forme

par exemple une tendance linéaire, et un cycle annuel,

Considérons la série suivante

> autoroute=read.table(
+ "http://freakonometrics.free.fr/autoroute.csv",
+ header=TRUE,sep=";")
>  a7=autoroute$a007
>  X=ts(a7,start = c(1989, 9), frequency = 12)

correspondant à du trafic sur l’autoroute (A7)

>  plot(X,xlim=c(1989,2000))

Pour estimer la tendance, on va faire une régression linéaire, sur le temps, à partir de la base suivante

>  T=time(X)
>  base=data.frame(Y=a7,X1=as.numeric(T))

Le modèle est ici

>  reg1=lm(Y~X1,data=base)

ce qui permet de faire une prévision sur les données passées, mais aussi futures,

>  tendance1=predict(reg1)
>  T=as.numeric(T)
>  lines(T,tendance1,col="red",lwd=3)
>  futur=seq(max(T)+1/12,length=36,by=1/12)
>  tendance2=predict(reg1,
+  newdata=data.frame(X1=futur))
>  lines(futur,tendance2,col="red")

Voilà pour la tendance. Pour le cycle, on va ajuste un modèle de la forme suivante

sur la série des résidus

>  base$Y2=residuals(reg1)

Cette fois, on va définir le mois de la date

> base$mois=round(base$X1-t
+ trunc(base$X1),digit=4)

Et on régresse le résidu sur le mois

>  reg2=lm(Y2~0+as.factor(mois),data=base)

Là encore, on peut faire une régression sur les données passées

>  cycle1=predict(reg2)
>  lines(T,tendance1+cycle1,col="blue",lwd=3)

mais aussi sur les données futures,

>  moisfutur=round(futur-trunc(futur),4)
>  cycle2=predict(reg2,newdata=
+  data.frame(mois=moisfutur))
>  lines(futur,tendance2+cycle2,col="blue")

Cela dit, pour les données futures, on peut aussi tenter de rajouter une région de confiance, en tenant compte de l’erreur d’ajustement,

>  E=residuals(reg2)
>  hist(E)
>  1.96*sqrt(var(E))
[1] 5637.897

On peut alors facilement ajouter notre région de confiance avec

> polygon(c(futur,rev(futur)),
+ c(tendance2+cycle2-1.96*sqrt(var(E)),
+ rev(tendance2+cycle2+1.96*sqrt(var(E)))),
+ col="yellow",border=NA)

Et voilà…. On peut quand même vérifier que le bruit est effectivement non-autocorrélé, et c’est le cas ici

>  plot(acf(E))

Comme on l’a vu en TP, ajuster une droite de régression c’est bien ici, mais peut être qu’on peut envisager un ajustement local. Par exemple prendre la moyenne autour de la date qui nous intéresse, +/- 6 mois (techniquement, si on rajoute 6 mois avant et après, on aura un soucis, alors on fait une petite correction au bord),

ce qui correspond à un lissage par moyenne mobile.

>  n=length(X)
>  F=rep(NA,n)
>  for(t in 7:(n-7)) F[t]=1/12*(sum(X[(t-5):(t+5)])+ (X[t-6]+X[t+6])/2)
>  lines(T,F,lwd=3,col="purple")

On voit ici – par notre ajustement non-paramétrique – que l’ajustement d’une tendance linéaire n’était pas si mauvais que ça

Par contre, pour des séries où la tendance n’est clairement pas linéaire, le lissage par moyenne mobile permet d’avoir une bonne estimation de la tendance de fond

> sncf=read.table(
+ "http://freakonometrics.free.fr/sncf.csv",
+ header=TRUE,sep=";")
> train=as.vector(t(as.matrix(sncf[,2:13])))
> X=ts(train,start = c(1963, 1), frequency = 12)
> plot(X)  
> n=length(X)
> F=rep(NA,n)
> for(t in 7:(n-7)) F[t]=1/12*(sum(X[(t-5):(t+5)])+(X[t-6]+X[t+6])/2)
> lines(as.numeric(time(X)),
+ F,lwd=3,col="purple")

On continuera en cours cette semaine autour du lissage, par moyenne mobile, et avec du lissage exponentiel ensuite. A suivre donc.

Confidence Regions for Parameters in the Simplex

Consider here the case where, in some parametric inference problem, parameter  is a point in the Simplex,

For instance, consider some regression, on compositional data,

> library(compositions)
>  data(DiagnosticProb)
>  Y=DiagnosticProb[,"type"]-1
>  X=DiagnosticProb[,c("A","B","C")]
>  model = glm(Y~ilr(X),family=binomial)
>  b = ilrInv(coef(model)[-1],orig=X)
>  as.numeric(b)
[1] 0.3447106 0.2374977 0.4177917

We can visualize that estimator on the simplex, using

>  tripoint=function(s){
+    p=s/sum(s)
+    abc2xy(matrix(p,1,3))
+  }

>  lab=LETTERS[1:3]
>  xl=c(-.1,1.25)
>  yl=c(-.1,1.15)
>  library(trifield)
>  A=abc2xy(matrix(c(1,0,0),1,3)) 
>  B=abc2xy(matrix(c(0,1,0),1,3))
>  C=abc2xy(matrix(c(0,0,1),1,3)) 
>  plot(0:1,0:1,col="white",
+  xlim=xl,ylim=yl,xlab="",ylab="",axes=FALSE)
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

If we want to compute a ‘confidence region’, we can either use Bayesian models (with a Dirichlet distribution as prior distribution), or use bootstrap. We will use here the second idea

>  MB=matrix(NA,1e4,2)
>  for(sim in 1:1e4){
+    idx=sample(1:nrow(DiagnosticProb),
+    size=nrow(DiagnosticProb),replace=TRUE)
+  Y=DiagnosticProb[idx,"type"]-1
+  X=DiagnosticProb[idx,c("A","B","C")]
+  model = glm(Y~ilr(X),family=binomial)
+  MB[sim,]=tripoint(as.numeric(
+    ilrInv(coef(model)[-1],orig=X)))}

To get some ‘confidence region’, we can then use the bagplot, to get either a region where 50% of the boostraped estimators are, or 95%,

>  library(aplpack)
> P1=bagplot(MB[,1],MB[,2], factor =1.96, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 
> P2=bagplot(MB[,1],MB[,2], factor =0.67, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 

Then we can easily plot those two regions,

>  plot(0:1,0:1,col="white")
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  polygon(P1$hull.loop,col="yellow",border=NA)
>  polygon(P2$hull.loop,col="orange",border=NA)
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

 

Multiattribute, Risk and Utility

As discussed already, in June 2016, with Olivier L’Haridon, we will organize a (small) conference, in Rennes, on risk models in a multi-attribute framework. Related to that conference, we have a working group on related topics. At the end of September, I gave a brief survey on multivariate distributions. And yesterday, Olivier gave the first part of survey on multivariate decision making. The second part will be in two weeks,

Collectionner les Figurines Lego avec ses Ami(e)s

Dans mon précédant billet, je regardais combien de figurines Lego je devais acheter pour avoir toute la série, avec ma fille. Pour la série de 16, il fallait, en moyenne, acheter 54 sachets, à 3 euros pièces. Dans les commentaires, @Guillaume et @Gaelle me faisaient noter que les échanges, ça pouvait aider.

Bon, tout d’abord j’aime pas qu’on donne des conseils à mes filles pour avoir raison contre moi ! Parce que oui, c’est sûr qu’en s’y mettant à plusieurs, le coût individuel sera plus faible… Ensuite, parce que les maths me semblent plus compliqués…. Cela dit, en simulant des parties, ça tombe assez vite…

Pour commencer, je fais une petite fonction (histoire d’alléger mes codes) pour incrémenter mes comptages de figurines, à chaque fois que j’en achète une

> f=function(mat,i){
+   mat[i]=mat[i]+1
+   return(mat)
+ }

Ensuite, on peut avoir une petite fonction, correspondant à une simulation, qui me renvoie le nombre de figurine de chaque type que j’ai du acheter pour avoir les 16 une fois. Ou plus. En fait, je veux avoir en tout k fois toutes les figurines (ce qui permet de partager parmi k ami(e)s)

> tirage=function(k=1){
+ test=FALSE
+ cpte=0
+ M=rep(0,16)
+ while(!test){
+   no=sample(1:16,size=1)
+   M=f(M,no)
+   cpte=cpte+1
+   if(min(M)==k) test=TRUE
+ }
+ return(list(compte=M,
+             nbre=cpte))
+ }

Maintenant on peut faire des simulations. Par exemple, pour retrouver le fait qu’il faut, en moyenne, acheter 54 sachet pour que ma fille (seule) ait toutes les figurines

> N=unlist(lapply(1:10000,function(x) tirage(1)$nbre))
> mean(N)
[1] 54.1854

Si maintenant on le fait avec mes deux filles,

> N=unlist(lapply(1:10000,function(x) tirage(2)$nbre))
> mean(N)
[1] 82.7961

ce qui est (beaucoup) moins que le double. Si au lieux de faire ça entre mes deux filles, mais que ma fille collectionne les figurines avec une amie à elle (et que les papas partagent les frais), on a en moyenne 41 sachets à acheter, chacun,

> mean(N)/2
[1] 41.39805

et pour trois, on passe à 36 sachets

> N=unlist(lapply(1:10000,function(x) tirage(3)$nbre))
> mean(N)/3
[1] 36.04363

etc

> V=unlist(lapply(1:10,function(k) mean(unlist(lapply(1:10000,function(x) tirage(k)$nbre)))))
> plot(1:10,V/1:10,type="b",pch=19)

On peut aller plus loin, en regardant les quantiles. Avec 4 ami(e)s, si on achète tous 45 sachets, on a 90% de chances que tout le monde ait les 16 figurines

D’ailleurs, si on regarde l’évolution du quantile (à 90%) sur la moyenne, on voit qu’il baisse vite avec la taille…. Avec 2 ami(e)s, il faut 35% de sachets en plus pour s’assurer d’avoir toutes les figurines dans 90% de cas (par rapport au nombre moyen) alors qu’il en faut juste 20% de plus avec 8 ami(e)s.

Maintenant, si quelqu’un a le courage (et le temps) de faire les maths, je suis preneur !

TP, Séries Temporelles

Deux bases pour illustrer le TP de séries temporelles. La première est du trafic sur autoroute, en France

> autoroute=read.table(
+ "http://freakonometrics.free.fr/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> X=ts(a7,start = c(1989, 9), frequency = 12)
> T=time(X)

et on utilisera aussi du traffic sur le réseau de la SNCF

> sncf=read.table(
+   "http://freakonometrics.free.fr/sncf.csv",
+   header=TRUE,sep=";")
> train=as.vector(t(as.matrix(sncf[,2:13])))
> X=ts(train,start = c(1963, 1), frequency = 12)

 

Dualité Individus / Variables en ACP

Dans le cours d’analyse des données, on va commencer l’analyse en composantes principales. Comme on l’a déjà évoqué à plusieurs reprises, il y a une dualité entre les individus dans l’espace des variables, et les variables dans l’espace des individus. Histoire de faire des dessins, on va se contenter de 3 individus, avec 3 variables,

> 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
>  M=M[1:3,1:3]
>  n=nrow(M)  
>  M_s=scale(M,scale=TRUE)

 

On peut commencer par visualiser les individus

Les directions principales sont ici

> U=eigen(t(M_s)%*%(M_s))$vectors

 

Les deux premières nous donne le plan principal

La projection sur le plan donne ici

Mais on peut aussi visualiser les données dans l’autre dimension.

Rappelons que les variables ont été centrées et réduites. La encore, l’analyse spectrale nous donne les vecteurs principaux (qui vont maximiser l’inertie en projection),

>  V=eigen((M_s)%*%t(M_s))$vectors

 

Là encore, on peut s’intéresser au plan principal

Là encore on peut regarder la projection dans cette dimension,

Comme on a centré et réduit nos variables, on voit apparaître le “cercle de corrélations“.

Le principe de dualité dit que l’on peut superposer les projections

 

ce qui permettra de faire une réelle analyse des données