Freakonometrics

To content | To menu | To search

Thursday, June 21 2012

A la recherche des groupes de trolls sur Twitter

Il y a quelques jours avait mis en ligne sur son blog http://olihb.com/ une très jolie carte permettant de visualiser les personnes (ou les comptes) impliquées sur Twitter sur les manifestations récentes au Québec (la carte était basée sur gephi). Beaucoup de monde twittait avec des hashtags, et l'idée était d'utiliser les plus classiques et les plus populaires pour identifier les comptes actifs (#ggi, #manifencours, etc). Mais la semaine passée j'avais l'impression que le débat s’était durci, avec beaucoup plus de trolls dans les discussions. 

Bref, le principe est que les trolls évoluent souvent en petits groupes, donc avec @3wen on a voulu voir si on pouvait identifier facilement les groupes de trolls (sans pour autant faire du sentiment analysis, juste visualiser des personnes qui se suivent beaucoup entre elles, les connivences entre twittos, i.e. du suivi mutuel).

La méthodologie est simple. On a deux bases à notre disposition. Un première contient tous les tweets contenant #ggi entre le 7 juin à 1 heure du matin et le 8 à minuit 45 (soit 14,700 tweets). La seconde contient tous les tweets contenant #ggi entre le 12 juin à 23 heures 30 et le 14 juin à 17 heures (soit 14,784 tweets). Ah oui, les heures sont GMT. Plus précisément, on a la fréquence d'arrivée suivante

Les durées d'observations sont courtes, mais l'idée est que si une personne a tweetté au moins une fois pendant ces périodes, on les a dans notre base, et ça suffit pour les suivre. On ne va pas étudier l'intensité de la mise en ligne de tweets. Je ne rentre pas dans le détail du code, mais c'est basé sur les fonctions présentées dans un précédant billet, et toujours avec la même librairie,

require("RJSONIO")

Sur des deux bases, i.e. 29,484 tweets, on récupère un ensemble de 6,556 twittos uniques (mais répartis dans deux listes: ceux qui ont tweeté le 7 et ceux qui ont tweeté le 14). Pour toutes ces personnes, on est allé voir qui les suivaient. Sur le principe, c'est simple, mais on fonctionne avec des API qui imposent des limites horaires. Bref, il faut bricoler. Pour une personne, on récupère avec la fonction suivante la liste de ses followers,

recuperefollowers=function(id,nbrequetes){
followers=try(scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=-1&user_id=",
id,sep=""),what = "character", encoding="latin1"))
if(is.null(attr(followers,"condition"))){
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=followers$id
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
plusDe5000=FALSE
while(nbrequetes<150 & next_cursor!="0" &
is.null(attr(followers,"condition"))){
plusDe5000=TRUE
followers=scan(paste(
"http://api.twitter.com/1/followers/ids.json?cursor=",
next_cursor,"&user_id=",id,sep="")
,what = "character", encoding="latin1")
followers=paste(followers[1:length(followers)],collapse=" ")
followers=fromJSON(followers, method = "C")
id_followers=c(id_followers,followers$id)
next_cursor=followers$next_cursor_str
nbrequetes=nbrequetes+1
}
if(plusDe5000 & nbrequetes>=149){
return(c(list(),nbrequetes,TRUE))
}else{
return(c(list(id_followers),nbrequetes,FALSE))
}}else{
nbrequetes=nbrequetes+1
print(paste("Probleme rencontre avec ",id,sep=""))
return(c(list("PROBLEME"),nbrequetes,FALSE))}}

On peut ensuite lancer cette fonction sur tous les identifiants qu'on a récupéré, et on stocke tout le monde dans une (petite) base

temp=try(recuperefollowers(
lesIDAParcourir[unID_index],nbrequetes))
lesFollowers=temp[[1]]
resul=cbind(rep(lesIDAParcourir[unID_index]
,length(lesFollowers)),lesFollowers)
nbrequetes=temp[[2]]
refaire=temp[[3]]
Sys.sleep(runif(1,1,1.5))
compte=compte+1
write.table(resul,paste("auto_resul_",
compte,".txt",sep=""),sep=";")

Si on veut que ça tourne, il faut juste faire patienter un heure après avoir fait 150 requêtes. Bref, dans la boucle, on met une petite fonction qui temporise,

if(nbrequetes>148){
nbrequetes=1
Sys.sleep(60*60+60*trunc(runif(1,2,3))) }

Une fois créées les 6,500 bases de nos 6,500 comptes, contenant la liste de tous les followers, on va les agréger (fort heureusement, on n'avait pas les  gros comptes avec plus d'un million de followers). La liste se récupère avec 

N=list.files("id/petits")

D'ailleurs, pour recoller les 6,500 bases ensemble, Ewen suggérait le code suivant, beaucoup beaucoup plus rapide qu'une gros boucle,

recupere_petits_fichiers=function(x){
temp=read.table(paste("id/petits/",x,sep=""),
header=TRUE,sep=";",comment.char="",
check.names=FALSE,colClasses=
c("character","character"))
colnames(temp)=list("id_twittos","id_follower")
head(temp)

On fait ensuite une récupération des followers pour la première base de tweets, et une autre pour ceux de la seconde. Les deux listes sont obtenues avec le code suivant

present_avant_temp=temp[which(
temp$id_twittos%in%lesID_avant),]
if(nrow(present_avant_temp)>0){
present_avant_temp=
present_avant_temp[which(
present_avant_temp$id_follower%in%lesID_avant),]}
else{
present_avant_temp=NULL}
 
present_apres_temp=
temp[which(temp$id_twittos%in%lesID_apres),]
if(nrow(present_apres_temp)>0){present_apres_temp=
present_apres_temp[
which(present_apres_temp$id_follower%in%lesID_apres),]}
else{present_apres_temp=NULL}
return(c(list(present_avant_temp),
list(present_apres_temp),
list(nbfollowers=nrow(temp))))

Pour constituer notre grosse base, la ruse est d'utiliser non pas une boucle sur tous les petites bases, mais

resul=lapply(N,recupere_petits_fichiers)
liens_avant=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[1]])),
stringsAsFactors=FALSE)
liens_apres=data.frame(do.call(
"rbind",lapply(resul,function(x) x[[2]])),
stringsAsFactors=FALSE)

