Tag Archives: computer

Interpréter un test de Wald

Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l’interprétation d’un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données

> p=.5

En l’occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face)

> n=20
> set.seed(1)
> echantillon=sample(c("pile","face"),
+    size=n,replace=TRUE,
+    prob=c(p,1-p))
> echantillon
[1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile"

Classiquement, on va recoder en https://latex.codecogs.com/gif.latex?\{0,1\}

> X=(echantillon=="pile")*1
> X
[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1

La vraisemblance est ici donnée par

> L=function(p) prod(dbinom(X,1,p))
> logL=function(p) sum(log(dbinom(X,1,p)))
> xp=seq(0,1,by=.01)
> yL=Vectorize(L)(xp)
> plot(xp,yL,type="l")

mais on va plutôt travailler sur la log-vraisemblance

> ylogL=Vectorize(logL)(xp)
> plot(xp,ylogL,type="l")

Continue reading Interpréter un test de Wald

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

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

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

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

Où ont été observées les fortes hausses du prix de l’essence ?

Vendredi, dans mon billet “le prix de l’essence en cas de pénurie“, j’essayais de visualiser l’explosion du prix de l’essence dans certaines stations services, ces dernières semaines. Mais plus généralement, on peut se demander où ont été observées les fortes hausses.

A partir de la base téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de récupérer les informations les prix, par jour, dans toutes les stations essence de France. Pour ça, on utilise la fonction suivante,

> spatial=function(dt){
+  base=NULL
+  for(no in 1:length(l)){  
+    prix=list()
+    date=list()
+    j=0
+    for(i in 1:length(l[[no]])){
+    v=names(l[[no]])
+    if(!is.null(v[i])){
+    if(v[i]=="prix"){
+    j=j+1
+    date[[j]]=as.character(l[[no]][[i]]["maj"])
+    }}
+  }
+ n=j
+ D=as.Date(substr(unlist(date),1,10),"%Y-%m-%d")
+ k=which(D==D[which.max(D[D<=dt])])
+ if(length(k)>0){
+ B=Vectorize(function(i) l[[no]][[k[i]]])
+   (1:length(k))
+ if("nom" %in%  rownames(B)){  
+   k=which(B["nom",]=="Gazole")
+   prix=as.numeric(B["valeur",k])/1000
+   if(length(prix)==0) prix=NA
+   base1=data.frame(indice=no,
+ lat=as.numeric(l[[no]]$.attrs["latitude"])
+    /100000,
+ lon=as.numeric(l[[no]]$.attrs["longitude"])
+    /100000,
+ gaz=prix,
+ cp=l[[no]]$.attrs["cp"])
+ base=rbind(base,base1)
+ }}}
+ return(base)}

Par exemple, on peut récupérer les prix fin mai

> B1=spatial(as.Date("2016-05-31"))

et début mai,

> B2=spatial(as.Date("2016-05-01"))

En fusionnant ces deux bases, on peut récupérer la variation du prix, sur un mois,

> names(B1)=c("indice","lat","long","fin","cp")
> names(B2)=c("indice","lat","long","debut","cp")
> B=merge(B1,B2)
> B$var=(B$fin-B$debut)/B$debut*100

Plus particulièrement, on va s’intéresser aux stations dans lesquelles le prix a augmenté de plus de 10%,

> idx=which((B$lon>(-10))&(B$lon<20)&
+             (B$lat>35)&(B$lat<55))
> B=B[idx,]
> B=B[!is.na(B$var),]
> B$Y=(B$var>10)*1

Notons déjà que dans plus de 90% des stations, le prix a augmenté entre début mai, et fin mai

> mean(B$var > 0)
[1] 0.9288522

mais les hausses de plus de 10% restent rares (moins de 5% des stations essence)

> mean(B$var > 10)
[1] 0.03305452

La distribution de la variation des prix (en %) est la suivante

> plot(density(B$var),xlim=c(-5,20))

Si on regarde maintenant où ont été observées ces fortes hausses, on observe qu’elles peuvent avoir lieu n’importe où

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.25),
+ rgb(1,0,0,.95))[1+(B$var>10*1)],cex=.5)

et c’est encore plus flagrant quand on regarde les hausses de plus de 5%

> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.35),
+ rgb(1,0,0,.95))[1+(B$var>5*1)],cex=.5)

Et si on regarde les endroits où la proportion des stations ayant fortement augmenté les prix est la plus forte, on obtient la carte suivante,

Autrement dit la hausse des prix a pu être observée un peu partout en France. On peut alors comparer avec la carte des pénuries, mise en ligne sur https://mon-essence.fr/

Si certaines régions avec des fortes hausses de prix correspondent à des pénuries d’essence, ce n’est pas le cas partout. Bref, comprendre ce qui s’est passé sur les prix ce mois de mai est loin d’être simple… à suivre donc….

Le prix de l’essence en cas de pénurie

Cet après midi, sur le site de 60 millions de consommateurs, Lionel Maugain revenait sur les stations services qui avaient de très fortes variations de prix pendant les jours de pénurie.

