Tag Archives: computer

Unicode, UTF-8… la base de l’encodage

Un petit billet rapide, sans aucune prétention. Juste pour revenir sur des trucs bizarres sur l’encodage… J’avais commencé à m’y intéresser il y a 3 ans, alors que je discutais avec des utilisateurs de R en Inde et au Japon. Mais revenons à la base… Si on regarde la page “unicode” de Wikipedia (ou plus précisément, si on regarde le code source, on peut avoir des informations (dans les métadonnées) : en l’occurrence, on voit que la page est encodée en UTF-8

C’est cohérent car de plus en plus de documents en ligne semblent encodé dans ce format. Si on revient à la source, l’encodage c’est lié au passage par les bytes 0/1. Plein de sites expliquent (en gros) la philosophie. Quand j’étais petit, on nous expliquait les caractères ASCII, qui utilisent 7 bits ce qui permet de disposer  128 caractères uniquement, numérotés de 0 à 127.  Genre “A” était le 65ème caractère, et “a” le 97ème. Il y avait des symboles exotiques (comme “@” qui est le 64ème) ou la ponctuation, mais pas d’accents, comme “à”. Si on veut les lettres avec des accents, il faut plus que 128 caractères. Est-alors arrivée ISO/CEI 8859, qui proposait d’encoder les caractères sur 8 bits (et pas 7). On parle aussi de ISO 8859-1 plus populaire sous le nom “latin-1” très populaires dans les pays latins, avec le “à“, et puis on a continué d’enrichir pour arriver au ISO 8859-15 (latin-9) avec dedans le symbole “€”. Pour les langues asiatiques, il a fallu enrichir encore. Bref, c’est un jeu sans fin. Dans les années 90, est arrivée la norme ISO 10646 correspondant au “jeu universel de caractères” ou UCS pour Universal Character Set (en anglais), qui donnera Unicode. Je crois que les puristes font la différence entre les jeux de caractères (l’ensemble de caractères auxquels on attribue à chacun un point de code unique, la liste des caractères ASCII par exemple) et l’encodage (qui est la façon de représenter ce point de code en mémoire). Mais on va faire simple…

Là où ça se complique c’est que les OS utilisaient (par défaut) des encodages différents. Sous Windows, UTF-16 et Unicode sont majoritairement utilisé. Sous Mac OS,il y avait un truc appelé MacRoman mais il semble que UTF-8 l’emporte aujourd’hui. Enfin, sous linux, historiquement c’était du latin-1 par défaut, mais progressivement l’UTF-8 semble prendre le dessus. Bref, tous ceux qui naviguent entre les OS le savent : les accents c’est l’enfer. Et le principal soucis c’est qu’il est impossible, quand on récupère un fichier d’avoir son encodage a priori. Regardons le fichier suivant BaseANSI.txt

On peut tenter de regarder s’il y a des sortes de méta-données dans le document…

mais ce n’est pas le cas. Le document est un document brut, avec juste les données, et aucune métadonnées. Rien sur l’encodage (en particulier, car c’est ce qui nous intéresse aujourd’hui). Pour l’anecdote, si au lieu d’ouvrir le fichier avec mon bloc-note sous Windows, on tente sur un mac, on a

Et quand je passe par le terminal DOS sous Windows, c’est la même chose (ce qui est un peu rassurant)

Regardons un autre fichier BaseUNICODE.txt. Si vous voyez une différence avec l’autre, dites le moi…

Pourtant il y en a une : la preuve, dans un terminal, on peut lire les accents du fichier,

Un dernier ? Regardons BaseUTF8.txt,Là encore, on a l’impression d’avoir le même fichier. Et pourtant non.

Commençons par récupérer les trois fichiers (via R)

download.file(url="http://freakonometrics.free.fr/BaseUTF8.txt",
destfile="BaseUTF8.txt")
download.file(url="http://freakonometrics.free.fr/BaseANSI.txt",
destfile="BaseANSI.txt")
download.file(url="http://freakonometrics.free.fr/BaseUNICODE.txt",
destfile="BaseUNICODE.txt")

Techniquement, ces documents sont différents, même si dans le bloc note, on ne voyait rien : il n’ont pas la même taille !

file.size("BaseUNICODE.txt")
[1] 322
file.size("BaseANSI.txt")
[1] 163
file.size("BaseUTF8.txt")
[1] 171

On peut essayer de les lire, pour voir…

read.table("BaseANSI.txt",header=TRUE,sep=";")
                       Département No Habitants
1                              Ain 01    631877
2                          Ardèche 07    324209
3                 Bouches-du-Rhône 13   2016622
4                     Corse-du-Sud 2A    152730
5 Côte-dOr;21;533147\nCôtes-dArmor 22    598357

On a un petit soucis avec l’apostrophe, mais rien de grave

read.table("BaseANSI.txt",header=TRUE,sep=";",quote="\"")
       Département No Habitants
1              Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4     Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

R arrive parfaitement a lire (sur mon ordinateur du moins) cette première base de données. Tentons une seconde

read.table("BaseUTF8.txt",header=TRUE,sep=";",quote="\"")
    ï..DÃ.partement No Habitants
1               Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4      Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

On a manifestement un soucis, un soucis d’encodage. Le soucis c’est que, d’ordinaire, on ne sait pas comment le fichier a été encodé. Et comme on l’a vu, il n’y a pas de métadonnées qui nous donne l’encodage. Comment faire ? Heureusement, il existe une fonction sous R qui devine l’encodage. Littéralement.

library(readr)
guess_encoding("BaseUTF8.txt", n_max = 1000)
# A tibble: 3 x 2
  encoding   confidence
 
1 UTF-8           1.00 
2 ISO-8859-1      0.620
3 ISO-8859-2      0.430

Autrement dit, le best-guess est ici un encodage UTF-8. Essayons pour voir,

read.table("BaseUTF8.txt",header=TRUE,sep=";",encoding="UTF-8",quote="\"")
  X.U.FEFF.Département No Habitants
1                  Ain 01    631877
2              Ardèche 07    324209
3     Bouches-du-Rhône 13   2016622
4         Corse-du-Sud 2A    152730
5            Côte-d'Or 21    533147
6        Côtes-d'Armor 22    598357

Ça marche ! ou disons que ça marche presque. J’ai un soucis dans mes libellés de variables. Mais disons que c’est anecdotique ici.

Tentons le dernier fichier

read.table("BaseUNICODE.txt",header=TRUE,sep=";",quote="\"")
  ÿþD
1  NA
2  NA
3  NA
4  NA
5  NA
6  NA
7  NA
8  NA

Comme ça ne marche pas, essayons d’avoir un best-guess sur l’encodage

guess_encoding("BaseUNICODE.txt", n_max = 1000)
# A tibble: 3 x 2
  encoding   confidence
 
1 UTF-16LE        1.00 
2 ISO-8859-1      0.530
3 ISO-8859-2      0.270

On en tient un !

read.table("BaseUNICODE.txt",header=TRUE,encoding="UTF-16LE")
  ÿþD
1  \n
2   B
3  \n
4   C
5  \n
6   C
7  \n
8   C

Damned, ça ne marche pas… Heureusement, Ewen est venu à ma rescousse. La solution est relativement étrange

read.table("BaseUNICODE.txt",header=TRUE,fileEncoding="UTF-16LE", sep=";", quote="")
       Département No Habitants
1              Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4     Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

Oui, on n’utilise pas l’encodage via “encoding” mais “fileEncoding” ! Subtil non… Cela dit, si on voulait être cohérent on devrait utiliser la fonction du package pour lire la base, non ? Car dans le package “readr” il y a une fonction “read_table”

read_table("BaseUNICODE.txt",locale=locale(encoding = "UTF-16LE"))
Error in guess_types_(datasource, tokenizer, locale, n = guess_max) : 
  Incomplete multibyte sequence

Sauf que quand on la lance, on a un message d’erreur… Bref, l’encodage c’est compliqué. C’est ce qui est expliqué dans le forum de discussion du package cela dit. Je ne suis pas sûr que mon billet serve à quoi que ce soit, mais je voulais garder une trace de tout ça !

Modéliser la consommation électrique

A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à ca

plot(electricite[passe,"Load"],type="l")

On peut tenter un modèle assez simple, où la consommation à la date Y_t est fonction d’une tendance linéaire a+bt, de la position dans l’année (sous une forme non-paramétrique), de la température T_t (sous une forme non-linéaire, disons un polynôme de degré 3), de la températeure en t-1, T_{t-1} et de l’indice de production industriel IPI_t (sous une forme linéaire pour faire simple).

model1 = lm(Load~1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

L’idée de la fonction polynomial pour l’impact de la température venait du graphique suivant (où on représente la série de la consommation, une fois enlevée la tendance linéaire)

On peut aussi supposer une forme autorégressive, où Y_{t} sera fonction de Y_{t-1}

model2 = lm(Load~1+Load1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

On peut alors tenter de faire de la prévision. Le soucis avec ce second modèle, c’est la composante autorégressive. Pour prévoir Y_{t+h}, il faut utiliser la prévision faite en t+h-1, \hat Y_{t+h-1}.

base_prevision=data.frame(
Load_prev = NA,  Load1=c(electricite[length(passe),"Load"],rep(NA,110),
Time = electricite[futur,"Time"],
NumWeek = electricite[futur,"NumWeek"],
Temp = MODELES[,"y4"], 
Temp1 = c(electricite[length(passe),"Temp"],
MODELES[1:110,"y4"]),
IPI = electricite[futur,"IPI"])
for(t in 1:110){
 p = predict(model2,newdata=base_prevision)[t]
 base_prevision[t,"Load_prev"] = p
 base_prevision[t+1,"Load1"] = p}

On peut alors comparer la prévision \hat Y_{t} avec l’observation Y_{t}.

plot(1:111,electricite[futur,"Load"],type="l")
lines(1:111,base_prevision[,"Load_prev"],col="red")

On est bon les mois d’été (on capture bien le pic de la première quinzaine du mois d’août) par contre, on sous-estime fortement la consommation l’hiver. Et c’est gênant…

Pour finir, on peut oublier le passage par les variables explicatives, et tenter directement un modèle de série temporelles.

plot(electricite[passe,"Load"],type="l")

Le soucis, c’est le côté hétéroscédastique de la série, avec une pente sur le minimum plus faible que sur le maximum.

n=length(passe)
plot(1:n,electricite[passe,"Load"],type="l")
m=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),min)
M=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),max)
vm=which(electricite[passe,"Load"]%in%m$x)
vM=which(electricite[passe,"Load"]%in%M$x)
points(sort((1:n)[vm]),m$x,pch=19,col="blue")
points(sort((1:n)[vM]),M$x,pch=19,col="blue")
xm=((1:n)[vm])
xM=((1:n)[vM])
regm=lm(m$x~xm,col="blue")
regM=lm(M$x~xM,col="blue")
abline(regm,lty=2)
abline(regM,lty=2)