On a ainsi récupéré une liste avec plusieurs centaines de milliers de comptes Twitter. On a retenu seulement ceux qui avaient twitté avec le hashtag #ggi. On a ainsi pu constituer, pour nos 6,556 twittos une base de personnes qui tweettent, et de personnes qui les suivent (et qui ont twitté pendant cette - courte - période avec le hastag #ggi). Bref, une grosse matrice de correspondance de qui suit qui. On a pu le faire sur les deux bases. C'est précisément la différence entre les deux bases que l'on a cherché à visualiser, et gephi semble permettre de comparer les bases dans le temps. A gauche, on a le nuage de connexions entre les comptes pour la première période, et à droite pour la seconde. 

Bon, le premier est plus compact, mais on voit qu'il y a des différences, en particulier sur les bords, que ce soit en haut

ou à droite du nuage,

Dans les deux cas, on voit que plein de petits comptes sont venus se greffer au nuage, souvent en suivant quelqu'un qui le suit en retour. Et parfois, ce sont les seuls connexions qu'ils ont (qui utilisent le hashtag #ggi). Ce sont ces comptes qui pourraient constituer les nouveaux trolls qu'on espérait identifier. L'avantage est qu'on a les noms de tous ces comptes (on ne les a pas mis sur le graphique), mais dans un second temps, on va aller regarder ce qu'on fait précisément ces comptes. Peut-être interprétons nous mal ces graphiques, qui représentent une arrivée de nouveaux comptes, et que cette dynamique est la même pour tous les nouveaux comptes, qui commencent par de faibles connexions, et ensuite vont s'étoffer au fur et à mesure. Ce qu'il faut que l'on comprenne aussi c'est pourquoi ils sont tous aussi rapprochés. Ils semblent être suivis par des personnes qui sont suivis pas les mêmes personnes. Ce qui laisse à penser que l'on a effectivement identifié une certaine catégorie de comptes.
On peut reprendre les graphs avec une résolution plus fine, avant

et après

ou alors en changeant les options, avant

à comparer avec après

On voit peut être encore plus clairement sur ces deux dernières images les nouveaux comptes se greffer sur le graph existant. Si la qualité ne suffit pas, des graphiques (très) haute résolution sont téléchargeables (15Mo chacun), avant et après (avec les identifiants cette fois). Sinon, en promenant la souris sur le dessin dessous, on voit encore mieux la différence,

Maintenant - pour être complètement honnête - sur la lecture des graphs, j'avoue ne pas avoir tout bien compris. Les représentations semblent claires, mais j'aurais bien aimé être certain d'avoir compris la construction proposée par gephi avant de les mettre en ligne. Car si on regarde sur moins ce compte, on comprend que l'algorithme permet effectivement de bien mettre en avant les réseaux. Parmi les illustrations les plus frappantes, la figure suivante montre des interconnections avec Facebook à gauche, et purement aléatoire à droite (illustration trouvé sur le site de gephi),

La grande difficulté est alors de représenter spatialement les groupes, comme le notent Gastner et Newman dans The Spatial Structure of NetworksWatts et Strogatz dans Collective dynamics of 'small-world' networks ou encore Nisha et Venkatesh dans Small worlds: How and why. Mais pour etre tout à fait complet, je renverrais aussi vers la thèse de doctorat de Jure Leskovec, intitulée Dynamics of large networks ainsi qu'au livre (intégralement téléchargeable) Networks, Crowds, and Markets: Reasoning About a Highly Connected World de David Easley et Jon Kleinberg. Bref, on a découvert un outils probablement très riche, mais il va nous falloir du temps pour comprendre ce qu'il fait vraiment.

Monday, June 11 2012

Do you still have time to sleep ?

Last week, @3wen (Ewen) helped me to write nice R functions to extract tweets in R and build datasets containing a lot of information. I've tried a couple of time on my own. Once on tweet contents, but it was not convincing and once on the activity on Twitter following an event (e.g. the death of someone famous). I have to admit that I am not a big fan of databases that can be generated using standard function to study tweets. For instance, we can only extract tweets, not re-tweets (which is also an important indicator of tweet-activity). @3wen suggested to use
require("RJSONIO") 
The first step is to extract some information from a tweet, and store it in a dataset (details can be found on https://dev.twitter.com/)
obtenir_ligne <- function(unTweet){
date_courante=unTweet$created_at
id_courant=unTweet$id_str
text=unTweet$text
nb_followers=unTweet$user$followers_count
nb_amis=unTweet$user$friends_count
utc_offset=unTweet$user$utc_offset
listeMentions=unTweet$entities$user_mentions
return(c(list(c(id_courant,date_courante,text,
nb_followers,nb_amis,utc_offset)),
list(do.call("rbind",lapply(listeMentions,
function(x,id_courant) c(id_courant,
x$screen_name),unTweet$id_str)))))
}
Now that we  have the code to extract information from one tweet, let us find several tweets, from one user, say my account,
nom="Freakonometrics"
The (small) problem here, is that we have a limitation: we can only get 100 tweets per call of the function
n=100
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,sep=""),what = "character",
encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")
Then, we use our function to build a database with 100 lines,
extracTweets <- lapply(tweets_courants,
obtenir_ligne)
mentions=do.call("rbind",lapply(extracTweets,
function(x) x[[2]]))
colnames(mentions)=list("id","screen_name")
res=t(sapply(extracTweets,function(x) x[[1]]))
colnames(res) <- list("id","date","text",
"nb_followers","nb_amis","utc_offset")
The idea then is simply to use a loop, based on the latest id observed
dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str
So, here we go,
compteurLimite=100
 
while(compteurLimite<4100){
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,"&max_id=",dernier_id,sep=""),
what = "character", encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")
 
extracTweets <- lapply(tweets_courants[
2:length(tweets_courants)],obtenir_ligne)
mentions=rbind(mentions,do.call("rbind",
lapply(extracTweets,function(x) x[[2]])))
res=rbind(res,t(sapply(extracTweets,function(x) x[[1]])))
t(sapply(extracTweets,function(x) x[[1]]))
dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str
compteurLimite=compteurLimite+100
}
 