En utilisant la base des prix des carburants à la pompe, téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de visualiser ces fortes variations de prix. Pour lire la base, rien de plus simple

> rm(list=ls())
> annee=year=2016
> fichier="PrixCarburants_annuel_2016.xml"
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

Ensuite, on transforme la grosse liste en une série temporelle, à l’aide de la fonction suivante

> time_series=function(no,type_gas="Gazole"){
+    prix=list()
+    date=list()
+    nom=list()
+    j=0
+    for(i in 1:length(l[[no]])){
+      v=names(l[[no]])
+      if(!is.null(v[i])){
+      if(v[i]=="prix"){
+      j=j+1
+      date[[j]]=as.character(l[[no]][[i]]["maj"])
+      prix[[j]]=as.character(l[[no]][[i]]["valeur"])
+      nom[[j]]=as.character(l[[no]][[i]]["nom"])
+     }}
+    }
+    id=which(unlist(nom)==type_gas)
+    n=length(id)
+    jour=function(j) as.Date(substr(date[[id[j]]],1,10),"%Y-%m-%d")
+    jour_heure=function(j) as.POSIXct(substr(date[[id[j]]],1,19), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
+    ext_y=function(j) substr(date[[id[j]]],1,4)
+    ext_m=function(j) substr(date[[id[j]]],6,7)
+    ext_d=function(j) substr(date[[id[j]]],9,10)
+    ext_h=function(j) substr(date[[id[j]]],12,13)
+    ext_mn=function(j) substr(date[[id[j]]],15,16)
+    prix_essence=function(i) as.numeric(prix[[id[i]]])/1000
+    base1=data.frame(indice=no,
+   id=l[[no]]$.attrs["id"],
+   adresse=l[[no]]$adresse,
+   ville=l[[no]]$ville,
+   lat=as.numeric(l[[no]]$.attrs["latitude"])
    /100000,
+   lon=as.numeric(l[[no]]$.attrs["longitude"])
    /100000,
+   cp=l[[no]]$.attrs["cp"],
+   saufjour=l[[no]]$ouverture["saufjour"], 
+   Y=unlist(lapply(1:n,ext_y)),
+   M=unlist(lapply(1:n,ext_m)),
+   D=unlist(lapply(1:n,ext_d)),
+   H=unlist(lapply(1:n,ext_h)),
+   MN=unlist(lapply(1:n,ext_mn)),
+   prix=unlist(lapply(1:n,prix_essence)))
+   base1=base1[!is.na(base1$prix),]
+   date_d=paste(year,"-01-01 12:00:00",sep="")
+   date_f=paste(year,"-06-01 12:00:00",sep="")
+   vecteur_date=seq(as.POSIXct(date_d, format =
+   "%Y-%m-%d %H:%M:%S"),
+   as.POSIXct(date_f, format = 
+   "%Y-%m-%d %H:%M:%S"),by="days")
+   date=paste(base1$Y,"-",base1$M,"-",base1$D,
+   " ",base1$H,":",base1$MN,":00",sep="")
+   date_base=as.POSIXct(date, format = 
+   "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   idx=function(t) sum(vecteur_date[t]>=date_base)
+      vect_idx=Vectorize(idx)(1:length(vecteur_date))
+      P=c(NA,base1$prix)
+      vp=P[1+vect_idx]
+      vp=vp*100/(vp[!is.na(vp)][1])
+      base2=ts(vp,start=year,frequency=365)
+      list(base=base1,ts=base2)
+    }

Dans l’article, la plus forte variation a été observée à Sevran. Pour récupérer la liste des stations services, on utilise

> ville=list()
> for(no in 1:length(l)) 
      ville[[no]]=l[[no]]$ville
> ville=unlist(ville)
> idx=which(ville=="Sevran")

Pour la première station, l’évolution du prix est relativement modérée

> plot(time_series(idx[1])$ts,col="red",type="l")
> abline(h=100,lty=2)

Même si le prix semble avoir augmenté de 5% pendant les dernières semaines. Mais ce n’est effectivement rien en comparaison de l’autre station,

> plot(time_series(idx[2])$ts,col="red",type="l")
> abline(h=100,lty=2)

avec un prix qui a augmenté de 50% en 2 mois, et presque 30% en quelques jours.

A Villeneuve-d’Asq, on retrouve un comportement similaire, mais moins important

> idx=which(ville=="Villeneuve-d'Ascq")
> plot(time_series(idx[3])$ts,col="red",type="l")
> abline(h=100,lty=2)

Et si je tweetais en anglais pour une mauvaise raison

Sur Twitter, assez naturellement, j’ai tendance à tweeter en anglais. J’ai toujours eu l’impression qu’en anglais, la limite des 140 caractères était moins contraignante. Qu’en anglais, la même chose pouvait être dite de manière plus concise.

Histoire de vérifier cette vague intuition, je suis allé chercher des livres sur Gutenberg. Le plus simple est de récupérer un livre en français, au format texte (txt), de compter le nombre de symboles, et de trouver le même livre en anglais. La fonction, à partir d’un code de livre, ressemble à ça

> cpte=function(no){
+     url0=paste(
+         "http://www.gutenberg.org/cache/epub/",
+         no,"/pg",no,".txt",sep="")
+     download.file(url0,destfile="u.txt")
+     s=scan("u.txt",what="character")
+     sum(nchar(s))
+ }

Par exemple, pour David Copperfield, de Dickens, le  code pour la version originale, en anglais, est

> gb=766

et en français, il faut récupérer 2 fichiers

> fr_a=17868
> fr_b=17869

Aussi,

> g=cpte(gb)
essai de l'URL 'http://www.gutenberg.org/cache/epub/766/pg766.txt'
Content type 'text/plain; charset=utf-8' length 1992535 bytes (1.9 MB)
=================================================
downloaded 1.9 MB

Read 288181 items
> f=cpte(fr_a)+cpte(fr_b)
essai de l'URL 'http://www.gutenberg.org/cache/epub/17868/pg17868.txt'
Content type 'text/plain; charset=utf-8' length 1161492 bytes (1.1 MB)
=================================================
downloaded 1.1 MB

Read 191311 items
essai de l'URL 'http://www.gutenberg.org/cache/epub/17869/pg17869.txt'
Content type 'text/plain; charset=utf-8' length 1135661 bytes (1.1 MB)
=================================================
downloaded 1.1 MB

Read 186515 items

le ratio du nombre de caractère en anglais, et en français, est de l’ordre de 90%. Autrement dit, la version anglaise contient 10% de caractères en moins que la version française,

> g/f
[1] 0.9090498

Si on regarde une pièce de théâtre, comme Othello, la version anglaise est 25% plus courte que la version française.

> fr=18179
> gb=1793
> g=cpte(gb)
essai de l'URL 'http://www.gutenberg.org/cache/epub/1793/pg1793.txt'
Content type 'text/plain; charset=utf-8' length 188634 bytes (184 KB)
=================================================
downloaded 184 KB

Read 24951 items
> f=cpte(fr)
essai de l'URL 'http://www.gutenberg.org/cache/epub/18179/pg18179.txt'
Content type 'text/plain; charset=utf-8' length 239158 bytes (233 KB)
=================================================
downloaded 233 KB

Read 36157 items
> g/f
[1] 0.746805

Plus généralement pour les pièces de théâtre de Shakespeare, on valide mon intuition : les versions en anglais sont toujours plus courtes que les versions en français,

Pour faire ce graphique, on peut utiliser la trentaine de livres que j’ai regardé, à l’aide du fichier

> download.file(
+ "http://freakonometrics.free.fr/livre.RData",
+ "livre.RData")
essai de l'URL 'http://freakonometrics.free.fr/livre.RData'
Content type 'text/plain' length 1095 bytes
=================================================
downloaded 1095 bytes

> load("livre.RData")

La base ressemble à ça

> base
           nom        FR      UK language   ratio
1  3 Mousquet.    1110255 1156011       FR 104.12
2  Dorian Grey     398653  404799       UK 101.54
3  Pierre-Jean     248299  228838       FR  92.16
4      Odyssée     691648  651044       GR 118.17
5       Iliade     550924  514953       GR  74.45
6  Jules César     157726  105685       UK  67.00
7      Macbeth     160571   96276       UK  59.95
8      Othello     189984  141881       UK  74.68
9       Hamlet     275629  145380       UK  52.74
10   Com. Mép.     115583   97704       UK  84.53
11  Ant.-Cléo.     185806  130857       UK  70.42
12   Rom.-Jul.     195273  126323       UK  64.69
13    Roi Lear     196185  124483       UK  63.45
14  Sganarelle      50408   63387       FR 125.74
15 OliverTwist     778914  770784       UK  98.95
16  D.Copperf.    1808770 1644262       UK  90.90
17    Assomoir     585237  565444       FR  96.61
18   T. Raquin     351525  353837       FR 100.65
19    20000 L.     745559  518932       FR  69.60
20    80 Jours     367388  329073       FR  89.57
21  5 Semaines     491510  483289       FR  98.32
22 CentreTerre     362617  370584       FR 102.19
23          93     585073  637499       FR 108.96
24  Notre Dame     885048  928256       FR 104.88
25   Jane Eyre     922233  919521       UK  99.70
26    Possédés    1241293 1293319       RU 104.19
27 Guerre-Paix    2409360 2739725       RU 113.71
28     Corbeau      27884   21886       UK  78.48
29  Mme Bovary     585237  565444       FR  96.61

Si on visualise tous les titres, et pas seulement les pièces de Shakespeare, on obtient

Autrement dit, la différence n’est pas aussi significative qu’elle y paraissait. Plus troublant, si on regarde les livres écrit en anglais (en rouge ci-dessous) on observe effectivement que la version française est plus longue. Par contre, pour les livres écrits en français (en bleu), la version anglaise est sensiblement de la même taille…

A suivre donc…

La parade des nains, ou la visualisation des inégalités

Dans mon précédant billet, sur la répartition de la population sur le territoire français, j’utilisais l’indice de Gini, comme mesure d’inégalité. Et je notais que l’indice était passé de 0.45 au début du XIXème siècle à 0.75 aujourd’hui. Mais force est de constater que cet indice n’est pas très visuel (à moins que l’on ne soit à l’aise avec les courbes de Lorenz). En 1970, Jan Pen proposa introduisit la “parade des nains” afin d’illustrer les inégalités. Pour décrire des inégalités de revenus, on représente des individus dont la taille est proportionnelle au revenu.

Il y a quelques (rares) géants, et beaucoup de nains. On peut tracer la courbe associée à cette parade des nains. Formelle, pour , on trace

où  est la fonction de répartition de la variable que l’on étudie. Rappelons que pour une loi lognormale , l’indice de Gini est

De manière duale, à partir d’un indice de Gini, on peut obtenir les paramètres d’une loi lognormale sous-jacente. Par exemple, pour une loi lognormale, associée à un indice de Gini de l’ordre de 0.25, on obtient

> library(ineq)
> gini=0.25
> sigma=qnorm((1+gini)/2)*sqrt(2)
> x=rlnorm(n,0,sigma)
> Pen(x)

Pour des indices de Gini de 0.50 et de 0.75, on obtient les courbes suivantes

> gini=0.5
> sigma=qnorm((1+gini)/2)*sqrt(2)
> n=1e4
> y1=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> gini=0.75
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y2=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)
> plot((1:n)/(n+1),y1,ylim=c(-.5,10))
> rect(0,0,1,1,col="light yellow",border=NA)
> polygon(c(0,(1:n)/(n+1),1),c(0,y1,0),
+ col=rgb(1,0,0,.5))
> polygon(c(0,(1:n)/(n+1),1),c(0,y2,0),
+ col=rgb(0,0,1,.5))
> abline(h=1,lty=2,col="grey")

La courbe jaune est une distribution parfaitement égalitaire. La courbe rouge est obtenue avec un indice de Gini de l’ordre de 0.50 alors que celle en bleu est obtenue avec un indice de 0.75.

Pour reprendre notre exemple de population, répartie sur le territoire, on peut aussi imaginer un territoire linéaire de 1 km de long. On suppose que la zone de gauche est la moins dense, et celle de droite est la plus dense.

  • avec une répartition uniforme, supposons que tout le monde puisse loger dans des immeubles de 5 étages.
  • avec un indice de Gini de 0.50,  les premiers 200m sont des immeubles de 1 étage, les 200m suivants, ce sont des immeubles de 2 étages, etc. A 700m, on obtient les derniers immeubles de 5 étages. Sur la fin, on a des immeubles de 15 étages sur les 50 derniers mètres, et avant un immeuble de 12 étages, et avant encore des immeubles de 10 étages.
  • dans le cas d’un indice de Gini de l’ordre de 0.75, il n’y a rien sur 25% du territoire, puis sur les 25% suivant, des immeubles de 1 étage, etc. 80% des logements ne dépassent pas 5 étages.  Par contre, sur la fin, on des tours relativement haute, dépassant 20 étages.
> n=19
> gini=0.5
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y1=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> gini=0.75
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y2=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> plot((1:n)/(n+1),y1*5,ylim=c(-.5,20))
> for(i in 1:n){
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,1*5,
+ col="light yellow",border=NA)
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,round(y1[i]*5),
+ col=rgb(1,0,0,.5),border=NA)
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,round(y2[i]*5),
+ col=rgb(0,0,1,.5),border=NA)
+ }
> abline(h=1*5,lty=2,col="grey")
> axis(2)