Une solution classique (et simple) est de passe au logarithme

plot(electricite[passe,"Load"],type="l",log="y")
y=electricite[passe,"Load"]
z = log(y)
plot(z,type="l")
B = data.frame(z=z,t=1:length(z))
abline(lm(z~t,data=B),lty=2,col="red")

Il faut ensuite enlever la tendance linéaire, pour stationnariser la série

z = residuals(lm(z~t,data=B))
Z = ts(z,start=1996,frequency=52)
modelz1 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
modelz1 = arima(Z,order = c(4,0,0), 
seasonal = list(order = c(1,0,1)))

Ce premier modèle est stationnaire, sans racine unité. On peut tenter d’introduire une racine unité saisonnière

modelz2 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(0,1,0)))
modelz2 = arima(Z,order = c(1,0,3), 
seasonal = list(order = c(1,1,0)))

Enfin, le dernier sera un peu plus simple

modelz3 = arima(Z,order = c(1,0,0), 
seasonal = list(order = c(2,0,0)))

On stocke alors toutes les prévisions dans une base

DONNEES_u = data.frame(
  u1 = predict(modelz1,n.ahead = 111)$pred,
  u2 = predict(modelz2,n.ahead = 111)$pred,
  u3 = predict(modelz3,n.ahead = 111)$pred)