resFreakonometrics=res=
data.frame(res,stringsAsFactors=FALSE)
All the information about my own tweets (and re-tweets) are stored in a nice dataset. Actually, we have even more, since we have extracted also names of people mentioned in tweets,
mentionsFreakonometrics=
data.frame(mentions)
We can look at people I mention in my tweets
gazouillis=sapply(split(mentionsFreakonometrics,
mentions$screen_name),nrow)
gazouillis=gazouillis[order(gazouillis,
decreasing=TRUE)]
 
plot(gazouillis)
plot(gazouillis,log="xy")
> gazouillis[1:20]
tomroud freakonometrics       adelaigue       dmonniaux
155              84              77              56
J_P_Boucher         embruns      SkyZeLimit        coulmont
42              39              35              31
Fabrice_BM            3wen          obouba          msotod
31              30              29              27
StatFr     nholzschuch        renaudjf        squintar
26              25              23              23
Vicnent        pareto35        romainqc        valatini
23              22              22              22 
If we plot those frequencies, we can clearly observe a standard Pareto distribution,

Now, let us spend some time with dates and time of tweets (it was the initial goal of this post)... One more time, there is a (small) technical problem that we have to deal with: language. We need a function to convert date in English (on Twitter) to dates in French (since I have a French version of R),
changer_date_anglais <- function(date_courante){
mois <- c("Jan","Fév", "Mar", "Avr", "Mai",
"Jui", "Jul", "Aoû", "Sep", "Oct", "Nov", "Déc")
months <- c("Jan", "Feb", "Mar", "Apr", "May",
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
jours <- c("Lun","Mar","Mer","Jeu",
"Ven","Sam","Dim")
days <- c("Mon","Tue","Wed","Thu",
"Fri","Sat","Sun")
leJour <- substr(date_courante,1,3)
leMois <- substr(date_courante,5,7)
return(paste(jours[match(leJour,days)]," ",
mois[match(leMois,months)],substr(
date_courante,8,nchar(date_courante)),sep=""))
}
So now, it is possible to plot the times where I am online, tweeting,
DATE=Vectorize(changer_date_anglais)(res$date)
DATE=sapply(resSkyZeLimit$date,
changer_date_anglais,simplify=TRUE)
 
DATE2=strptime(as.character(DATE),
"%a %b %d %H:%M:%S %z %Y")
lt= as.POSIXlt(DATE2, origin="1970-01-01")
heure=lt$hour+lt$min/60
plot(DATE2,heure)

On this graph, we can see that I am clearly not online almost 6 hours a day (or at least not on Twitter). It is possible to visualize more precisely the period of the day where I might be on Twitter,

hist(heure,breaks=0:24,col="light green",proba=TRUE)
X=c(heure-24,heure,heure+24)
d=density(X,n = 512, from=0, to=24,bw=1)
lines(d$x,d$y*3,lwd=3,col="red")

or, if we want to illustrate with some kind of heat plot,

Note that we did it for my Twitter account, but we can also run the code on (almost) anyone on Twitter. Consider e.g. @adelaigue. Since Alexandre is tweeting in France, we have to play with time-zones,
res=extractR("adelaigue")
DATE=Vectorize(changer_date_anglais)(res$date) DATE2=strptime(as.character(DATE), "%a %b %d %H:%M:%S %z %Y",tz = "GMT")+2*60*60

or I can also look at @skythelimit who's usually twitting from Singapore (I am in Montréal). I can seen clearly when we might have overlaps,

res=extractR("skythelimit")

Nice isn't it. But it is possible to do much better... for instance, for those who do not ask specifically not to be Geo-located, we can see where they do tweet during the day, and during the night... I am quite sure a dozen posts with those functions can be written...

Thursday, May 24 2012

French dataset: population and GPS coordinates

A short post today based on recent work by @3wen (Ewen Gallic, graduate Student in Rennes, spending a year in Montreal). Since we were working on a detailed French dataset (per commune), we needed a dataset containing a list all communes, with population and location. GPS coordinates were extracted from Google, using the following php file, inspired by http://www.andrew-kirkpatrick.com/ on Google geocoding api with php webpage. Population was interpolated from INSEE's datasets, i.e. http://www.insee.fr/ (since data are over a 35 year period, from 1975 to 2010, changes have been taken into account as carefully are possible - e.g. merges and splits of cities - based on that description). A spline model has been used for all cities (with three degrees of freedom, and null and negative interpolation became one, since we'll be using loglinear models afterwards). Names are from that dataset, still on INSEE's website, http://www.insee.fr/.

A zipped file can be downloaded here popfr19752010.zip, but it is also possible to use the code below (it is a 24Mo dataset). Since it was hard to find such a dataset online (different files can be found, but we found none with population and location), we have decided to upload that dataset. Please let us know if there are problems with those data...

> base=read.csv(
+
"http://freakonometrics.free.fr/popfr19752010.csv", + header=TRUE)

Using that code, it is possible to locate all the communes in France (metropolitan), for instance

> library(maps)
> map("france")
> points(base$long,base$lat,cex=.1,col="red",pch=19)
> points(base$long,base$lat,cex=2*base$pop_2010/
+ max(base$pop_2010),col="blue",pch=19)

Several additional lines of code on that dataset (and also others) will be uploaded, soon.

Cette oeuvre est mise à disposition sous licence Paternité - Partage à l'Identique 3.0 non transposé. Pour voir une copie de cette licence, visitez http://creativecommons.org/. Date : 24 mai 2012, par Ewen GALLIC. Sources : INSEE, API Google Maps v3 et GeoHack (coordonnées GPS), propres calculs (estimation de population à partir des données INSEE).
  • reg : code region INSEE (character)
  • dep : code departement INSEE (character, corse 201 et 202 au lieu de 2A et 2B)
  • com : code commune INSEE (character)
  • article : article du nom de la commune (character)
  • com_nom : nom de la commune (character)
  • long : longitude (numeric)
  • lat : latitude (numeric)
  • pop_i : estimation de la population à la date i (ramenée à 1 si <=0), i=1975,...,2010 (numeric)

Wednesday, December 21 2011

Maps with R, and polygon boundaries

With R, it is extremely easy to draw maps. Let us start with something simple, like French regions. Baptiste mentioned on his blog that shapefiles can be downloaded from http://ign.fr/ website. Hence, if you extract the zip file, it is possible to get claims frequency per region (as done in the course ACT2040),

> library(maptools)
> library(maps)
> departements<-readShapeSpatial("DEPARTEMENT.SHP")
> region<-tapply(baseFREQ[,"nbre"],
+ as.factor(baseFREQ[,"region"]),sum)/
+ tapply(baseFREQ[,"exposition"],
+ as.factor(baseFREQ[,"region"]),sum)
> depFREQ=rep(NA,nrow(departements))
> names(depFREQ)=as.character(
+ departements$CODE_REG)
> for(nom in names(region)){
+ depFREQ[names(depFREQ)==nom] =
+ region[nom]}
> plot(departements,col=gray((depFREQ-.05)*20))
> legend(166963,6561753,legend=seq(1,0,by=-.1)/20+.05,
+ fill=gray(seq(1,0,by=-.1)),cex=1.25, bty="n")

Another application is on earthquakes. It is possible to use shapefiles of tectonic plates contour, and to relate earthquakes to plates. Shapefiles can be found on http://www.colorado.edu/ (here).

http://freakonometrics.blog.free.fr/public/perso4/plate-tekto.gif

First, we can extract the shapes of the tectonic plates

> plates = readShapePoly("plates.shp",
+ proj4string=CRS("+proj=longlat"))
> PP=SpatialPolygons2PolySet(plates)

Consider Montreal,

> montreal=c(-73.600,45.500)

Given that specific location, it is possible to use the following code to get the associated plate,

> PLATE.loc=function(pt){
+ K=NA
+ for(k in 1:17){
+ c=point.in.polygon(pt[1], pt[2],
+ PP[PP$PID==k,c("X")],PP[PP$PID==k,c("Y")],
+ mode.checked=FALSE)
+ if(c>0){K=k}
+ }
+ return(K)}
> abline(v=montreal[1],col="red")
> abline(h=montreal[2],col="red")
> PLATE.loc(montreal)
[1] 1

and then to plot the associated tectonic plate very easily

> PLATE=function(k0){
+ library(maps)
+ map("world")
+ polygon(PP[PP$PID==k0,c("X")],PP[PP$PID==k0,c("Y")],
+ col="red")
+ for(k in (1:17)[-k0]){polygon(PP[PP$PID==k,c("X")],
+ PP[PP$PID==k,c("Y")],col="light blue")}
+ map("world",add=TRUE)}
> PLATE(PLATE.loc(montreal))

Those code were used in the paper written with Mathieu, and that will be presented on January 30th at the Geotop seminar.

Tuesday, May 3 2011

Playing with robots

My son would be extremely proud if I tell him I can spend hours building robots. Well, my robots are not as fancy as Dr Tenma's, but they usually do what I ask them to do. For instance, it is extremely simple to build a robot with R, to extract data from websites. I have mentioned it here (one tennis matches), but it failed there (on NY Marathon). To illustrate the use of robots, assume that one wants to build his own dataset to study prices of airline tickets. First, we have to choose a departure city (e.g. Paris) and an arrival city (e.g. Montreal). Then, one wants to look at all possible dates from April first (I ran it last month) till the end of December (so we create a vector with all leaving dates, namely a vector for the day, one for the month, and one for the year). Then, we choose a return date (say 3 days after).
DEP="Paris"
ARR="Montreal"
DATE1D=rep(c(1:30,1:31,1:30,1:31,1:31,1:30,1:31,1:30,
1:31,1:31,1:29),3)
DATE1M=rep(c(rep(4,30),rep(5,31),rep(6,30),rep(7,31),
rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31),
rep(1,31),rep(2,29)),3)
DATE1Y=rep(c(rep(2011,30+31+30+31+31+30+31+
30+31+31+28),rep(2012,31+29)),3)
k=3
DATE3D=c((1+k):30,1:31,1:30,1:31,1:31,1:30,1:31,
1:30,1:31,1:31,1:29,1:k)
DATE3M=c(rep(4,30-k),rep(5,31),rep(6,30),rep(7,31),rep(8,31),
rep(9,30),rep(10,31),rep(11,30),rep(12,31),rep(1,31),rep(2,29),
rep(3,k))
DATE3Y=c(rep(2011,30+31+30+31+31+30+31+30+31+
31+28-k),rep(2012,31+29+k))