ou, si on zoome, on obtient

La courbe rouge, c’est la répartition de la population, en France au XIXème siècle, ou en Allemagne aujourd’hui (tel que décrit dans mon précédant billet). La courbe bleue représente la France aujourd’hui.

Non-Uniform Population Density in some European Countries

A few months ago, I did mention that France was a country with strong inequalities, especially when you look at higher education, and research teams. Paris has almost 50% of the CNRS researchers, while only 3% of the population lives there.

It looks like Paris is the only city, in France. And I wanted to check that, indeed, France is a country with strong inequalities, when we look at population density.

Using data from sedac.ciesin.columbia.edu, it is possible to get population density on a small granularity level,

> rm(list=ls())
> base=read.table(
+      "/home/charpentier/glp00ag.asc",
+      skip=6)
> X=t(as.matrix(base,ncol=8640))
> X=X[,ncol(X):1]

The scales for latitudes and longitudes can be obtained from the text file,

> #ncols         8640
> #nrows         3432
> #xllcorner     -180
> #yllcorner     -58
> #cellsize      0.0416666666667

Hence, we have

> library(maps)
> world=map(database="world")
> vx=seq(-180,180,length=nrow(X)+1)
> vx=(vx[2:length(vx)]+vx[1:(length(vx)-1)])/2
> vy=seq(-58,85,length=ncol(X)+1)
> vy=(vy[2:length(vy)]+vy[1:(length(vy)-1)])/2