On rajoute alors la tendance linéaire à la prévision du résidu

reg = lm(z~t,data=B)
dfutur = data.frame(t=futur)
DONNEES_z = data.frame(
  z1 = DONNEES_u$u1+predict(reg,newdata=dfutur),
  z2 = DONNEES_u$u2+predict(reg,newdata=dfutur),
  z3 = DONNEES_u$u3+predict(reg,newdata=dfutur)

On a besoin de la variance pour la prévision. En effet, classiquement, la prévision sera l’espérance de notre grandeur d’intérêt. Ici, on a fait un modèle linéaire sur \log Y, i.e. \log Y \sim\mathcal{N}(\mu,\sigma^2) et donc \mathbb{E}[Y]=\exp\left(\mu+\sigma^2/2\right)

sigma=summary(reg)$sigma
DONNEES_seu = data.frame(
  seu1=sqrt(predict(modelz1,n.ahead = 111)$se^2+sigma^2),
  seu2=sqrt(predict(modelz2,n.ahead = 111)$se^2+sigma^2),
  seu3=sqrt(predict(modelz3,n.ahead = 111)$se^2+sigma^2))

On fait ici l’hypothèse que les estimateurs des prévisions des deux modèles (la tendance linéaire et celle du modèle autorégressif) sont indépendants et donc on peut sommer les termes de variance. Aussi, les prévisions pour Y sont

DONNEES_y = data.frame(
  y1 = exp(DONNEES_z$z1+1/2*DONNEES_seu$seu1^2),
  y2 = exp(DONNEES_z$z2+1/2*DONNEES_seu$seu2^2),
  y3 = exp(DONNEES_z$z3+1/2*DONNEES_seu$seu3^2))

On peut représenter les prévisions des trois modèles (avec les observations)

plot(futur,electricite[futur,"Load"],type="l")
lines(futur,DONNEES_y$y1,col="red")
lines(futur,DONNEES_y$y2,col="blue")
lines(futur,DONNEES_y$y3,col="purple")

On peut comparer avec la prévision précédente,

lines(futur,base_prevision[,"Load_prev"],
col="orange")

On y perd un peu l’été, par contre, on s’est un peu amélioré l’hiver.

Using convolutions (S3) vs distributions (S4)

Usually, to illustrate the difference between S3 and S4 classes in R, I mention glm (from base) and vglm (from VGAM) that provide similar outputs, but one is based on S3 codes, while the second one is based on S4 codes. Another way to illustrate is to manipulate distributions.

Consider the case where we want to sum (independent) random variables. For instance two lognormal distribution. Let us try to compute the median of the sum.

The distribution function of the sum of two independent (positive) random variables is F_{S_2}(x)=\int_0^x F_{X_1}(x-y)dF_{X_2}(x)

pSum2 = function(x) integrate(function(y) 
plnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value

Let us visualize that cumulative distribution function

vx=seq(0.1,50,by=.1)
vy=Vectorize(pSum2)(vx)
plot(vx,vy,type="l",ylim=c(0,1))
abline(h=.5,lty=2)

Let us find an upper bound to compute (in a decent time) quantiles

pSum2(350)
[1] 0.99195

and then use the uniroot function to inverse that function

qSum = function(u) uniroot(function(x) 
Vectorize(pSum2)(x)-u, interval=c(0,350))$root
vu=seq(.01,.99,by=.01)
vv=Vectorize(qSum)(vu)

The median is here

qSum(.5)
[1] 14.155

Why not consider the sum of three (independent) distributions ? Its cumulative distribution function can be writen using our previous function F_{S_3}(x)=\int_0^x F_{S_2}(x-y)dF_{X_3}(x)

pSum3 = function(x) integrate(function(y) 
pSum2(x-y)*dlnorm(y,2,2),0,x)$value

If we look at some values we good

pSum3(4)
[1] 0.015624
pSum3(5)
Error in integrate(function(y) plnorm(x - y, 1, 2) * 
dlnorm(y, 2, 1),  : 
  maximum number of subdivisions reached

So obviously, there are computational issues here.

Let us consider the following alternative expression F_{S_3}(x)=\int_0^x F_{X_3}(x-y)dF_{S_2}(x). Of course, it is necessary here to compute the density of the sum of two variables

dSum2 = function(x) integrate(function(y) 
dlnorm(x-y,1,2)*dlnorm(y,2,1),0,x)$value
pSum3 = function(x) integrate(function(y) 
dlnorm(x-y,2,2)*dSum2(y),0,x)$value

Again, let us compute some values

pSum3(4)
[1] 0.0090285
pSum3(5)
[1] 0.01186

This one seems to work quite well. But it is just an illusion.

pSum3(9)
Error in integrate(function(y) dlnorm(x - y, 1, 2) *
 dlnorm(y, 2, 1),  : 
  maximum number of subdivisions reached

Clearly, with those S3-type functions, it wlll be complicated to run computations with 3 variables, or more.

Let us consider distributions in the S4-type format of the following package

library(distr)
X1 = Lnorm(mean=1,sd=2)
X2 = Lnorm(mean=2,sd=1)
S2 = X1+X2

To compute the median, we simply have to use

distr::q(S2)(.5)
[1] 14.719

We can also visualize it easily

plot(q(S2))

which looks (very) close to what we got, manually.  But here, it is also possible to work with the sum of 3 (independent) random variables

X3 = Lnorm(mean=2,sd=2)
S3 = X1+X2+X3

To compute the median, use

distr::q(S3)(.5)
[1] 33.208

The function is here

plot(q(S3))

(Advanced) R Crash Course, for Actuaries

The fourth year of the Data Science for Actuaries program started this morning. I will be there for the introduction to R. The slides are available online (created with slidify, the .Rmd file is also available)

A (standard) markdown is also available (as well as the .Rmd file). I have to thank Ewen for his help on slidify (especially for the online quizz, and the integration of leaflet maps or the rgl animated graph….)

Justice “actuarielle”, algorithmes… et données

Il y a un peu plus d’un an, Virginie Gautron m’envoyait plein de documents sur la “justice actuarielle”, concept que je découvrais alors. Pour comprendre un peu de quoi il s’agit, je peux renvoyer vers pénalité et gestion des risques : vers une justice « actuarielle » en Europe ? qui dresse un état de l’art, en francais. J’avoue avoir un peu mis de côté ensuite (faute de temps), et puis l’autre jour, j’ai dévoré surveillé et punir de Nicolas Bourgoin, qui contient un chapitre sur le sujet.

Sans vouloir tomber dans le cliché des deux cultures (à la CP Snow), il faut bien admettre que pour un économètre, lire des écrits de juristes parlait de modèles prédictifs (ou actuariels) est une expérience marquante. Et ce qui me trouble (j’ai d’ailleurs le même soucis avec la psychologie – c’est un sujet proche) c’est que l’accent est souvent mis sur le modèle, l’algorithme, ce truc qui semble fasciner les juristes, alors qu’il me semble que le coeur du problème, ce sont les données. On nous dit, par exemple, que les algorithmes sont racistes (c’est un peu la thèse de Cathy O’Neil dans Weapons of Math Destruction), mais le problème ce n’est pas le modèle. C’est le fait qu’il y a un énorme biais de sélection dans les données ! Et selon l’adage “garbage in, garbage out“… dit de manière plus policée, je pourrais tenter un “purgamentum init, exit purgamentum” pour parler comme un juriste. Ou dit autrement, si on entre des mauvaises données en entrée, faut pas espérer grand chose en sortie…

Pour comprendre le biais de sélection, on peut relire Mostly Harmless Econometrics de Joshua Angrist et Jörn-Steffen Pischke. En particulier, dans l’introduction, ils posent la question simple suivant : vaut-il mieux, suite à une opération, renvoyer les gens chez eux, ou bien les laisser en chambre de convalescence, à l’hôpital ? Le soucis est que les données sont biaisées : oui, si on regarde les données brutes, ceux qui sont rentré chez eux se sont remis plus vite… mais peut être a-t-on autorisé les gens à rentrer justement parce qu’ils allaient bien ? Et on retrouve ce biais de sélection un peu partout. Probablement même encore davantage dans les problèmes relatifs à la justice.

Prenons un exemple simple. On va s’intéresser à la probabilité d’être coupable d’un crime, et d’être condamné par la justice. C’est ma variable Y. Pour avoir un modèle économétrique, il me faut des variables explicatives. J’ai une variable X_3 et une autre X_4, la première influençant vraiment la probabilité de commettre un crime, et l’autre étant juste du bruit. Maintenant, entrons dans le coeur du modèle. Supposons que j’ai une variable corrélée à ma richesse \tilde{X}_2, comme le lieu d’habitation, un critère racial, ou pour faire simple le diplôme. Notons X_1 le niveau de diplôme de la personne (mais encore une fois, on peut imaginer n’importe quoi, l’idée est que cette variable est corrélée à la richesse). Et surtout, cette variable n’est pas du tout une variable causale dans mon modèle. La probabilité de commettre un crime n’est pas influencée par X_1. Une personne ayant un diplôme de second cycle universitaire a autant de chance de commettre un crime qu’une personne ayant raté le brevet. On va même supposer que la probabilité de commettre un crime soit indépendante de la richesse. En fait, la richesse me permet de payer un bon avocat. Et avoir un bon avocat permet de diminuer la probabilité d’être condamné. Pour faire simple, j’ai deux niveaux de richesse, qui me donne soit un avocat commis d’office, soit un avocat qui s’impliquera peut être davantage (je simplifie à des fins pédagogiques). Bref, mon modèle aura la forme suivante : la probabilité d’être condamné pour un crime s’écrit p=(0.025+0.05\cdot X_3)\cdot[k\cdot \mathbb{1}(X_2=R)+\mathbb{1}(X_2=P)] Autrement dit, ma probabilité est une fonction croissante (linéaire) de X_3, et si je suis riche (ce qu’on note X_2=R), j’ai moins de chance d’être condamné, avec un facteur k\in[0,1], toutes choses étant égales par ailleurs, ceteris paribus, histoire de parler latin comme les juristes.

Pour simuler un jeu de données, on considère une population d’un million de personnes, et je suppose que ma richesse est très corrélée avec mon niveau de diplôme. Pour info, la classe A ce sont les gens les plus diplômés, et je suppose qu’un bon avocat divise par trois la chance d’être condamné.

 n=1000000
 r=.95
 k=1/3
 S=matrix(c(1,r,r,1),2,2)
 library(mnormt)
 set.seed(1)
 vectX=rmnorm(n,varcov=S)
 X1=cut(vectX[,1],breaks=c(-10,1,-1,10), c("C","B","A"))
 X1=relevel(X1,"A")
 X2=cut(vectX[,2],breaks=c(-10,1,10),c("P","R"))
 X3=runif(n)
 X4=runif(n)
 Y1=rbinom(n,size = 1,prob=(.025+.05*X3)*(k*(X2=="R")+(X2=="P")))
 B0=data.frame(X1,X2,X3,X4,Y1)

On peut faire un peu de statistiques pour regarder les caractéristiques des gens condamnés (uniquement)

B1=data.frame(X1=X1[Y1==1],X2=X2[Y1==1],X3=X3[Y1==1],X4=X4[Y1==1])

Par exemple si je regarde la proportion de gens les plus diplômés: ils sont 15% dans la population totale, et 5% parmi les personnes condamnées

 mean(B0$X1=="A")
[1] 0.158734
 mean(B1$X1=="A")
[1] 0.05016072
 t.test(B0$X1=="A",B1$X1=="A")
data:  B0$X1 == "A" and B1$X1 == "A"
t = 219.5, df = 1283500, p-value < 2.2e-16
95 percent confidence interval:
 0.1076038 0.1095428
sample estimates:
 mean of x  mean of y 
0.15873400 0.05016072

Autrement dit, on a une différence significative entre ces deux pourcentages. Pareil si je regarde la valeur moyenne de X_3, qui était de 50% dans la population et de 58% dans la population condamnée

 t.test(B0$X3,B1$X3)
data:  B0$X3 and B1$X3
t = -163.46, df = 844940, p-value < 2.2e-16
95 percent confidence interval:
 -0.08455430 -0.08255058
sample estimates:
mean of x mean of y 
0.5000094 0.5835618

alors que pour X_4, on a 50% dans la population totale, et 50% dans la population condamnée. Mais on s’y attendait….

 t.test(B0$X4,B1$X4)
data:  B0$X4 and B1$X4
t = -0.81161, df = 48645, p-value = 0.417
95 percent confidence interval:
 -0.003869879  0.001603452
sample estimates:
mean of x mean of y 
0.5001326 0.5012659

Passons à la seconde phase. Il semble y avoir un vrai intérêt dans la littérature sur les modèles actuariels pour la construction de modèles liés à la probabilité de récidive. Là encore, je vais faire une hypothèse forte : faire de la prison n’influence en rien la probabilité de récidiver. Autrement dit, pour cette population déjà condamnée, je vais refaire tourner mes simulations, avec une probabilité de commettre un crime qui dépend seulement de X_3, mais avec toujours notre effet richesse qui permet de payer un avocat qui fera baisser la probabilité d’être condamné.

 set.seed(2)
 n1=nrow(B1)
 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))
 B2=cbind(B1,Y2)
 reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

J’ai construit ici mon modèle permettant de calculer une probabilité de récidive. Le modèle est construit sur le niveau de diplôme X_1 (qui n’influence en rien, encore une fois, la probabilité de commettre un crime), la variable X_3 (qui effectivement joue un rôle) et la variable de bruit X_4.

 set.seed(2)
 n1=nrow(B1)
 Y2=rbinom(n1,size = 1,prob=(.025+.05*B1$X3)*(k*(B1$X2=="R")+(B1$X2=="P")))
 B2=cbind(B1,Y2)
 reg=glm(Y2~X1+X3+X4,data=B2,family=binomial)

Je fais ici un score logistique car c’est le plus logique.

 summary(reg)
Coefficients:
            Estimate Std. Error z value  Pr(|z|)    
Inter.) -4.04753    0.11592 -34.916    2e-16 ***
X1C      0.61311    0.10914   5.618 1.93e-08 ***
X1B      0.68807    0.09977   6.896 5.34e-12 ***
X3       1.00610    0.08052  12.495    2e-16 ***
X4      -0.13562    0.07286  -1.861   0.0627 .

Premier élément à observer, la variable X_4 est à la limite d’être significative… Certains auteurs sur ces thèmes diraient que c’est un signal faible. Pourtant, en l’occurence, c’est du bruit ! Juste du bruit !! Et parmi les signaux forts, très clairement significatifs, on retrouve X_1, notre niveau de diplôme. Ici, si on est faiblement qualifié (niveaux B ou C) alors on a une chance de récidive beaucoup plus importante ! C’est très fortement significatif ! Alors que pour rappel, la probabilé de commettre un crime, et de récidiver n’est aucunement causé par cette variable ! C’est juste une variable corrélée avec la variable richesse, qui permet de faire baisser la probabilité d’être condamné.

Plus intéressant, si je regarde la valeur du coefficient, en fonction de k, autrement dit la valeur d’un bon avocat (que l’on peut s’offrir si on est riche) on obtient la courbe suivante.

Autrement, sur la partie droite, même si un bon avocat ne baisse que de 5% la chance d’être condamné, le niveau de diplôme semble être significativement liée à la probabilté de récidive. Avec un pic autour de 30%. Si les avocats ont un pouvoir incroyable, par exemple en divisant par 20 la probabilité d’être condamné, étrangement, la variable de diplôme redevient très peu significative, mais je pense que c’est parce qu’on n’a alors plus personne de très diplômé dans la base.

On pourrait me reprocher d’avoir une corrélation trop forte entre richesse et niveau de diplôme, ici. En mettant une corrélation à 70% (dans mon modèle Gaussien latent), j’obtiens

qui donne un effet moindre, mais là encore, si un bon avocat permet de réduire la probabilité d’être condamné par un facteur compris entre 5% et 90%, alors la variable de diplôme semblera significative.

Mon modèle n’est pas faux ici. C’est essentiellement que j’ai un biais dans mes données. Et ne pas en tenir compte fait que le modèle n’a aucun sens, et ne peut pas être interprété sans commettre d’erreurs (et ici je cumule les âneries, en voyant un signal faible alors qu’il y a du bruit, et en voyant un signal fort pour une variable qui n’a aucun lien causal). J’insiste lourdement, mais le diplôme n’intervient aucunement ici. Et on peut remplacer le diplôme par toute variable corrélée à la richesse, comme locataire et propriétaire, etc.

Bref, l’utilisation de modèles économétrique dans des problèmes judiciaires n’est pas simple ! Je parle ici de biais de sélection, mais la formulation qui me semble la plus juste serait une mauvaise spécification du modèle (on n’utilise pas les bonnes variables, mais ces dernières sont peut-être non-observables – ou difficiles à justifier éthiquement : peut-on demander le revenu de la personne pour faire un score de récidive, et avoir un score qui dépende effectivement du revenu ?). Il y a aussi un problème d’identification, car on a du mal à savoir ce qui influence le fait de commettre un crime, et la condamnation. Ce qu’on observe, ce sont les condamnations, pas les actes (je ne parle pas des condamnations à tort ici, que j’avais pu évoquer dans un autre billet, sur le principe de précaution). Mais j’essayerais de revenir sur ce point après les vacances….

Dernière anecdote pour finir. Après avoir découvert ces modèles via Virginie Gautron, j’avais discuté avec une amie, présidente de tribunal, qui me parlait de la nullité des experts qu’elle avait sous la main. Et elle voulait faire un peu de ménage. Bref, elle me demandait si je pouvais faire un modèle attribuant un score à un expert. Et là encore, ma réponse a été que oui, faire un modèle c’est facile, à condition d’avoir des données. Comment puis-je avoir un score Y d’un expert? Qui sait ce qu’est un bon expert? C’est quelqu’un contre qui on ne demande pas de contre-expertise? Malheureusement, la justice ne fonctionne pas ainsi, et la demande n’a rien à voir avec la qualité de l’expert. Et la qualité, vu du point de vue de la juge, est une notion qui peut être très différente de celle des deux parties qui s’opposent… Bref, je pourrais facilement faire un modèle. Mais en avoir un qui ait du sens, c’est une autre paire de manches !

Optimal Portfolios Exam

Pour l’examen du cours de Gestion d’Actifs (optimal portfolios) un cas pratique avec des vraies données, disponibles sur

loc="http://freakonometrics.free.fr/base_portefeuille.csv"
base=read.csv2(loc,sep=";")
base[,1]=as.Date(as.character(base[,1]),"%Y-%m-%d")
for(i in 2:10) base[,i]=as.numeric(as.character(base[,i]))

Pour le rendu, je demande un markdown dont la trame est en ligne ici,  exam_portfolio.Rmd, et la sortie est exam_portfolio.html. Il faut télécharger le fichier Rmd, l’ouvrir avec RStudio, puis me renvoyer le fichier rempli (ainsi que le fichier html compilé) pour midi, à arthur.charpentier@univ-rennes1.fr.

The myth of interpretability of econometric models

There are important discussions nowadays about data modeling, to choose between the “two cultures” (as mentioned in Breiman (2001)), i.e. either econometrics models or machine/statistical learning models. We did discuss this issue recently in Econométrie et Machine Learning (so far only in French) with Emmanuel Flachaire and Antoine Ly. One argument often used by econometricians is the interpretability of econometric models. Or at least the attempt to get an interpretable model.

We also have this discussion in actuarial science, for instance in ratemaking (or insurance pricing). Machine learning based models usually perform better (for some a priori chosen metric), but actuaries claim that econometric models are more easily interpretable. In actuarial literature, we assume that claim frequency Y is driven by some non-observable risk factor \Theta, and therefore, we do have heterogeneous risks in our portfolio. And, it can be seen as legitimate to differentiate prices. Assume that this risk factor \Theta is strongly correlated with X_1, the age of the driver. Because in our portfolio, old drivers tend to have more accidents. Here, we could pretend to have a “causal story” (as defined in Freedman (2009)) because of a possible interpretation of the model. So it is natural here to consider a regression model of Y on X_1 to derive our actuarial pricing model. But assume that, possibly, risk factor \Theta is also strongly correlated with X_2, that can be related to spatial features (say latitude, which denoted a north/south position). Because in our portfolio, drivers living in the south tend to have more accidents (reads are known to be more dangerous there). Here, we could pretend to have a second “causal story”.

Of course, since \Theta is strongly correlated with X_1 and X_2, it means that X_1 and X_2 are strongly correlated. Here also, this correlation can be interpreted (not in a causal way as previously, but still), since we know that old people like to live in southern regions. So, what should we do here ? Let us run some simulations to  illustrate.

 set.seed(123)
 n=1e5
 Theta=rnorm(n)
 X1=Theta+rnorm(n)/8
 X2=Theta+rnorm(n)/8
 L=exp(-3+Theta)
 Y=rpois(n,L)
 B=data.frame(Y,X1,X2)

Our first idea was to consider a model where Y is “explained” by the first variable X_1,

 g1=glm(Y~X1,data=B,family=poisson)
 summary(g1)
 
Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
(Inter.) -2.97778    0.01544 -192.88   <2e-16 ***
X1        0.97926    0.01092   89.64   <2e-16 ***

As expected, our variable is “significant”, but also, probably more interesting, X_2, has no impact on the residuals

 B$e1=residuals(g1,type="pearson")
 g1e=lm(e1~X2,data=B)
 summary(g1e)
 
Coefficients:
          Estimate Std. Error t value Pr(>|t|)
(Inter.) 0.0003618  0.0031696   0.114    0.909
X2       0.0028601  0.0031467   0.909    0.363

The interpretation is that once we corrected claim frequency for the age of the drivers, there is no spatial effect here. So, a good model should be based only on the age of the drivers.

But we can also consider the other story. We can consider a model where Y is “explained” by the second variable X_2,

 g2=glm(Y~X2,data=B,family=poisson)
summary(g2)
 
Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
(Inter.) -2.97724    0.01544 -192.81   <2e-16 ***
X2        0.97915    0.01093   89.56   <2e-16 ***

Here also we have a valid model, that can be interpreted, and here also X_1, has no impact on the residuals

 B$e2=residuals(g2,type="pearson")
 g2e=lm(e2~X1,data=B)
 summary(g2e)
 
Coefficients:
          Estimate Std. Error t value Pr(>|t|)
(Inter.) 0.0004863  0.0031733   0.153    0.878
X1       0.0027979  0.0031504   0.888    0.374

The story is similar here. If we correct from the spatial pattern, claims frequency does not depend on the age of the driver.

So, what should we do now? We do have two models, and each of them is as interpretable as the other one. Note that we can not use any statistical tool to distinguish the two: they are comparable

 AIC(g1)
[1] 51013.39
 AIC(g2)
[1] 51013.15

Why not incorporate the two explanatory variables X_1 and X_2, at the same time, in our regression model, and let “the model” decide what to do…?

 g=glm(Y~X1+X2,data=B,family=poisson)
 summary(g)
 
Coefficients:
         Estimate Std. Error  z value Pr(>|z|)    
(Inter.) -2.98132    0.01547 -192.723    2e-16 ***
X1        0.49310    0.06226    7.920 2.38e-15 ***
X2        0.49375    0.06225    7.931 2.17e-15 ***

It looks like we completely lost the interpretability of the model, since our two explanatory variables are (strongly) correlated. Actually, instead of saying “use one, and drop the other one (since it brings no further information)”, it says “use both, each one will explain half of the variable”. Strange interpretation, isn’t it?  So why not try some LASSO here?

library(glmnet)
fit=glmnet(x=as.matrix(B[,c("X1","X2")]), 
    y=B$Y,family="poisson")
plot(fit,xvar="lambda")

Here also, it says that we either keep both, or none. So it cannot be used for variable selection (which is an important motivation to use LASSO technique). So, what should be do if we several interpretable models, but no way to choose? Because usually, we claim that we prefer to use a model with an interpretation. But what should be done here?

Optimal Portfolios, or sort of…

Last week, we got our first class on portfolio optimization. We’ve seen Markowitz’s theory where expected returns and the covariance matrix are given,

> download.file(url="http://freakonometrics.free.fr/portfolio.r",destfile = "portfolio.r")
> source("portfolio.r")
> library(zoo)
> library(FRAPO)
> library(IntroCompFinR)
> library(rrcov)
> data( StockIndex )
> pzoo = zoo ( StockIndex , order.by = rownames ( StockIndex ) )
> rzoo = ( pzoo / lag ( pzoo , k = -1) - 1 ) * 100
> Moments <- function ( x , method = c ( "CovClassic" , "CovMcd" , "CovMest" , "CovMMest" , "CovMve" , "CovOgk" , "CovSde" , "CovSest" ) , ... ) {
method <- match.arg ( method )
ans <- do.call ( method , list ( x = x , ... ) ) + return ( getCov ( ans ) )} > covmat=Moments(as.matrix(rzoo),"CovClassic")
> (covmat=round(covmat,1))
SP500 N225 FTSE100 CAC40 GDAX HSI
SP500   17.8 12.7 13.8 17.8 19.5 18.9
N225    12.7 36.6 10.8 15.0 16.2 16.7
FTSE100 13.8 10.8 17.3 18.8 19.4 19.1
CAC40   17.8 15.0 18.8 30.9 29.9 22.8
GDAX    19.5 16.2 19.4 29.9 38.0 26.1
HSI     18.9 16.7 19.1 22.8 26.1 58.1
> er=apply(as.matrix(rzoo),2,mean)
> (er=round(er,1))
SP500 N225 FTSE100 CAC40 GDAX HSI
0.6 -0.2 0.4 0.5 0.8 1.0
> ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50)