It is also possible (for a nice robot), to skip all prior dates

skip=max(as.numeric(Sys.Date()-as.Date("2011-04-01")),1)

Then, we need a website where requests can be written nicely (with cities and dates appearing explicitly). Here, I cannot not mention the website that I used since it is stated on the website that it is strictly forbidden to run automatic requests... Anyway, consider a loop create a url address (actually I chose the value of the date randomly, since I had been told that those websites had memory: if you ask too many times for the same thing during a short period of time, prices would go up),

URL=paste("http://www.♦♦♦♦/dest.dll?qscr=fx&flag=q&city1=",
DEP,"&citd1=",ARR,"&",
"date1=",DATE1D[s],"/",DATE1M[s],"/",DATE1Y[s],
"&date2=",DATE3D[s],"/",DATE3M[s],"/",DATE3Y[s],
"&cADULT=1",sep="")

then, we just have to scan the webpage, looking for ticket prices (just looking for some specific names)

page=as.character(scan(URL,what="character"))
I=which(page%in%c("Price0","Price1","Price2"))
if(length(I)>0){
PRIX=substr(page[I+1],2,nchar(page[I+1]))
if(PRIX[1]=="1"){PRIX=paste(PRIX,page[I+2],sep="")}
if(PRIX[1]=="2"){PRIX=paste(PRIX,page[I+2],sep="")}

Here, we have to be a bit cautious, if prices exceed 1000. Then, it is possible to start a statistical study. For instance, if we compare to destination (from Paris), e.g. Montréal and New York, we obtain the following patterns (with high prices during holidays),

It is also possible to run the code twice (here it was run last month, and a couple of days ago), for the same destination (from Paris to Montréal),


Of course, it would be great if I could run that code say every week, to build up a nice dataset, and to study the dynamic of prices...

The problem is that it is forbidden to do this. In fact, on the website, it is mentioned that if we want to extract data (for an academic purpose), it is possible to ask for an extraction. But if we do tell that we study specific prices, data might be biased. So the good idea would be to use several servers, to make several requests, randomly, and to collect them (changing dates and destination). But here, my computing skills - unfortunately - reach a limit....

Saturday, January 15 2011

Itération de calculs de quantiles

Un billet rapide pour répondre à un commentaire de Claire (ici) sur les quantiles, et les régressions quantiles. En probabilité, il est possible d'itérer pour obtenir un quantile. Par exemple, le quantile à 75% est tel que

http://freakonometrics.blog.free.fr/public/perso/quantit1.png
(en supposant la loi continue, afin de simplifier l'exposé). Mais une autre façon de le définir est de prendre la médiane des observations qui dépassent la médiane, i.e. rechercher http://freakonometrics.blog.free.fr/public/perso/quantit3.png tel que
http://freakonometrics.blog.free.fr/public/perso/quantit2.png
Notons que sur une estimation de quantile(s) on peut éventuellement avoir une petite différence.
> s1=.5
> s2=1-(1-.75)/(1-s1)
> X=rexp(100000)
> quantile(X[X>quantile(X,s1)],s2)
50%
1.392084
 
> quantile(X,.75)
75%
1.392062
Regardons maintenant le cas de la régression quantile. Par exemple une régression quantile de niveau 90%
library(quantreg)
plot(cars)
abline(rq(dist~speed,tau=.9,data=cars),
col="red",lwd=2)

On peut tenter de faire une régression quantile à 80% sur les observations dépassant le quantile à 50%
plot(cars)
abline(rq(dist~speed,tau=.9,data=cars),
col="red",lwd=.7)
abline(rq(dist~speed,tau=.5,data=cars),
col="blue",lwd=2)
Xp=predict(rq(dist~speed,tau=.5,data=cars))
I=(cars$dist>Xp)
points(cars[I,],col="blue",pch=19)
abline(rq(dist~speed,tau=.8,data=cars[I,]),
col="red",lwd=2)

ou une régression quantile à 50% sur les observations dépassant le quantile à 80%
plot(cars)
abline(rq(dist~speed,tau=.9,data=cars),
col="red",lwd=.7)
abline(rq(dist~speed,tau=.8,data=cars),
col="blue",lwd=2)
Xp=predict(rq(dist~speed,tau=.8,data=cars))
I=(cars$dist>Xp)
points(cars[I,],col="blue",pch=19)
abline(rq(dist~speed,tau=.5,data=cars[I,]),
col="red",lwd=2)

Manifestement itérer n'est pas neutre (mais comme derrière se cachent des problèmes d'optimisation, on ne devrait pas être surpris). On peut aussi en tenter plein, pour voir tout ce qu'on pourrait obtenir,

plot(cars)
Y=cars$dist
X=cars$speed
for(s1 in seq(.01,.89,by=.01)){
s2=1-(1-.9)/(1-s1)
reg1=rq(Y~X,tau=s1)
I1=(Y>=predict(reg1))
Y1=Y[I1];X1=X[I1]
reg2=rq(Y1~X1,tau=s2)
abline(reg2,col="red")
}

Il est probablement possible d'utiliser le fait que l'estimateur change pour robustifier la régression quantile... mais c'est un autre sujet (ou un sujet trop long pour un billet du samedi soir)....

Tuesday, January 11 2011

Maps with R, part... n+1

Following the idea posted on James Cheshire's blog (here), I have tried to play a little bit with R and Google. And it works ! Consider for instance life expectancy at birth (that can be found - and downloaded - here). Using the following code, it is possible to generate an R graph in a web page,

install.packages("googleVis")
library(googleVis)
input=read.csv("UNdata_Export_20110111_101634359.txt",
sep=";",header=TRUE)
Map<- data.frame(input$Country.or.Area, input$Value)
names(Map)<- c("Country", "Age")
Geo=gvisGeoMap(Map, locationvar="Country", numvar="Age",
options=list(height=350, dataMode='regions'))
plot(Geo)
XX(ok, there are still some problems with some countries - e.g. Russia and the US - but it should not be that difficult to fix). Thanks James, this is awesome !

Wednesday, November 10 2010

Generating a quasi Poisson distribution, version 2

Here and there, I mentioned two codes to generated quasiPoisson random variables. And in both of them, the negative binomial approximation seems to be wrong. Recall that the negative binomial distribution is

http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg01.png
where
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg02.png
and in R, a negative binomial distribution can be parametrized using two parameters, out of the following ones
  • the sizehttp://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg03.png
  • the probabilityhttp://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg04.png
  • the mean
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg05.png
Recall also that the variance for a negative binomial distribution should be
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg06.png
Here, we consider a distribution such that http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg07.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg08.png. In the previous posts, I used
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg09.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg10.png
i.e.
rqpois = function(n, lambda, phi) {
mu = lambda
k = mu/(phi * mu - 1)
r1 = rnbinom(n, mu = mu, size = k)
r2 = rnbinom(n, size=phi*mu/(phi-1),prob=1/phi)
k = mu/phi/(1-1/phi)
r3 = rnbinom(n, mu = mu, size = k)
r4 = rnbinom(n, size=mu/phi/(1-1/phi),prob=1/phi)
r = cbind(r1,r2,r3,r4)
return(r)
}
but as we can see below, none of those two functions work,
> N=rqpois(1000000,2,4)
> mean(N[,1])
[1] 2.001992
> mean(N[,2])
[1] 8.000033
> var(N[,1])/ mean(N[,1])
[1] 7.97444
> var(N[,2])/ mean(N[,2])
[1] 4.002022
with the first one, the expected value is correct, while it is the overdispersion parameter for the second. Now, if we do the maths correctly, it comes
http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg11.png

http://perso.univ-rennes1.fr/arthur.charpentier/latex/bineg12.png
which should now work,
> mean(N[,3])
[1] 2.001667
> mean(N[,4])
[1] 2.002776
> var(N[,3])/ mean(N[,3])
[1] 3.999318
> var(N[,4])/ mean(N[,4])
[1] 4.009647
So, finally it is better when we do the maths well.

Friday, November 5 2010

Pretty R code in the blog

David Smith (alias @revodavid, see also on the Revolutions blog, here) pointed out that my R code was not easy to read (not only due to my computing skills, but mainly because of the typography I use). He suggested that I use the Pretty R tool (here). And I will...

So, just to answer quickly to a question I received by email (a few weeks ago, sorry for the delay), here is the code to get the following nice plot

library(evd); data(lossalae) 
x <- lossalae$Loss; y <- lossalae$ALAE
xhist <- hist(log(x), plot=FALSE)
yhist <- hist(log(y), plot=FALSE)
top <- max(c(xhist$counts, yhist$counts))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE),
c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
plot(x, y, xlab="", ylab="",log="xy",col="red")
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes=FALSE, ylim=c(0, top),
space=0,col="light green")
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes=FALSE, xlim=c(0, top),
space=0, horiz=TRUE,col="light blue")
I guess it is more convenient to read it, isn't it...