If we plot our density, as in a previous post, on Where People Live,

> I=seq(1,nrow(X),by=10)
> J=seq(1,ncol(X),by=10)
> image(vx[I],vy[J],log(1+X[I,J]),
+ col=rev(heat.colors(101)))
> lines(world[[1]],world[[2]])

we can see that we have a match, between the big population matrix, and polygons of countries.

Consider France, for instance. We can download the contour polygon with higher precision,

> library(rgdal)
> fra=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds",
+ "fr.rds")
> Fra=readRDS("fr.rds")
> n=length(Fra@polygons[[1]]@Polygons)
> L=rep(NA,n)
> for(i in 1:n) L[i]=nrow(Fra@polygons[[1]]@Polygons[[i]]@coords)
> idx=which.max(L)
> polygon_Fr=
+       Fra@polygons[[1]]@Polygons[[idx]]@coords
> min_poly=apply(polygon_Fr,2,min)
> max_poly=apply(polygon_Fr,2,max)
> idx_i=which((vx>min_poly[1])&(vx<max_poly[1]))
> idx_j=which((vy>min_poly[2])&(vy<max_poly[2]))
> sub_X=X[idx_i,idx_j]
> image(vx[idx_i],vy[idx_j],
+       log(sub_X+1),col=rev(heat.colors(101)),
+       xlab="",ylab="")
> lines(polygon_Fr)

