Tag Archives: computer

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)

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.

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)

 

Manipuler des Séries Chronologiques

Ce matin en cours, on a vu un peu comment manipuler des séries temporelles, et en particulier comme jouer avec l’échelle de temps (passer de données hebdomadaire à de données mensuelles, par exemple).

Pour commencer, travaillons un peu sur les données extraites de google.fr/trends/ par exemple sur le mot clé fièvre,

report = read.csv(
"http://freakonometrics.free.fr/report-fievre.csv",
  skip=6,header=FALSE,nrows=626)
 X=report$V2
 plot(X,type="l")

Au lieu d’avoir un vecteur, on peut créer un format de série temporelle (hebdomadaire )

 X_h=ts(X,start=c(2004,3),frequency=52)
 plot(X_h)

On va essayer d’agréger les données pour les rendre mensuelles

 date1=substr(as.character(report$V1),1,10)
 date2=substr(as.character(report$V1),14,23) 
 d1=as.Date(date1,"%Y-%m-%d")
 d2=as.Date(date2,"%Y-%m-%d")
 d1[2]-d1[1]
 Vect_J=Vect_M=Vect_X=Vect_A= NULL
 for(i in 1:length(d1)){
   jours = seq(d1[i],d2[i],length=7)
   mois = as.POSIXlt(jours)$mon+1
   annees = as.POSIXlt(jours)$year+1900
   donnees = rep(X[i],7)
   Vect_J = c(Vect_J,jours)
   Vect_M = c(Vect_M,mois)
   Vect_A = c(Vect_A,annees)
   Vect_X = c(Vect_X,donnees)
 }
 base=data.frame(J=Vect_J,M=Vect_M,
                 A=Vect_A,X=Vect_X)
 head(base)
 base$AM = base$A*100+base$M
 head(base)
 Y=tapply(base$X,as.factor(base$AM),mean)
 Z=ts(as.numeric(Y[1:144]),
 start=c(2004,1),frequency=12)
 plot(Z)

On peut superposer les deux d’ailleurs (il y a peut être un léger décalage sur la série hebdomadaire)

plot(X_h,col="blue")
lines(as.numeric(time(Z)),as.numeric(Z),
col="red",lwd=2)

Continue reading Manipuler des Séries Chronologiques

Analyse des Données et Cartes

Vendredi, on continuera en cours la classification (non supervisée) et en particulier les méthodes hiérarchiques, et les arbres. On va utiliser la base suivante, avec les résultats des élections présidentielles (premier tour) de 2012,

elections2012 <- read.csv("http://komodo.regardscitoyens.org/public/presidentielles2012/resultats/resultats_departements_final_T1.csv",
sep=";",header=TRUE,dec=",")

La base qu’on garde contient juste les voix exprimées, en pourcentage

voix <- which(substr(names(elections2012),1,12)=="X..Voix.Exp.")
X <- as.matrix(elections2012[,voix])
colnames(X) <- substr(names(elections2012)[voix],13,nchar(names(elections2012)[voix]))
rownames(X) <- elections2012[,3]

Parfois, il convient de normaliser les données. Ici, sur les département, ça n’est pas nécessaire. Par contre, sur les candidats, ça le serait (comme pour l’ACP en fait).  Pour visualiser les distances (entre départements), on peut utiliser

heatmap(X)

On peut faire une CAH, sur notre matrice de distance entre rangs

cah <- hclust(dist(X))
plot(cah,cex=.6)

et si on souhaite garder 5 groupes (par exemple) on utilise

rect.hclust(cah,k=5)
groups.5 <- cutree(cah,5)

On peut aussi visualiser les classes avec la fonction

library(dendroextras)
plot(colour_clusters(cah,k=5))

Pour comprendre qui se trouve dans nos groupes, on peut utiliser

aggregate(X,list(groups.5),mean)

qui va nous renvoyer les “votes  moyens” pour chaque candidat, dans chaque groupe. Maintenant, pour en finir avec la visualisation, on peut utiliser le code suivant (un peu long, certes) qui permettra de visualiser sur un carte les différents groupes

carte_classe <- function(groupes){
library(stringr)
elections2012$dep <- elections2012$Libellé.du.département
elections2012$dep <- tolower(elections2012$dep)
elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
library(maps)
france<-map(database="france")
france$dep <- france$names
france$dep <- tolower(france$dep)
france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
corresp_noms <- elections2012[, c("Libellé.du.département", "dep")]
corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
col2001<-groupes+1
names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms$Libellé.du.département)]
color <- col2001[match(france$dep, names(col2001))]
map(database="france", fill=TRUE, col=color)
}

