Tag Archives: computer

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 ?

Taux d’intérêt négatifs et explosion des valeurs de rentes

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrage évaluation du préjudice corporel, de Max le Roy. Comme l’indiquait le descriptif du livre (à l’époque) j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs,

“La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

En 2011, pour une mise à jour du livre, j’avais refait les calculs avec des tables mises à jours. Oui, à l’époque je trouvais que la priorité était d’arrêter d’utiliser les tables TD et TV qui dataient de 1988. Mais on était resté avec des taux à 5%. J’avais d’ailleurs été surpris (je peux le dire maintenant, il y a prescription) lorsque j’avais tenté de discuter le choix du taux d’intérêt, et qu’il m’avait été répondu “ah, ça change quelque chose ?“.

Oui, les taux d’actualisation impactent les valeurs des rentes. Et je n’avais pas idée à quel point. Le code utilisé pour construire les valeurs de rentes est assez simple.

> PRIX=function(annee=2011,age,sex="HOM",
+ taux=0.04, duree,C=1000){
+  file  = paste("http://freakonometrics.free.fr/",
+  sex,"-table-SPLx.csv",sep="")
+  B  = read.table(file,header=TRUE,sep=",")
+  an    = annee-age; if(an>2005){an=2005}
+  nom   = paste("X",an,sep="")
+  L     = B[,nom]
+  Q     = L[(age+1):length(L)]/L[(age+1)]
+  actualisation = (1+taux)^(0:min(duree,120-age))
+  prixsup = sum(Q[2:(min(duree,120-age)+1)]/
+  actualisation[2:(min(duree,120-age)+1)] )
+   prixinf = sum(Q[1:(min(duree,120-age))]/
+  actualisation[1:(min(duree,120-age))] )
+  return(C*c(prixsup,prixinf))}

Je renvois ici deux “prix” de rente, suivant que la rente est versée en fin ou en début d’année. Par exemple, avec un taux d’actualisation de 5% (comme dans le livre), pour un homme de 50 ans, si on considère une rente versée jusqu’à sa mort (j’ai mis ici une durée maximal de 150 années), la valeur actualisée espérée d’une rente (en tenant compte des probabilités de décès) de 1000 euros par an (encore une fois, tant que la personne est en vie) est de l’ordre de 17,000 euros.

> PRIX(age=50,sex="HOM",duree=150,taux=.05)
[1] 16435.31 17435.30

Avec un taux d’actualisation de 0%, on atteinte plus du double.

> PRIX(age=50,sex="HOM",duree=150,taux=0)
[1] 38822.48 39822.46

Et c’est encore pire pour les jeunes. Si on calcule pour plusieurs âges des rentes versées jusqu’au décès, avec différents niveaux de taux (5%, 2%, 0% et -2%, oui, je prends un taux négatif)

> vage=seq(15,65,by=5)
> vp5=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.05)[1])(vage)
> vp2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0.02)[1])(vage)
> vp0=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=0)[1])(vage)
> vpn2=Vectorize(function(a) PRIX(age=a,sex="HOM",duree=150,taux=-0.02)[1])(vage)
> plot(vage,vpn2,type="b",col="red")
> lines(vage,vp0,type="b")
> lines(vage,vp2,type="b",col="blue")
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"))

on observe une explosion des rentes pour les jeunes (ce qui n’est pas non plus surprenant)

Pour un jeune de 20 ans, la valeur actualisée probable d’une rente est presque 10 fois plus élevée, si on passe d’un taux de 5% à -2%

> PRIX(age=15,sex="HOM",duree=150,taux=.05)
[1] 19427.94 20427.93
> PRIX(age=15,sex="HOM",duree=150,taux=-.02)
[1] 200778.3 201770.3

Ce qu’on peut visualiser sur le graphique ci-dessous (avec, par âge, 100 qui correspond à la valeur obtenue avec un taux de 5%)

> plot(vage,vpn2/vp5*100,type="b",col="red")
> lines(vage,vp0/vp5*100,type="b")
> lines(vage,vp2/vp5*100,type="b",col="blue")
> abline(h=100,lty=2)
> legend("topright",c("taux -2%","taux 0%","taux 2%"),col=c("red","black","blue"),pch=1)

Passer de 5% à 2% donnait des rentes plus important, en passant de 100 à 200. En passant à 0% on passe de 100 à 400 pour les jeunes. Mais en passant à des taux négatifs, on passe de 100 à 1000 !

Dans un scénarios où les taux restent encore bas longtemps, on peut imaginer que les indemnités versées aux victimes de préjudices corporels vont exploser dans les mois à venir…

Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that  has a similar distribution… so there is no reason to test for a Gamma model for  before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if  has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution  is Gaussian but also the non-conditional distribution of .

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of  is also Gaussian

> library(nortest)
> ad.test(Y)

	Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

	Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here  is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value}
 
> plot(2:100,Vectorize(pv)(2:100),type="l")
> abline(h=.05,col="red")

So here, it could be possible to have also a Gaussian distribution, for . As least to accept that assumption, statistically.

In the context of a Poisson regression, it is well know that it’s not possible to have at the same time  that is Poisson distributed (that’s a Poisson regression) and also  that is Poisson distributed. That simply comes from the fact that

while

and because of the conditional Poisson distribution, then

Thus,

So  cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for .

More generally, it is very difficult to have a distribution family for   that is also the distribution of the non-conditional variable . In the context of a finite mixture ( takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference).

Hence, as a keep saying, before running a regression model on with some given family, it is never a good idea to check if the non-conditional distribution  has the same distribution. Because there is no reason, usually, to remain in the same family.