We are now able to extract information about population for France, only (actually, it is only mainland France, islands are not considered here… to avoid complicated computations

> library(sp)
> xy=expand.grid(x = vx[idx_i], y = vy[idx_j])
> dim(xy)
[1] 65730     2

Here, we have 65,730 small squares, in France.

> pip=point.in.polygon(xy[,1],xy[,2],
+     polygon_Fr[,1],polygon_Fr[,2])>0
> dim(pip)=dim(sub_X)
> Fr=sub_X[pip]
> sum(Fr)
[1] 58105272

Observe that the total population within the French polygon is close to 60 million people, which is consistent with actual figures. Now, if we look more carefully at repartition over the French territory

> library(ineq)
> Gini(Fr)
[1] 0.7296936

Gini coefficient is rather high (over 70%), but it is also possible to visualize Lorenz curve,

> plot(Lc(Fr))

Observe that in 5% of the territory, we can find almost 54% of the population

> 1-min(LcF$L[LcF$p>.95])
[1] 0.5462632

In order to compare with other countries, consider the

> LC=function(rds="fr.rds"){
+ Fra=readRDS(rds)
+ n=length(Fra@polygons[[1]]@Polygons)
+ L=rep(NA,n)
+ for(i in 1:n) 
L[i]=nrow(Fra@polygons[[1]]@Polygons[[i]]@coords)
+ idx=which.max(L)
+ polygon_Fr=
+      Fra@polygons[[1]]@Polygons[[idx]]@coords
+ min_poly=apply(polygon_Fr,2,min)
+ max_poly=apply(polygon_Fr,2,max)
+ idx_i=which((vx>min_poly[1])&(vx<max_poly[1]))
+ idx_j=which((vy>min_poly[2])&(vy<max_poly[2]))
+ sub_X=X[idx_i,idx_j]
+ xy=expand.grid(x = vx[idx_i], y = vy[idx_j])
+ dim(xy)
+ pip=point.in.polygon(xy[,1],xy[,2],
+     polygon_Fr[,1],polygon_Fr[,2])>0
+ dim(pip)=dim(sub_X)
+ Fr=sub_X[pip]
+ return(list(gini=Gini(Fr),LC=Lc(Fr))
+ }
> FRA=LC()

For instance, consider Germany, or Italy

> deu=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/DEU_adm0.rds","deu.rds")
> DEU=LC("deu.rds")
> ita=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/ITA_adm0.rds","ita.rds")
> ITA=LC("ita.rds")

It is possible to plot Lorenz curve, together,

> plot(FRA$LC,col="blue")
> lines(DEU$LC,col="black")
> lines(ITA$LC,col="red")

Observe that France is clearly below the other ones. Compared with Germany, there is a significant difference

> FRA$gini
[1] 0.7296936
> DEU$gini
[1] 0.5088853

More precisely, if 54% of French people live in 5% of the territory, only 40% of Italians, and 32% of the Germans,

> 1-min(FRA$LC$L[FRA$LC$p>.95])
[1] 0.5462632
> 1-min(ITA$LC$L[ITA$LC$p>.95])
[1] 0.3933227
> 1-min(DEU$LC$L[DEU$LC$p>.95])
[1] 0.3261124

How long could it take to run a regression

This afternoon, while I was discussing with Montserrat (aka @mguillen_estany) we were wondering how long it might take to run a regression model. More specifically, how long it might take if we use a Bayesian approach. My guess was that the time should probably be linear in , the number of observations. But I thought I would be good to check.

Let us generate a big dataset, with one million rows,

> n=1e6
> X=runif(n)
> Y=2+5*X+rnorm(n)
> B=data.frame(X,Y)

Consider as a benchmark the standard linear regression,

> lm_freq = function(n){
+   idx = sample(1:1e6,size=n)
+   reg = lm(Y~X,data=B[idx,])
+   summary(reg)
+ }

Here the regression is a subset of smaller size. We can do the same with a Bayesian approach, using stan,

> stan_lm ="
+ data {
+ int N;
+ vector[N] x;
+ vector[N] y;
+ }
+ parameters {
+ real alpha;
+ real beta;
+ real tau;
+ }
+ transformed parameters {
+ real sigma;
+ sigma <- 1 / sqrt(tau);
+ }
+ model{
+ y ~ normal(alpha + beta * x, sigma);
+ alpha ~ normal(0, 10);
+ beta ~ normal(0, 10);
+ tau ~ gamma(0.001, 0.001);
+ }
+ "

Define then the model

> library(rstan)
> system.time( 
  stanmodel <<- stan_model(model_code = stan_lm))
utilisateur     système      écoulé 
      0.043       0.000       0.043

We want to see how long it might take to run a regression,

> lm_bayes = function(n){
+   idx = sample(1:1e6,size=n)
+   fit = sampling(stanmodel,
+       data = list(N=n,
+                   x=X[idx],
+                   y=Y[idx]),
+       iter = 1000, warmup=200)
+   summary(fit)
+ }

We use the following package to see how long it takes

> library(microbenchmark)
> time_lm = function(n){
+  M = microbenchmark(lm_freq(n),
+      lm_bayes(n),times=50)
+  return(apply( matrix(M$time,nrow=2),1,mean))
+ }

We can now compare the time it took with ten, one hundred, on thousand, and ten thousand observations,

> vN = c(10,100,1000,10000)
> T = Vectorize(time_lm)(vN)

we can then plot it

> plot(vN,T[2,]/1e6,log="xy",col="red",type="b",
+      xlab="Number of Observations",ylab="Time")
> lines(vN,T[1,]/1e6,col="blue",type="b")

It looks like (if we forget about the very small sample) that the time it takes to run a regression is linear, with the two techniques (the frequentist and the Bayesian ones).

And actually, the same story olds for logistic regressions. Consider the following dataset

> n=1e6
> X=runif(n)
> S=-3+2*X+rnorm(n)
> Y=rbinom(n,size=1,prob=exp(S)/(1+exp(S)))
> B=data.frame(X,Y)

The frequentist version of the logistic regression is

> glm_freq = function(n){
+   idx = sample(1:1e6,size=n)
+   reg = glm(Y~X,data=B[idx,],family=binomial)
+   summary(reg)
+ }

and the Bayesian one, using stan,

> stan_glm = "
+ data {
+ int N;
+ vector[N] x;
+ int<lower=0,upper=1> y[N];
+ }
+ parameters {
+ real alpha;
+ real beta;
+ }
+ model {
+ alpha ~ normal(0, 10);
+ beta ~ normal(0, 10);
+ y ~ bernoulli_logit(alpha + beta * x);
+ }
+ "
> stanmodel = stan_model(model_code = stan_glm) )
> glm_bayes = function(n){
+   idx = sample(1:1e6,size=n)
+   fit = sampling(stanmodel,
+        data = list(N=n,
+        x = X[idx],
+        y = Y[idx]),
+        iter = 1000, warmup=200)
+   summary(fit)
+ }

Again, we can see how long it takes to run those regression models

> time_gl m= function(n){
+   M = microbenchmark(glm_freq(n),
+   glm_bayes(n),times=50)
+   return(apply( matrix(M$time,nrow=2),1,mean))
+ }

 

Where People Live, part 2

Following my previous post, I wanted to use another dataset to visualize where people live, on Earth. The dataset is coming from sedac.ciesin.columbia.edu. We you register, you can download the database

> base=read.table("glp00ag15.asc",skip=6)

The database is a ‘big’ 1440×572 matrix, in each cell (latitude and longitude) we have the population

>  X=t(as.matrix(base,ncol=1440))
>  dim(X)
[1] 1440  572

The dataset looks like

> image(seq(-180,180,length=nrow(X)),
+ seq(-90,90,length=ncol(X)),
+ log(X+1)[,ncol(X):1],col=rev(heat.colors(101)),
+ axes=FALSE,xlab="",ylab="")

Now, if we keep only place where people actually live (i.e. removing cold desert and oceans) we get

> M=X>0
> image(seq(-180,180,length=nrow(X)),
+ seq(-90,90,length=ncol(X)),
+ M[,ncol(X):1],col=c("white","light green"),
+ axes=FALSE,xlab="",ylab="")

Then, we can visualize where 50% of the population lives,

> Order=matrix(rank(X,ties.method="average"),
+ nrow(X),ncol(X))
> idx=cumsum(sort(as.numeric(X),
+ decreasing=TRUE))/sum(X)
> M=(X>0)+(Order>length(X)-min(which(idx>.5)))
> image(seq(-180,180,length=nrow(X)), + seq(-90,90,length=ncol(X)), + M[,ncol(X):1],col=c("white",
+ "light green",col="red"), + axes=FALSE,xlab="",ylab="")

50% of the population lives in the red area, and 50% in the green area. More precisely, 50% of the population lives on 0.75% of the Earth,

> table(M)/length(X)*100
M
         0          1          2 
69.6233974 29.6267968  0.7498057

And 90% of the population lives in the red area below (5% of the surface of the Earth)

> M=(X>0)+(Order>length(X)-min(which(idx>.9)))
> table(M)/length(X)*100
M
        0         1         2 
69.623397 25.512335  4.864268 
> image(seq(-180,180,length=nrow(X)),
+ seq(-90,90,length=ncol(X)),
+ M[,ncol(X):1],col=c("white",
+ "light green",col="red"),
+ axes=FALSE,xlab="",ylab="")

Supports de cours en R

Lundi, second cours de R pour la formation Actuaire: Data Science. Pour la première fois, je quitte mon format de slides (au format pdf) pour autre chose. Plus pratique pour présenter du code R. Une solution est de passer par R Markdown, et plus spécifiquement Slidify (qui permet de créer des slides)

On a un document, dans l’éditeur de RStudio, et la structure est assez simple. On a un préambule, et ensuite du texte,

C’est du Markdown, tout simplement. Mais on peut aussi mettre du code R, pour générer des tableaux, des sorties, des graphes,

Ensuite, au  lieu de faire run pour exécuter le code, on utilise Knit HTML, qui va permettre de générer les slides

Plus précisément, ces derniers vont s’ouvrir dans une fenêtre RStudio, mais qu’on peut aussi visualiser dans un browser,

C’est  avec ces fonctions que les slides ont été générés.

Une autre solution est de passer par Jupyter. L’installation prend un peu de temps (mais c’est relativement bien documenté). En lançant Jupyter, on crée un notebook en R,

Par défaut, on a un premier bloc de lignes de commandes

Pour en rajouter avant, ou après, rien de plus simple

On a alors un second blog, après l’affichage d’une sortie de régression, afin de générer un graphique

On peut soit garder les instructions, telles quelles, soit les exécuter, afin de visualiser les sorties

qui peuvent être des sorties de régression, ou des graphiques

On peut aussi insérer des blocs qui seront un Markdown (du texte, quoi), et pas seulement du code

On peut par exemple insérer un ‘titre’ de section

Et là encore, on exécute, pour afficher.

Cette année, pour mon cours de R, j’ai opté pour Slidify, mais peut être que Jupyter pourrait être une alternative intéressante, pour rendre le cours plus interactif. On en reparlera…

Prévision de l’intérêt pour le chocolat, avec des modèles ARIMA

Hier en cours, on s’était lancé sur la prévision de la série obtenue sur Google Trends pour le mot clé “chocolat“. La série est la suivante

> url = "http://freakonometrics.free.fr/report-chocolat.csv"
> report = read.csv(url,skip=6,
+ header=FALSE,nrows=636)
> X=report$V2
> plot(X,type="l")

Comme on l’avait évoqué, c’est plus simple de construire un modèle ARIMA sur des séries mensuelles qu’hebdomadaire. Donc on va mensualiser la série. Ah oui, on va aussi enlever les deux dernières années pour faire une comparaison de modèles (en confrontant les prévision avec les observations)

> 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")
> 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)
> base$AM = base$A*100+base$M
> Y = tapply(base$X,as.factor(base$AM),mean)
> Z = ts(as.numeric(Y[1:(146-24)]),
+ start=c(2004,1),frequency=12)
> plot(Z)