carte_classe(groups.5)

Avec ces fonctions, on devrait pouvoir tester différentes méthodes pour constituer des groupes.

Mensualiser une série hebdomadaire

Pour le devoir de série temporelles, les données fournies par Google trend sont hebdomadaires, ce qui peut rendre la modélisation compliquée. Comme on l’a évoqué en cours, il peut être plus simple de travailler sur des données mensuelles. La petite fonction suivante permet de transformer les données pour avoir des données mensuelles (avec des moyennes par mois, ce qui fait qu’un mois de 28 jours et un mois de 31 jours sont comparables). Histoire d’illustrer, prenons un mot au hasard, disons…. épilation. On commence par sauver le fichier csv, et on le lit sous R,

EPILATION=read.table("/home/charpentier/report-epilation.csv", skip=4,header=TRUE,sep=",",nrows=426)

La petite fonction suivant va nous aider à convertir la base en une série temporelle mensuelle,

H2M=function(BASE){ 
X=BASE[,2] 
date1=substr(as.character(BASE$Semaine),1,10) 
date2=substr(as.character(BASE$Semaine),14,23) 
D1=as.Date(date1,"%Y-%m-%d") 
D2=as.Date(date2,"%Y-%m-%d") 
vm=vy=N=NA 
for(t in 1:length(D1)){ 
mois=seq(D1[t],D2[t],length=7) 
vm=c(vm,as.POSIXlt(mois)$mon+1) 
vy=c(vy,as.POSIXlt(mois)$year+1900) 
N=c(N,rep(X[t],7))} N=N[-1]; 
vm=vm[-1]; 
vy=vy[-1] 
YM=vy*100+vm 
Z=tapply(N,as.factor(YM),mean) 
Zts=ts(Z,start=c(2004,1),frequency=12) 
return(Zts)}

Si on travaille sur la série hebdomadaire, on a la série suivante

hebdo=ts(EPILATION$épilation,start=2004, frequency=52)

La fonction d’autocorrélation est alors

acf(hebdo,160)

Maintenant, on peut travailler sur les données mensuelles. On utilise alors

mensuel=H2M(EPILATION)

Le graphique est alors le suivant

La fonction d’autocorrélation est cette fois la suivante

acf(mensuel,40)

On retrouve le comportement cyclique, avec une saisonnalité annuelle.

Inequalities, course 3

Tomorrow, we will discuss inequality indices, from a statistical perspective, and also an axiomatic point of view. In order to illustrate, we will use to following dataset,

> income <- read.csv("http://www.vcharite.univ-mrs.fr/pp/lubrano/cours/fes96.csv",sep=";",header=FALSE)$V1

Slides can be found online. Since it is the first year I give this course, all comments are welcome…

Analyse des Correspondances, suite

Hier soir, on voyait comment faire une analyse des correspondances à partir du tableau de contingence. Mais on peut avoir le problème autrement, à partir des individus. Certes, ces derniers ne sont pas observés (vraiment), mais peu importe.

> data(HairEyeColor)
> N = HairEyeColor[,,"Male"] + HairEyeColor[,,"Female"]
> Hair = rep(rep(rownames(N),ncol(N)),
+ as.vector(N))
> Eye =  rep(colnames(N),apply(N,2,sum))
> df = data.frame(Ind=1:sum(N),Hair,Eye)
> tail(df)
    Ind  Hair   Eye
587 587 Blond Green
588 588 Blond Green
589 589 Blond Green
590 590 Blond Green
591 591 Blond Green
592 592 Blond Green

Continue reading Analyse des Correspondances, suite

Subjective Ways of Cutting a Continuous Variables

You have probably seen @coulmont‘s maps. If you haven’t, you should probably go and spend some time on his blog (but please, come back afterwards, I still have a story to tell you). Consider for instance the maps we obtained for a post published in Monkey Cage, a few months ago,

The codes were discussed on a blog post (I spent some time on the econometric model, not really on the map, by that time).

My mentor in cartography, Reka (aka @visionscarto) taught me that maps were always subjective. And indeed.

Continue reading Subjective Ways of Cutting a Continuous Variables

Mais que s’est-il passé pendant la Première Guerre Mondiale?