Wednesday, September 29 2010

Detecting distributions with infinite mean

In a post I published a few month ago (in French, here, based on some old paper, there), I mentioned a statistical procedure to test if the underlying distribution http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO02.png of an i.i.d. sample http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO01.png had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO03.png=200 observations generated from a Pareto distribution
http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO04.png
(here http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, as a start)
> a=2
> X=runif(200)^(-1/a)
Here, we will use the package developped by Mathieu Ribatet,
> library(RFA)
Le chargement a nécessité le package : tcltk
Chargement de Tcl/Tk... terminé
Message d'avis :
le package 'RFA' a été compilé avec la version R 2.10.1
  • Using Generalized Pareto Distribution (and LR test)
A first idea is to fit a GPD distribution on observations that exceed a threshold http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO11.png >1.
Since we would like to test http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png (against the assumption that the expected value is finite, i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO07.png), it is natural to consider the likelihood ratio
http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO08.png
Under the null hypothesis, the distribution of http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO09.png should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let  us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill's estimator) to consider http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO10.png, and thus, to fit a GPD on the http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png largest values. And then to plot all that on a graph (like the Hill plot)
> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }
Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model (http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png the likelihood ratio statistics has a chi-square distribution).
If http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png largest values (and a confidence interval for the estimator),
If http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance), we have
i.e. for those two models, we clearly reject the assumption of infinite mean (even if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png gets too close from 1, we should consider thresholds large enough). On the other hand, if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),
  • Using Hill's estimator