Compte tenu de l’amplitude croissante des cycles, il peut être intéressant, ici, de passer au logarithme de la série. Et comme on observe une rupture dans la tendance (linéaire au début, et relativement constante ensuite), on peut ajuste une régression linéaire avec une rupture (mais continue)

> plot(as.numeric(Z),type="l")
> X=log(as.numeric(Z))
> T=1:length(X)
> db=data.frame(X,T)
> trend=lm(X~T+I((T-80)*(T>80)),data=db)
> plot(X,type="l")
> lines(predict(trend),col="red")

C’est le résidu de cette série qu’on va chercher à modéliser,

> Y=residuals(trend)
> plot(Y,type="l")

Pour modéliser cette série, on peut commencer par regarder sa série d’autocorrélation

>  plot(acf(Y,lag=36),lwd=5)

et d’autocorrélation partielles,

> plot(pacf(Y,lag=36),lwd=5)

La série est stationnaire, mais il y a une composante cyclique forte (peut être une racine unité saisonnière ?). On peut tenter un premier modèle AR, ou ARMA (avec un AR les résidus ne sont pas un bruit blanc).

> model1 = arima(Y,order=c(12,0,12),
+ seasonal = list(order = c(0, 0, 0), period = 12))

La série des résidus est ici un bruit blanc