We can now visualize the efficient frontier (and admissible portfolios) below

> u=c(12,ef$sd,12,12)
> v=c(5,ef$er,-1,5)
> plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5)
> points(sqrt(diag(covmat)),er,pch=19,col="blue")
> text(sqrt(diag(covmat)),er,names(er),pos=4, col="blue",cex=.6)
> polygon(u,v,border=NA,col=rgb(0,0,1,.3))

https://freakonometrics.hypotheses.org/files/2017/11/image-voronoi-post-026-1.png

That was the starting point of our class. We did also mention that something important was actually hard to visualize on that graph : the correlation between returns. It is not in the points (which are univariate, with expected return and standard deviation), but in the efficient frontier. For instance, here is our correlation matrix

> (cormat=covmat/(sqrt(diag(covmat) %*% t(diag(covmat)))))
SP500 N225 FTSE100 CAC40 GDAX HSI
SP500   1.00 0.50 0.79 0.76 0.75 0.59
N225    0.50 1.00 0.43 0.45 0.43 0.36
FTSE100 0.79 0.43 1.00 0.81 0.76 0.60
CAC40   0.76 0.45 0.81 1.00 0.87 0.54
GDAX    0.75 0.43 0.76 0.87 1.00 0.56
HSI     0.59 0.36 0.60 0.54 0.56 1.00