An alternative could be to use Hill's estimator (with Alexander McNeil's package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill's estimator,
> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)
Again, if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we clearly rejct the assumption of infinite mean,
and similarly, if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance)
Here the test is more robust than the one based on the GPD. And if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean), again we accept http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png,
Note that if http://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=0.7, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,
(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

Friday, September 24 2010

EM and mixture estimation, with R (part 2)

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM01.png
In order to describe how EM algorithm works, assume first that both http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png are perfectly known, and the mixture parameter is the only one we care about.
  • The simple model, with only one parameter that is unknown
Here, the likelihood is
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM04.png
so that we write the log likelihood as
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM05.png
which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (that cannot be observed), taking value when http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM07.png is drawn from http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and 0 if it is drawn from http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png. More generally (especially in the case we want to extend our model to 3, 4, ... mixtures), http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png.
With that notation, the likelihood becomes
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM10.png
and the log likelihood
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM11.png
the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png is known, denoted http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png. Then I can predict the value of http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png (i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM08.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM09.png) for all observations,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM14.png
So I can inject those values into my log likelihood, i.e. in
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM15.png
having maximum (no need to run numerical tools here)
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM16.png
that will be denoted http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM17.png. And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png is the best predictor of http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM19.png given my observations (as well as my belief in http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png). Then comes a maximization (M) step, where using http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM06.png, I can estimate probability http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM13.png.
  • A more general framework, all parameters are now unkown
So far, it was simple, since we assumed that http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png and  http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM30.png. Recall that we had http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png.
If http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png, instead of being in the segment http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM31.png was in http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM32.png, then we could have considered mean and standard deviations of observations such that http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=0, and similarly on the subset of observations such that http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png=1.
But we can't. So what can be done is to consider http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png as the weight we should give to observation i when estimating parameters of http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM02.png, and similarly, 1-http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM18.png would be weights given to observation i when estimating parameters of http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM03.png.
So we set, as before
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM33.png
and then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM34.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM35.png
and for the variance, well, it is a weighted mean again,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM36.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM37.png

and this is it.
  • Let us run the code on the same data as before
Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png and then http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }

Then I get http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM12.png, or http://perso.univ-rennes1.fr/arthur.charpentier/latex/mixEM99.png. So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we're doing much better ! Maybe we should also go for a Gibbs sampling procedure ?... next time, maybe....

Thursday, September 23 2010

Optimization and mixture estimation, with R

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised...

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.
The density is here http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-01.png proportional to
http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-02.png
The true model is http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-03.png, and http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
sd2 <- param[5]
+  -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 - p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+  }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])

Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+   fn = logvrainegX)

Actually, since I might have identification problems, I take either http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png or http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-06.png, depending whether http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-07.png or http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-08.png is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-05.png), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: the black vertical segments. Obviously, when the sample is not enough heterogeneous (i.e. http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-09.png and http://perso.univ-rennes1.fr/arthur.charpentier/latex/mix-ml-04.png rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical  starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again... so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad.... but the probability to be far away from the tru value is not small at all... except when the difference between the two means exceeds 3...
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad... and the two algorithm perform similarly (maybe the first one being a bit closer to the true parameter).

Thursday, September 9 2010

Le million ! le milllion !

Hier soir (ou ce matin, je suis perdu avec ce décalage horaire) Christelle me demandait de parler un peu de prévision avec R (ici). Au lieu de renvoyer vers l'aide en ligne, penons un exemple pratique (et simple, voire si possible intéressant): la fréquentation d'un blog (en l'occurrence http://blogperso.univ-rennes1.fr/arthur.charpentier/). Considérons le nombre de pages vues, par jour, selon Google Analytics.
> base=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/million.csv",
+ sep=";",header=TRUE
)> X=base$nombre
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X)
> plot(D,X)
> plot(acf(X,lag=90))
La série a l'allure suivante (oui, le compteur n'a été installé qu'il y a deux ans),