> E=residuals(model1)
> plot(acf(E,lag=36),lwd=5)

On peut alors utiliser ce modèle pour faire des prévisions, sur la série initiale (avant de prendre le logarithme, et d’ajuster une tendance)

> futur=length(X)+1:24
> Yp=predict(model1,n.ahead=24)$pred+
+    predict(trend,newdata=data.frame(T=futur))
> Xp=exp(as.numeric(
+   Yp+.5*predict(model1,n.ahead=24)$se^2))
> plot(as.numeric(Z),type="l",xlim=c(1,150))
> lines(futur,Xp,col="red")
> Y2=tapply(base$X,as.factor(base$AM),mean)
> obs_reel=as.numeric(Y2[(146-24)+1:24])
> lines(futur,obs_reel,col="blue")    

avec en rouge notre modèle, et en bleu les vraies observations. On peut alors calculer la somme des carrés des erreurs, sur ces 24 observations,

> sum( (obs_reel-Xp)^2 )
[1] 190.9722

Mais on peut tenter d’autres modèles, par exemple en changeant la tendance (on ne va pas le faire) ou en changeant le modèle ARIMA (en intégrant une racine unité saisonnière)

> model3=arima(Y,order=c(0,0,2),
+ seasonal = list(order = c(0, 1, 1), period = 12))
> E=residuals(model3)
> model3