We can actually change the correlation between FT500 and FTSE100 (which is here .786)

courbe=function(r=.786){
R=cormat
R[1,3]=R[3,1]=r
covmat2=(sqrt(diag(covmat) %*% t(diag(covmat))))*R
ef <- efficient.frontier(er, covmat2, alpha.min=-2.5, alpha.max=2.5, nport=50)
plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return",
xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5)
points(sqrt(diag(covmat)),er,pch=19,col=c("blue","red")[c(2,1,2,1,1,1)])
text(sqrt(diag(covmat)),er,names(er),pos=4,col=c("blue","red")[c(2,1,2,1,1,1)],cex=.6)
polygon(u,v,border=NA,col=rgb(0,0,1,.3))
}

for instance, with a correlation of 0.6, we get the following efficient frontier

> courbe(.6)

and with a stronger correlation

> courbe(.9)

So clearly, correlation does matter. A lot. But more important, one should keep in mind that expected returns and covariances are not given, but estimated. Previously, we did use the standard estimator for the variance matrix. But another (more robust) estimator can be considered

covmat=Moments(as.matrix(rzoo),"CovSde")
er=apply(as.matrix(rzoo),2,mean)
ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50)
plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return",xlim=c(3.5,11),ylim=c(0,2.5),col="red",lwd=1.5)
points(sqrt(diag(covmat)),er,pch=19,col="blue")
text(sqrt(diag(covmat)),er,names(er),pos=4,col="blue",cex=.6)
polygon(u,v,border=NA,col=rgb(0,0,1,.3))