Les autocorrélations montre une forte saisonnalité hebdomadaire, avec moins de consultations en semaine que le week-end,

La question que l'on cherche à résoudre est "ai-je une chance d'atteindre le million de pages vues d'ici la fin de l'année ?".

On peut traduire cette question de deux manières
  • quelle est la probabilité que le 1er janvier, j'ai atteint le million de pages vues,
  • quelle est la probabilité que la date où le million de pages vues sera atteint soit avant le 1er janvier
Une fois formalisée la question, reste à faire un peu d'économétrie.
  • modélisation économétrique
En faisant simple et rapide, afin de prendre en compte la corrélation forte avec la semaine précédente, et le fait que l'on s'intéresse à la somme cumulée, on peut considérer un modèle ARIMA
  • avec un retard d'ordre 7 pour les composantes moyennes mobiles et autorégressives,
  • avec une série intégrée à l'ordre 1,
L'ajustement se fait de la manière suivante,
> X=cumsum(base$nombre)
> model  <- arima(X,c(7 ,   # partie AR
+                                             1,    # partie I
+                                             7))   # partie MA
et la prévision se fait via
> forecast <- predict(model,200)

Ensuite, ce n'est qu'une représentation graphique,
> u=max(D)+1:200
> polygon(  c(u,rev(u)),  c(forecast$pred - 1.96*forecast$se,
+  rev(forecast$pred + 1.96*forecast$se)), col = "yellow", border=NA)
> lines(u,forecast$pred,col="blue",lwd=2)
> lines(u,forecast$pred- 1.96*forecast$se,col="blue",lty=2)
> lines(u,forecast$pred+ 1.96*forecast$se,col="blue",lty=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"),col="red")
> abline(h=1000000,col="red")


Afin de répondre à la question posée, on peut étudier les différentes probabilités envisagées,
  • probabilité que le 1er janvier, j'ai atteint le million de pages vues
Dans un premier temps, on utilise la normalité des prédictions (en supposant une normalité du bruit) pour obtenir la loi de la prédiction à une date quelconque,
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> x=seq(800000,1100000,by=100)
> y=dnorm(x,mean=forecast$pred[k],sd=forecast$se[k])
> plot(x,y,type="l",lwd=2)
> x0=x[x>=1000000]
> polygon(  c(x0,rev(x0)), c(y[x>=1000000],rep(0,length(x0))), col = "yellow",border=NA)
> lines(x,y,type="l",lwd=2)
> abline(v=1000000)

Le cas de la probabilité d'atteindre plus d'un million de visiteur le 1er janvier est alors
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])
[1] 0.2604821

  • problème dual, 
Dans un second temps, on peut envisager une autre approche, consistant à se demander quelle pourrait être la loi de la date du jour où le million de pages vues sera atteint,
> P=rep(NA,300)
> for(k in 1:300){
+ P[k]=1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])}
> plot(max(D)+1:300,P,type="l",lwd=2)
> x=max(D)+1:300
> x0=x[x<=as.Date("01/01/2011","%d/%m/%Y")]
> polygon(  c(x0,rev(x0)), c(P[x<=as.Date("01/01/2011","%d/%m/%Y")],rep(0,length(x0))),  col = "yellow", border=NA)
> lines(max(D)+1:300,P,type="l",lwd=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"))

La probabilité que cette date soit antérieure au 1er janvier est alors
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.2604821
(ce qui, fort heureusement, correspond à la probabilité calculée par le problème primal). Bref, j'ai 1 chance sur 4 d'atteindre le million de pages vues avant la nouvelle année....
  • sensibilité à un changement de modèle
C'est bien joli tout ça, mais ces calculs sont largement soumis à la contrainte du choix de modèle que j'ai fait arbitrairement au début... on peut se demander si en changeant de modèle, les résultats changent sensiblement, ou pas. Au lieu de tenter un autre modèle ARIMA (voire SARIMA), j'ai préféré changer la série de référence... et me focaliser sur 2010 uniquement. D'un côté j'enlève les premières semaines où le niveau de fréquentation était très faible, de l'autre, je donne un poids très important aux vacances d'été, i.e. la période juin-août, pendant laquelle les internautes semblent moins sensibles à la modélisation économétrique.
Si l'on modélise la fréquentation pour 2010, seulement, on obtient

avec les distributions suivantes, tout d'abord pour la densité du nombre de visiteur atteint au 1er janvier 2011,

et pour la fonction de répartition de la date où sera atteint le million de pages vues,
soit une probabilité de l'ordre de 35%, dans les deux cas,
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.3531648
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast2$pred[k],sd=forecast2$se[k])
[1] 0.3531648