Coefficients:
         ma1     ma2     sma1
      0.2246  0.3034  -0.9999
s.e.  0.0902  0.0925   0.3503

sigma^2 estimated as 0.002842:  log likelihood = 152.37,  aic = -296.75

On vérifie que les résidus sont effectivement un bruit blanc

> plot(acf(E,lag=36),lwd=5)
> Box.test(E,lag=12)

	Box-Pierce test

data:  E
X-squared = 6.326, df = 12, p-value = 0.8988

On peut alors faire une prévision sur notre série brute,

> futur=length(X)+1:24
> Yp=predict(model3,n.ahead=24)$pred+
+   predict(trend,newdata=data.frame(T=futur))
> Xp=exp(as.numeric(
+   Yp+.5*predict(model3,n.ahead=24)$se^2))
> plot(as.numeric(Z),type="l",xlim=c(1,150))
> lines(futur,as.numeric(Xp),col="red")
> Y2=tapply(base$X,as.factor(base$AM),mean)
> obs_reel=as.numeric(Y2[(146-24)+1:24])
> lines(futur,obs_reel,col="blue")

La somme des carrés des erreurs est ici un peu plus faible

> sum( (obs_reel-Xp)^2 )
[1] 173.8138

C’est à dire que ce second modèle est meilleur que le précédant, sur les deux dernières années. C’est probablement celui que l’on devrait utiliser pour faire des prévisions pour les années à venir (à moins qu’en continuant l’analyse on n’en trouve un meilleur).

R Crash Course, Data Science for Actuaries, Year 2

This Monday, we will start the second year of the Actuary: Data Science (ADS) program, supported by the (French) Institute of Actuaries. I will be there on monday morning for the opening, and we will start the R & Datamining course. The slides are now online,

In order to get nice slides, I have been using slidify. The R code is available online, as well as some pdf version of the slides.

Dynamique de la Pyramide des Ages

Très joli billet sur blog.revolutionanalytics.com avec un code de @kyle_e_walker permettant, très simplement (moyennant une inscription pour avoir une clé permettant d’utiliser l’API du census) de construire une pyramide des âges dynamiques.

> devtools::install_github('walkerke/idbr')
> library(idbr)
> library(ggplot2)
> library(animation)
> library(dplyr)
> library(ggthemes)
> idb_api_key("mykey1239F2f324zf9GGZgege32R2ii4")

On importe alors les données pour les hommes et les femmes,

> male <- idb1('FR', 2010:2050, sex = 'male') %>%
+   mutate(POP = POP * -1,
+   SEX = 'Male')
> female <- idb1('FR', 2010:2050, sex = 'female') %>% mutate(SEX = 'Female')

et on stocke le tout

> france <- rbind(male, female) %>%
+   mutate(abs_pop = abs(POP))

Ensuite, on crée l’animation,

> saveGIF({
+   
+   for (i in 2010:2050) {
+     
+     title <- as.character(i)
+     
+     year_data <- filter(france, time == i)
+     
+     g1 <- ggplot(year_data, aes(x = AGE, y = POP, fill = SEX, width = 1)) +
+       coord_fixed() + 
+       coord_flip() +
+       annotate('text', x = 98, y = -800000, 
+       label = 'Data: US Census Bureau IDB; idbr R package', size = 3) + 
+       geom_bar(data = subset(year_data, SEX == "Female"), stat = "identity") +
+       geom_bar(data = subset(year_data, SEX == "Male"), stat = "identity") +
+       scale_y_continuous(breaks = seq(-1000000, 1000000, 500000),
+       labels = paste0(as.character(c(seq(1, 0, -0.5), c(0.5, 1))), "m"), 
+       limits = c(min(france$POP), max(france$POP))) +
+       theme_economist(base_size = 14) + 
+       scale_fill_manual(values = c('#ff9896', '#d62728')) + 
+       ggtitle(paste0('Population structure of France, ', title)) + 
+       ylab('Population') + 
+       xlab('Age') + 
+       theme(legend.position = "bottom", legend.title = element_blank()) + 
+       guides(fill = guide_legend(reverse = TRUE))
+     print(g1) 
+   }
+ }, movie.name = 'france_pyramid.gif', interval = 0.1, ani.width = 700, ani.height = 600)

Et le résultat est vraiment joli, non ?