It did influence (horizontal) position of points, since variances are now different, as well as the efficient frontier, with clearly much lower variances that can be reached.

And to illustrate a last point, to illustrate the fact that we do have estimators based on observed returns, what if we had observed different ones? A way to get an idea of what might happened is to use bootstrap, e.g. of daily returns.

> covmat=Moments(as.matrix(rzoo),"CovClassic")
> er=apply(as.matrix(rzoo),2,mean)
> ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) > a=sqrt(diag(covmat))
> b=er
> k=1
> plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="white",lwd=1.5)
> polygon(u,v,border=NA,col=rgb(0,0,1,.3))
> for(i in 1:100){
+ id=sample(nrow(rzoo),replace=TRUE)
+ covmat=Moments(as.matrix(rzoo)[id,],"CovClassic")
+ er=apply(as.matrix(rzoo)[id,],2,mean)
+ points(sqrt(diag(covmat))[k],er[k],cex=.5)
+ }

or for another asset

Here is what we got on the (estimated) efficient frontier

> covmat=Moments(as.matrix(rzoo),"CovClassic")
> er=apply(as.matrix(rzoo),2,mean)
> ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50) > plot(ef$sd,ef$er,type="l",xlab="Standard Deviation",ylab="Expected Return", xlim=c(3.5,11),ylim=c(0,2.5),col="white",lwd=1.5)
> points(sqrt(diag(covmat)),er,pch=19,col="blue")
> text(sqrt(diag(covmat)),er,names(er),pos=4, col="blue",cex=.6)
> polygon(u,v,border=NA,col=rgb(0,0,1,.3))
> for(i in 1:100){
+ id=sample(nrow(rzoo),replace=TRUE)
+ covmat=Moments(as.matrix(rzoo)[id,],"CovClassic")
+ er=apply(as.matrix(rzoo)[id,],2,mean)
+ ef <- efficient.frontier(er, covmat, alpha.min=-2.5, alpha.max=2.5, nport=50)
+ lines(ef$sd,ef$er,col="red")
+ }