La réponse courte est que des gens sont morts. Beaucoup. Cela étant dit, on ne dit pas grand chose. On peut comparer les pyramides des âges pour mieux comprendre ce qui a pu se passer. Juste avant la guerre (en 1913), la pyramide des âges ressemblait à ça (en utilisant les données de mortality.org)

> EXPO  <- read.table(
+ "http://freakonometrics.free.fr/Exposures-France.txt", header=TRUE,skip=2)
> EM=EXPO$Male
> EF=EXPO$Female
> Y= EXPO$Year
> A= EXPO$Age
> I=which(A=="110+")
> base=data.frame(Female=EF,Male=EM,Y=Y,Ages=A)
> base=base[-I,]
> France1913=base[base$Y==1913,]
> France1919=base[base$Y==1919,]
> France1913$Ages=as.numeric(
+ as.character(France1913$Ages))
> France1919$Ages=as.numeric(
+ as.character(France1919$Ages))
> France1913=France1913[,c("Male","Female",
+ "Ages")]
> library(pyramid)
> plot(c(0,100), c(0,100), type="n", 
+ frame=FALSE, axes=FALSE, xlab="", ylab="",
+ main="Pyramide des Ages, France, 1913")
> pyramidf(France1913, frame=c(10, 75, 0, 90), 
+ Clab="", Lcol="skyblue", Rcol="pink",
+ Cstep=10, Laxis=0:4*60000, AxisFM="d")

En revanche, juste après la guerre (en 1919), la pyramide des âges des âges ressemblait à celle là Continue reading Mais que s’est-il passé pendant la Première Guerre Mondiale?

Formation R à la CIMA, au Gabon

Je vais bientôt descendre une semaine au Gabon, au siège de la CIMA, la Conférence Interafricaine des Marchés d’Assurance, pour une formation à R.

Souhaitant un cours aussi interactif que possible, je n’ai pas prévu de transparents. Je mettrais toutefois en ligne des documents, si nécessaire. J’avais fait un billet, la semaine passée, pour expliquer comment importer les tables de mortalité. On travaillera sur ces bases pour apprendre à manipuler les fonctions de base de R, lundi.

Sinon, en guise d’introduction, je mentionnerais quelques livres ou notes de cours, en pdf,

  • “R pour les débutants” par Emmanuel Paradis, (pdf)
  • “Introduction à la programmation en R, Quatrième édition” par Vincent Goulet (pdf)
  • “Brise Glace-R (ouvrir la voie aux pôles statistiques)” par Andrew Robinson et Arnaud Schloesing (pdf)
  • “Introduction à R” par Julien Barnier (pdf)
  • “Aide mémoire R” par Mayeul Kauffmann (pdf)
  • “Lire ; Compter ; Tester… avec R”  (pdf)
  • “L’Actuariat avec R” par Arthur Charpentier et Christophe Dutang (pdf)

et je peux aussi mentionner les slides d’Ewen Gallic,

(pour les derniers transparents, je doute toutefois qu’on parlera de ggplot2, je pense qu’on se limitera aux fonctions de base). A suivre…

Shapefiles from Isodensity Curves

Recently, with @3wen, we wanted to play with isodensity curves. The problem is that it is difficult to get – numerically – the equation of the contour (even if we can easily plot it). Consider the following surface (just for fun, in order to illustrate the idea)

> f=function(x,y) x*y+(1-x)*(1-y)
> u=v=seq(0,1,length=21)
> v=seq(0,1,length=11)
> f=outer(u,v,f)
> persp(u,v,f,theta=angle,phi=10,box=TRUE,
+ shade=TRUE,ticktype="detailed",xlab="",
+ ylab="",zlab="",col="yellow")

For instance, assume that we want to locate areas where the density exceed 0.7 (here in the lower left corner, SW, and the upper right corner, NE)

> image(u,v,f)
> contour(u,v,f,add=TRUE,levels=.7)

Continue reading Shapefiles from Isodensity Curves

Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction is now online, on hal.archives-ouvertes.fr/.

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and difficult to implement. We provide a simple technique — based of properties of Gaussian kernels — to efficiently compute weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. We illustrate the use of that technique to visualize hot spots of car accidents and campsite locations, as well as location of bike thefts.

There are new applications, and new graphs, too

Most of the codes can be found on https://github.com/ripleyCorr/Kernel_density_ripley (as well as datasets).