On obtient des probabilités relativement proches avec les deux modèles, et j'aurais envie de croire que l'objectif est envisageable.... "challenge accepted" comme dirait mon ami Barney Stinson (oui, c'est mon ami). Reste juste à trouver des sujets qui attireront du monde en cette fin d'année...

Tuesday, August 3 2010

Importer une base SAS sous R

Depuis pas mal de temps, je reçois des bases envoyées par différentes personnes, pour faire des stats, et si la plupart ont la bonne idée de me les envoyer directement en format csv, malheureusement, je reçois parfois des bases sas dont je ne sais trop que faire (car je n'ai pas sas).
La library(foreign) de R propose d'importer des bases SAS qui sont au format xport (ici). Mais généralement, les données SAS sont au format sas7bdat. Pendant l'année, j'avais l'habitude d'aller squatter la salle informatique qui dispose de SAS, pour ouvrir SAS et exporter les données en csv. Mais avec la fermeture estivale de la faculté, j'ai été un peu bloqué. J'ai alors découvert la version online de SAS sur l'ent de l'université. Malheureusement, il faut que les données soient dans un répertoire sur le réseau, ce qui nécessite au préalable d'envoyer des données sur le réseau au lieu de pointer sur un répertoire local (ce qui peut être loin, voir impossible si la base est trop grosse). Fort heureusement, l'autre jour, j'ai découvert,

On peut en effet télécharger gratuitement sur le site de www.sas.com un lecteur de base, qui permet ensuite d'exporter la base en csv... et ensuite de l'importer facilement sous R... Trop facile !

Thursday, March 18 2010

Les problèmes de mémoire avec R

J'avais plusieurs fois hésité à évoquer ce point, car mes connaissance sur le sujet sont assez limitées. Mais comme plusieurs personnes ont soulevé la question dans des mails, je vais faire un billet spécial. Commençons par la mémoire vive. Pour connaître la capacité actuelle, en megabytes, Mbs.
> memory.limit()
[1] 1015
que l'on peut ajuster
> memory.limit(size=2000)
[1] 2000
correspondant à 2Go de mémoire. On peut essayer d'aller plus loin, mais R semble alors irrité par notre méconnaissance de la gestion de la mémoire,
> memory.limit(size=200000)
Erreur dans memory.size(size) :
  ne soyez pas stupide ! Votre machine a une limite de mémoire adressable de 4 Go

Pour regarder un peu plus en détails ce qui se trouve dans la mémoire,
> memory.profile()
       NULL      symbol    pairlist     closure environment     promise    language     special
          1        5852       87097        2224         340        4295       24058          46
    builtin        char     logical     integer      double     complex   character         ...
        687      117878        3207        5313        2485           1       14449           0
        any        list  expression    bytecode externalptr     weakref         raw          S4
          0        1421          41           0         515         122         123         642
En particulier, on peut regarder dans la garbage collection,
> gc()
          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  270863  7.3     531268 14.2   407500 10.9
Vcells 6276691 47.9    7494157 57.2  6379649 48.7
Beaucoup d'objets peuvent coexister dans la mémoire,
> objects()
  [1] "a"             "A"             "a2"            "age"           "agebar"      
  [6] "al"            "an"            "au"            "b"             "B2"        

mais on peut tous les effacer, éventuellement
> rm(list=ls())
> objects()
character(0)
On notera que la mémoire a effectivement diminé,
> memory.profile()
       NULL      symbol    pairlist     closure environment     promise    language     special
          1        5855       71383        1753         296        4251       17600          46
    builtin        char     logical     integer      double     complex   character         ...
        687        6705        3037        4627         563           1       12463           0
        any        list  expression    bytecode externalptr     weakref         raw          S4
          0        1060           1           0         514         122         123         642
Si j'importe une grosse base de données, on voit qu'elle prend effectivement de la place,.
> SINISTRES = read.table("D:\ -data\\sinistres.txt",sep=";",header=FALSE)
> CONTRATS  = read.table("D:\ -data\\contrats.txt",sep=";",header=FALSE)
> head(SINISTRES)
  no nocontrat garantie    cout
1  1      3710      2DO 1291.76
2  2      3996      3VI   87.36
3  3      1552      1RC  995.20
4  4   1010996      2DO   77.33
5  5   4016701      1RC    0.00
6  6   4075803      2DO  512.74
On peut regarder la mémoire prise par ces deux objets,
> object.size(CONTRATS)
62382688 bytes
> object.size(SINISTRES)
2109728 bytes

Et si on fusionne des bases, on obtient un objet assez gros,
> object.size(BASECOUT)
11392528 bytes
qui est "proportionel" à la taille de la base
> nrow(BASECOUT)
[1] 105435
> nrow(CONTRATS)
[1] 678019
L'espace mémoire s'est considérablement rempli,
> gc()
          used (Mb) gc trigger  (Mb) max used  (Mb)
Ncells  140368  3.8     647796  17.3   813173  21.8
Vcells 9666689 73.8   25870145 197.4 25862489 197.4

Dans les bases, l'âge du conducteur est ici un entier,
> storage.mode(CONTRATS$agecond)
[1] "integer"

mais si je le change en réel, on peut noter que la base occupe davantage de place en mémoire,
> storage.mode(CONTRATS$agecond)="double"
> gc()
           used (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   140569  3.8     647796  17.3   813173  21.8
Vcells 10005749 76.4   28685834 218.9 26617390 203.1
> object.size(CONTRATS)
65094760 bytes
ou pour creuser un peu plus loin,
> gc(verbose=TRUE)
Garbage collection 67 = 11+3+53 (level 2) ...
3.8 Mbytes of cons cells used (22%)
76.4 Mbytes of vectors used (35%)
           used (Mb) gc trigger  (Mb) max used  (Mb)
Ncells   140588  3.8     647796  17.3   813173  21.8
Vcells 10005733 76.4   28685834 218.9 26617390 203.1
Bref, à forcer de jouer, je remplis la mémoire de R,
> memory.size()
[1] 86.15
>
memory.limit()
[1] 2000
J'ai un peu de mal à retrouver ces valeurs par ailleurs...
Mais sinon, il peut exister un autre problème avec R, problème que connaissent les chercheurs (et les élèves) de Rennes 1. Comme le disait Benoît, "tu devrais mettre cela sur ton blog : "Les chercheurs de la faculté des sciences économiques de l'université de rennes 1 n'ont qu'1 GO pour sauvegarder leur données: ils arrêtent leur recherche au bout d'une semaine faute de mémoire"". Mine de rien, c'est un problème vicieux, qui n'est pas un problème de mémoire vive, mais de mémoire physique sur l'espace de travail.
En fait, c'est un problème assez stupide, mais comme tous les problèmes stupides, ils font rapidement perdre une heure, voire davantage. Le but était d'exporter une base de données à un coauteur. Et le code R était le suivant
> dim(C)
[1] 3660  113
> write.table(C,file="MyDataBase.csv", sep=";")
> Mydata2=read.table("MyDataBase.csv",head=TRUE, sep=" ")
> dim(Mydata2)
[1] 2782    1

Autrement dit la base initiale contenait 3660 lignes, et 113 colonnes. On exporte la base, puis on la rimporte, et là, on n'a que 2782 lignes, un colonne.... et aucun message d'erreur (ou même de warning). En fait, le problème venait d'un disque plein.... on a mis du temps à comprendre, mais je me suis souvenu avoir le même problème lors d'un TD, où j'avais sauvé mon code sur mon espace réseau, sauf que l'espace était nul: R ne s'est jamais plein de ne pas pouvoir sauver le code faute de place. Je me suis retrouvé avec un fichier de taille 0 ko, entièrement vide... Trois heures de code perdues....

- page 1 of 2