Thus, it is somehow rather difficult to assess wheter a portfolio is optimal, or not… At least from a statistical perspective….

Optimal Portfolios #1

This afternoon, I will start a crash course on financial portfolio optimization, with application in R. This week, we start with simple things, with the theoretical setup, without and with a risk free asset. We will discuss then the problem of estimating parameters, in a robust way. Then we introduce the idea of consider a more general criteria to quantify risk than the variance (but it means more general distributions… this point will be discussed further next time). The slides are available here, and R codes from there (in a Markdown)

Traffic Flow of Kota Kinabalu (with R)

This morning, we had our first practicals on network flows, using  an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map,

m=matrix(c(0,5.995910, 116.105520,
1,5.992737, 116.093718,
2,5.992066, 116.109883,
3,5.976947, 116.095760,
4,5.985766, 116.091580,
5,5.988940, 116.080112,
6,5.968318, 116.080764,
7,5.977454, 116.075460,
8,5.974226, 116.073604,
9,5.969651, 116.073753,
10,5.972341, 116.069270,
11,5.978818, 116.072880),3,12)

we can be visualized below

library(OpenStreetMap)
map = openmap(c(lat= 6.000, lon= 116.06),
c(lat= 5.960, lon= 116.12))
map=openproj(map)
plot(map)
points(t(m[3:2,]),col="black", pch=19, cex=3 )
text(t(m[3:2,]),c("s",1:10,"t"),col="white")

If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least).

To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper

library(devtools)
install_github("ropensci/tabulizer")
library(tabulizer)
location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'
out <- extract_tables(location)

with Windows, it seems to be necessary to download another package first

library(devtools)
install_github("ropensci/tabulizerjars")
install_github("ropensci/tabulizer")
library(tabulizer)
location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'
out <- extract_tables(location)

Now we can get out data frame with capacities

B1=as.data.frame(out[[2]])
B2=as.data.frame(out[[3]])
E=data.frame(from=B1[3:20,"V3"],
to=B1[3:20,"V4"])
E=E[-c(6,8),]
capacity=as.character(B2$V3[-1])
capacity[6]="843"
capacity[4]="2913"
E$capacity=as.numeric(capacity)

We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read)

plot(map)
points(t(m[3:2,]),col="black", pch=19, cex=3 )
B=data.frame(i=as.character(c("s",paste("V",1:10,sep=""),"t")),
x=m[3,],y=m[2,])
for(i in 1:nrow(E)){
i1=which(B$i==as.character(E$from[i]))
i2=which(B$i==as.character(E$to[i]))
segments(B[i1,"x"],B[i1,"y"],B[i2,"x"],B[i2,"y"],lwd=3)
}
text(t(m[3:2,]),c("s",1:10,"t"),col="white")

To get the graph with capacities, an alternative is to use

library(igraph)
g=graph_from_data_frame(E)
E(g)$label=E$capacity
plot(g)

but it does not respect geographical locations of nodes. It can actually be done using

plot(g, layout=as.matrix(B[,c("x","y")]))

To get a better understanding of the capacities of the road, use

plot(g, layout=as.matrix(B[,c("x","y")]),
edge.width=E$capacity/200)

From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using

> (m=max_flow(graph=g, source="s", target="t"))
$value
[1] 2571

$flow
[1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422
[16] 1149

Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers  max flow min cut theorem to… and application of the Shortest Path… (“the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables.

E$flux1=m$flow
E(g)$label=E$flux1
plot(g, layout=as.matrix(B[,c("x","y")]),
edge.width=E$flux1/200)

That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value

E$flux2=c(1422,1149,1422,1149,0,0,0,0,
1149,1422,1149,0,0,1149,1422,1149)
E(g)$label=E$flux2
plot(g, layout=as.matrix(B[,c("x","y")]),
edge.width=E$flux2/200)

Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu.

location <- 'http://www.worldresearchlibrary.org/up_proc/pdf/999-150486366625-30.pdf'
out <- extract_tables(location)
dim(out[[3]])
B1=as.data.frame(out[[3]])
E=data.frame(from=B1[2:61,"V2"],
to=B1[2:61,"V3"],
capacity=B1[2:61,"V4"])
E$capacity=as.numeric(
as.character(E$capacity))
library(igraph)
g=graph_from_data_frame(E)
m=max_flow(graph=g,
source="S",
target="T")
E$flux1=m$flow
E(g)$label=E$flux1
plot(g,
edge.width=E$flux1/200,
edge.arrow.size=0.15)

Here the value of the maximal flow is 4017, just as they found in the original paper