Freakonometrics

To content | To menu | To search

Sunday, November 18 2012

In less than 50 days: 2013, year of statistics

A couple of days ago, Jeffrey sent me an e-mail, encouraging me to write about the international year of statistics, hoping that I would participate to the event. The goals of Statistics2013 are (as far as I understood) to increase public awareness of the power and impact of statistics on all aspects of society, nuture Statistics as a profession and promote creativity and development in Statistics.

So I will try to write more posts (yes, I surely can try) on statistical methodology, with sexy applications... So see you in a few days on this blog, to celebrate statistics ! I can even invite guests to write here... all contributors are welcome !


 

Wednesday, November 14 2012

I think I've already heard that tune...

Maybe I just wanted to talk about Vince Guaraldi in Movember (let's try this challenge for this month: mention only people wearing a moustache). This week-end, the CD of Charlie Brown and Chrismas (see e.g. on youtube) was playing at home (I know, it is a bit early). And during the second track, I said to myself "I know that tune". For some reasons, 5 or 6 notes reminded me of a song of Jacques Brel... What are the odds I thought (at first) ?

The difficult point is that calculating the probability that a given sequence appears over some runs might be compiicated (I am not that good when it is time to compute probabilities). But it is still possible to run some codes to estimate the probability to get a sequence over a fixed number of runs. For instance, consider the case there where there are only two notes (call them "H" and "T", just to illustrate). 

> runappear=function(seqA,nrun){
+ win=FALSE
+ seq=sample(c("H","T"),size=
+ length(seqA)-1,replace=TRUE)
+ i=length(seqA)-1
+ while((win==FALSE)&(i<nrun)){
+ i=i+1
+ seq=c(seq,sample(c("H","T"),size=1))
+ n=length(seq)
+ if(all(seq[(n-length(seqA)+1):n]==seqA)){
+ win=TRUE}}
+return(win)}

For instance, we can generate 100,000 samples to see how many times we have the tune (call it a tune) "THT" over 10 consecutive notes

> runappear(c("T","H","T"),10)
[1] TRUE
> simruns=function(seqA,nrun,sim){
+ s=0
+ for(i in 1:sim){s=s+runappear(seqA,nrun)};
+ return(s/sim)}
> simruns(c("T","H","T"),10,100000)
[1] 0.65684

But the most interesting point (from my point of view...) about those probabilities is that it yields a nontransitive game. Consider two players. Those two players - call them (A) and (B) - select tunes (it has to be three notes). If (A)'s tune appears before (B)'s, them (A) wins. The funny story is that there is no "optimal" strategy. More specificaly, if (B) knows what (A) has chosen, then it is always possible for (B) to choose a tune that will appear before (A)'s with more than 50% chance. Odd isn't it. This is Penney's game, as introduced by Walter Penney, see e.g. Nishiyama and Humble (2010).

It is possible to write a code where we run scenarios until one player wins,

> winner=function(seqA,seqB){
+ win="N"
+ seq=sample(c("H","T"),size=2,replace=TRUE)
+ while(win=="N"){
+ seq=c(seq,sample(c("H","T"),size=1))
+ n=length(seq)
+ if(all(seq[(n-2):n]==seqA)){win="A"}
+ if(all(seq[(n-2):n]==seqB)){win="B"}};
+ return(win)}
> A=c("H","H","H");B=c("H","T","H")
> winner(A,B)
[1] "A"

If we run one thousand games, we see here that (B) wins more frequently (with those tunes)

> game=function(seqA,seqB,n=1000){
+ win=rep(NA,n)
+ for(i in 1:n){win[i]=winner(seqA,seqB)}
+ return(table(win)/n)
+ }
> game(A,B)
win
A     B
0.403 0.597 

Let us look now at all possible tunes (8 if we consider a sequence of 3 notes), for players (A) and (B),

> A1=c("T","T","T")
> A2=c("T","T","H")
> A3=c("T","H","T")
> A4=c("H","T","T")
> A5=c("T","H","H")
> A6=c("H","T","H")
> A7=c("H","H","T")
> A8=c("H","H","H")
> B=data.frame(A1,A2,A3,A4,A5,A6,A7,A8)
> PROB=matrix(NA,8,8)
> for(i in 1:8){
+ for(j in 1:8){
+ PROB[i,j]=game(B[,i],B[,j],1000)[1]
+ }}

We have the following probabilities (that (A) wons over (B))

> PROB
[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,] 1.000 0.529 0.414 0.137 0.411 0.437 0.317 0.496
[2,] 0.493 1.000 0.674 0.242 0.664 0.625 0.480 0.701
[3,] 0.615 0.331 1.000 0.516 0.492 0.503 0.365 0.583
[4,] 0.885 0.735 0.518 1.000 0.494 0.513 0.339 0.580
[5,] 0.603 0.342 0.487 0.500 1.000 0.465 0.762 0.862
[6,] 0.602 0.360 0.527 0.506 0.507 1.000 0.329 0.567
[7,] 0.695 0.511 0.642 0.664 0.288 0.635 1.000 0.510
[8,] 0.502 0.279 0.428 0.382 0.117 0.406 0.487 1.000

It is possible to visualize the probabilities on the graph below

with below, the optimal strategy, that (B) should choose, given the tune chosen by (A),

(actually it is possible to compute explicitely those probabilites). Here, it is possible for (B) to choose a tune so that the probability that (A) wins is less than 50% (even less than 35%).

But we are a bit far away from the problem mentioned earlier... Furthermore, of course, assuming independence and randomness (in music) is extremely simplistic (see a short discussion here, a few weeks ago, in French). We might expected some Markov chain behind, in order to have something which will not hurt our ears, as stated in Brooks et al. (1957) for instance (quite a long time ago actually).

To conclude, about Jacques Brel's song: actually, Amsterdam is inspired (somehow) from Greensleeves, the traditional English folk song. And Vince Guaraldi played it on that CD...

So yes, there is a connection between all those songs. And no need to do maths to get it...

Saturday, November 3 2012

Normality versus goodness-of-fit tests

In many cases, in statistical modeling, we would like to test whether the underlying distribution from an i.i.d. sample lies in a given (parametric) family, e.g. the Gaussian family

where

Consider a sample

> library(nortest)
> X=rnorm(n)

Then a natural idea is to use goodness of fit tests (natural is not necessarily correct, we'll get back on that later on), i.e.

for some  and . But since those two parameters are unknown, it is not uncommon to see people substituting estimators to those two unknown parameters, i.e.

Using Kolmogorov-Smirnov test, we get

> pn=function(x){pnorm(x,mean(X),sd(X))};
> P.KS.Norm.estimated.param=
+ ks.test(X,pn)$p.value

But since we choose parameters based on the sample we use to run a goodness of fit test, we should expect to have troubles, somewhere. So another natural idea is to split the sample: the first half will be used to estimate the parameters, and then, we use the second half to run a goodness of fit test (e.g. using Kolmogorov-Smirnov test)

> pn=function(x){pnorm(x,mean(X[1:(n/2)]),
+ sd(X[1:(n/2)]))}
> P.KS.Norm.out.of.sample=
+ ks.test(X[(n/2+1):n],pn)$p.value>.05)

As a benchmark, we can use Lilliefors test, where the distribution of Kolmogorov-Smirnov statistics is corrected to take into account the fact that we use estimators of parameters

> P.Lilliefors.Norm=
+ lillie.test(X)$p.value 

Here, let us consider i.i.d. samples of size 200 (100,000 samples were generated here). The distribution of the -value of the test is shown below,

In red, the Lilliefors test, where we see that the correction works well: the -value is uniformly distributed on the unit inteval. There is 95% chance to accept the normality assumption if we accept it when the -value exceeds 5%. On the other hand,

  • with Kolmogorv-Smirnov test, on the overall sample, we always accept the normal assumption (almost), with a lot of extremely large -values
  • with Kolmogorv-Smirnov test, with out of sample estimation, we actually observe the opposite: in a lot of simulation, the -value is lower then 5% (with the sample was from a  sample).

The cumulative distribution function of the -value is

I.e., the proportion of samples with -value exceeding 5% is 95% for Lilliefors test (as expected), while it is 85% for the out-of-sample estimator, and 99.99% for Kolmogorov-Smirnov with estimated parameters,

> mean(P.KS.Norm.out.of.sample>.05)
[1] 0.85563
> mean(P.KS.Norm.estimated.param>.05)
[1] 0.99984
> mean(P.Lilliefors.Norm>.05)
[1] 0.9489

So using Kolmogorov-Smirnov with estimated parameters is not good, since we might accept  too often. On the other hand, if we use this technique with two samples (one to estimate parameter, one to run goodness of fit tests), it looks much better ! even if we reject  too often. For one test, the rate of first type error is rather large, but for the other, it is the rate of second type error...

Tuesday, September 18 2012

Open question on tail index

A short post today on a rather surprising result. I can't figure out how to prove it (I don not know if the result is valid, or known, so I claim it is an open question). Consider an i.i.d. sequence http://latex.codecogs.com/gif.latex?\{X_1,\cdots,X_n\} with common survival function http://latex.codecogs.com/gif.latex?\overline{F}. Assume further that http://latex.codecogs.com/gif.latex?\overline{F} is regularly varying at http://latex.codecogs.com/gif.latex?\infty, written http://latex.codecogs.com/gif.latex?\overline{F}\in%20RV_{\alpha} , i.e.

http://latex.codecogs.com/gif.latex?\lim_{t\rightarrow\infty}\frac{\overline{F}(tx)}{\overline{F}(t)}=x^{\alpha}

with http://latex.codecogs.com/gif.latex?\alpha\in(-\infty,0). Then the tail of http://latex.codecogs.com/gif.latex?X is Pareto type, and we can use Hill's estimator to estimate http://latex.codecogs.com/gif.latex?\alpha. An interesting result is that if http://latex.codecogs.com/gif.latex?\overline{F}\in%20RV_{\alpha} and http://latex.codecogs.com/gif.latex?\overline{G}\in%20RV_{\beta}, then http://latex.codecogs.com/gif.latex?\overline{F}\star\overline{G} is also regularly varying, where http://latex.codecogs.com/gif.latex?\star denote the convolution operator. More precisely, when http://latex.codecogs.com/gif.latex?\alpha=\beta, then (see Feller (), section VIII.8), then http://latex.codecogs.com/gif.latex?\overline{F}\star\overline{G} is regularly varying with the same index.

This can be visualized below, on simulation of Pareto variables, where the tail index is estimated using Hill's estimator. First, we start we one sample, either of size 20 (on the left) or 100 (on the right),

> library(evir)
> n=20
> set.seed(1)
> alpha=1.5
> X=runif(n)^(-1/alpha)
> hill(X)
> abline(h=1.5,col="blue")

If we generate two (independent) samples, and then look at the sum, Hill's estimator does not perform very well (the sum of two independent Pareto variates is no longer Pareto, but only Pareto type) we have

> set.seed(1)
> alpha=1.5
> X=runif(n)^(-1/alpha)
> Y=runif(n)^(-1/alpha)
> hill(X+Y)
> abline(h=1.5,col="blue")

The idea is then to use a Jackknife strategy in order to (artificially) increase the size of our sample. Thus, consider sums on all pairs of all http://latex.codecogs.com/gif.latex?X_i's, i.e. sample

http://latex.codecogs.com/gif.latex?\{X_1+X_2,\cdots,X_1+X_{n},X_2+X_3,\cdots,X_{n-1}+X_n\}

Let us use Hill's estimator on this (much larger) sample.

> XC=NA
> for(i in 1:(n-1)){
+ for(j in (i+1):n){
+ XC=c(XC,X[i]+X[j])
+ }}
> XC=XC[-1]
> hill(XC)
> abline(h=1.5,col="blue")
> abline(h=1.5*2,col="blue",lty=2)
> abline(h=1.5^2,col="blue",lty=3)

Here, with 20 observations from the initial sample, it looks like the tail index is http://latex.codecogs.com/gif.latex?2\alpha (on the left). With 100 observations, it looks like it is http://latex.codecogs.com/gif.latex?\alpha^2 (on the right). On the graphs above, I have plotted those two horizontal lines. It's odd, isn't it ? Of course, I did not really expect http://latex.codecogs.com/gif.latex?\alpha since we do not have an i.i.d. sample. Identically distributed yes, with a regularly varying survival function from what we've mentioned before. But clearly not independent. So Hill's estimator should not be a good estimator of http://latex.codecogs.com/gif.latex?\alpha, but it might be a good estimator of some function of http://latex.codecogs.com/gif.latex?\alpha

If we go one step further, an consider all triplets,

http://latex.codecogs.com/gif.latex?\{X_1+X_2+X_3,\cdots,X_1+X_2+X_{n},X_2+X_3+X_4,\cdots,X_{n-2}+X_{n-1}+X_n\}

(observe that now, the sample size is huge, even if we start with only 20 points). Then, it looks like the tail index should be http://latex.codecogs.com/gif.latex?\alpha^3 (at least, we can observe a step of Hill's plot around http://latex.codecogs.com/gif.latex?\alpha^3). 

> XC=NA
> for(i in 1:(n-2)){
+ for(j in (i+1):(n-1)){
+ for(k in (j+1):n){
+ XC=c(XC,X[i]+X[j]+X[k])
+ }}}
> XC=XC[-1]
> hill(XC)
> abline(h=1.5,col="blue")
> abline(h=1.5*3,col="blue",lty=2)
> abline(h=1.5^3,col="blue",lty=3)

My open question is whether there is general result behind, or if it was just a coincidence that those values appear so clearly.

Monday, June 18 2012

Date of death, birthday and Elvis Presley

10 days ago, a study published on http://www.annalsofepidemiology.org/ mentioned that "Death has a preference for birthdays" (as claimed in the title). The conclusion of the paper is that, in general, birthdays do not evoke a postponement mechanism but appear to end up in a lethal way more frequently than expected (“anniversary reaction”). Well, this is not new, and several previous articles have mentioned that point, e.g. Angermeyer et al. (1987).

I found the idea interesting since in demography, there is a large literature trying to extrapolate death rates from discrete to continuous time. Extrapolation are usually extremely smooth. But none of them integrate that aspect of mortality precisely on the birthday. The problem is that it is rather difficult to say something since datasets with individual observations are rare, online.

But yesterday, @coulmont sent me a tweet mentioning a website. I do not know if this is legal (even if some explanations are given), but I will mention courtesy of http://ssdmf.info/. It is a so-called Social Security Death Master File, containing individual informations about deaths in the US, as well as geographic information (as described on http://www.ssa.gov/), for people having a social security number.

With R, it is possible to work on those files (even they are huge, with tens of millions observations). For instance, we can check who is inside.

> elvis=scan("ssdm2",skip=22371720,n=1,what="character",sep=",")
> elvis
[1] " 409522002PRESLEY         ELVIS     0800197701081935  " 

If you believe that Elvis is dead, you might agree that this database can be accurate (or at least, not too bad). And further, we can see here how to read the result: Elvis was born on January 8, 1935 (8 last digits), and died on August 16, 1977 (8 digits before). Obviously here, there are some problems with the dataset (we do not have the day of the death of Elvis). So here, we remove all the observations that do not give us proper dates. Then, the idea is to assume that the person died in 2000 (or any year since the point is to focus on days and months). Then, we count the number of days between the day of death and the birthday in 2001 (that would have been after) and the one in 2000 (that was either before or after the death), so that we can derive the number of days after the birthday,

dates=substr(base,66,81)
death=as.Date(substr(dates,1,8),"%m%d%Y")
birth=as.Date(substr(dates,9,16),"%m%d%Y")
indice=is.na(death)|is.na(birth)
mean(indice)
mdeath=substr(dates,1,2)
ddeath=substr(dates,3,4)
mbirth=substr(dates,9,10)
dbirth=substr(dates,11,12)
indice=which(ddeath!="00")
birth1=as.Date(paste(mbirth[indice],
dbirth[indice],"2000",sep=""),"%m%d%Y")
birth2=as.Date(paste(mbirth[indice],
dbirth[indice],"2001",sep=""),"%m%d%Y")
death=as.Date(paste(mdeath[indice],ddeath[indice],
"2000",sep=""),"%m%d%Y")
k=length(indice)
diffday=cbind((as.numeric(death-birth1))[1:k],
(as.numeric(death-birth2))[1:k])
DIFF=apply(diffday,1,function(x) {min(x[x>=0])})

What we have here is the number of days following the previous birthday. If we look at the distribution of that number of days, we obtain

counts=table(DIFF)
plot(as.numeric(names(counts)),
as.numeric(counts))
counts["0"]/(mean(counts[100:200]))
> counts["0"]/(mean(counts[100:200]))
0
1.121261 

Thus, the death excess on the day of birth was around 12%, which is rather close to the one obtained from the Swiss mortality statistics 1969–2008 (in Ajdacic-Gross et al. (2012)). Note that here, we just play with a small subset of the entire dataset,

That database is probably extremely interesting, except that it suffers a huge selection bias, since only dead people are in that database. So it might be useless if we wish to study life expectancy of people named Bill versus people named Georges (that was something I wanted to investigate initially). But we'll see what else we can do with it (since Ewen have been able to write some code to go through that huge dataset).

Thursday, June 7 2012

La pyramide des âges, au Canada

L'autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la "pyramide" des âges (évoquée par exemple sur http://www.lapresse.ca/). Bon, j'ai des données assez proches, tirées de http://www.mortality.org/).

> pop=read.table(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,]

On peut faire assez facilement des pyramides des ages, par année, avec le code suivant,

> library(plotrix)
> seuils=seq(0,110,by=10) > pop$tranche=cut(pop$Age,seuils, right = FALSE) > an=1955 > base=pop[pop$Year==an,] > women=aggregate(base$Female, + list(base$tranche),sum)[,2] > men=aggregate(base$Female, + list(base$tranche),sum)[,2] > nom=as.character(unique(pop$tranche)) > pyramid.plot(men/sum(men)*100, + women/sum(women)*100,labels=nom,gap=2, + lxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red"), + rxcol=c("blue","blue","purple","purple","purple", + "purple","red","red","red","red","red"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

Au lieu de faire des barres horizontales, on peut cumuler, par tranche d'age, et regarder l'évolution dans le temps des proportions "moins de 20 ans" (bleu), "entre 20 et 60 ans" (mauve) et "plus de 60 ans" (rouge).

> seuils=c(0,20,60,110)
> pop$tranche2=cut(pop$Age,seuils, right = FALSE)
> YEAR=1921:2008
> totaux=matrix(NA,length(YEAR),3)
> for(i in 1:length(YEAR)){
+ base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total,
+ list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu'après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l'isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d'interprétation). A suivre...

Monday, May 21 2012

Présence policière et arrestations

Demain, c'est la grosse manifestation prévue à Montréal, (comme déjà le 22 mars, et le 22 avril). Avec les derniers débordements policiers de ces derniers soirs, il risque d'y avoir du monde. A propos de manifestations et de service d'ordre, une chroniqueur montréalais (je n'ai pu retrouvé la source) disais qu'il était surpris du nombre (important) d'arrestations lors de manifestations, et faisait un parallèle avec les grandes manifestations des années 70, beaucoup moins violentes, qui avaient donné lieu à beaucoup moins d'arrestations. 

En France, les services de communication de la Préfecture de Police de la ville de Paris ont eu la gentillesse de bien vouloir partager quelques statistiques sur les manifestations, et les interpellations. En particulier, un tableau reprend deux informations, avec le nombre de policiers mobilisés (correspondant à une participation aux évènements de voie publique encadrés par la DOPC) et le nombre d'interpellations,

Si on représente graphique ces deux séries (que ne peuvent être normées par le nombre de manifestations, à cause du problème d’agrégation annuelle... je renvois ici, sur la largeur des trottoirs en France, ou pour un point méthodologique, sur le comptage des manifestants)

Je trouve dommage de ne pas avoir d'historique plus long mais je vais continuer de creuser. En fait, au delà de ce rapide graphique (qui ne me permet pas franchement de faire d'analyse statistique intéressante), cette recherche de données m'a permis d'entrer en contact avec bon nombre de chercheurs, dont Fabien Jobard qui m'a offert une mise en perspective passionnante (que je vais essayer de remettre en forme ici: les erreurs de retranscription sont ici miennes).

Car la question me semblait simple: peut-on avoir le nombre de personnes arrêtées, ou interpellées, lors de manifestations ? Tout d'abord, il faudrait qu'il existe une nomenclature policière qui permette d'isoler des faits de violence en manifestation et, comme le notait Fabien Jobard, pour tout un tas de raison, cette nomenclature n'existe pas: on peut être arrêté pour violence sur policier, destructions, dégradations, violences ou voies de faits diverses... mais la conséquence pour les statisticiens est que chacun de ces actes est noyé dans un ensemble (trop) vaste qui ne permet pas d'isoler les actes spécifiques aux manifestations. Au delà de ces statistiques policières, il semble existe en France une statistique judiciaire, beaucoup plus précise, mais qui, selon Fabien, souffre des biais de sélection habituels de toute statistique judiciaire. Enfin, et c'est peut être le plus intéressant, c'est qu'il est délicat de distingue ce que les policiers appellent, à l'issue d'une manifestation, des "interpellations", entre celles qui en sont vraiment, et celles qui n'en sont pas. Pour reprendre des exemples donnés par Fabien Jobard, on peut imaginer que les policiers attrapent une personne, lui font une remontrance (j'utilise un vocable plus doux que celui évoqué dans le courriel de Fabien), et le relâche. Il peuvent aussi l’arrêter, le conduire devant un officier de police judiciaire, qui le relâche tout de suite (avec ou sans la remontrance mentionnée auparavant). Efin, il est possible d'attraper un individu suspect, de le conduire devant un officier de police judiciaire, qui transmet au Parquet (dans ce cas en revanche, il conviendra d'éviter les écarts de conduite). Seul ce dernier acte est une "interpellation". Pour ceux qui souhaitent aller plus loin, je renvois à vers Jobard (2010), qui a analysé les manifestations de 2005.

En vertu de la loi 78 en vigueur depuis la fin de semaine, je tiens à rappeler que je n'appelle personne à manifester. Mais de la même manière que j'encourage toute personne qui manipule des sondages à en faire passer, je vais essayer de me rapprocher (à nouveau) de la SPVM afin de comprendre ce qui se passe à Montréal. A suivre très probablement dans les jours à venir...

Sunday, May 20 2012

Qui ne voit pas la présence du diable dans la loi ?

Je suis de plus en plus dérouté par le mélange des genres dans l'espace médiatique, de ces prétendus spécialistes qui nous expliquent des choses qu'ils ont du découvrir 10 minutes avant de commencer à taper un article. Ce matin, des éditorialistes se sentaient autoriser à fabuler sur un sondage conduit sur internet auprès de 800 personnes, et j'avoue que cela m'a agacé... J'ignore si j'ai davantage de légitimité, mais je me suis dit que la loi 78 adoptée vendredi (au Québec) pourrait être un prétexte intéressant pour écrire sur un sujet d'actualité (chose que j'hésite d'ordinaire à faire, car ce n'est pas mon métier, loin de là). Pourquoi ne pourrais-je pas, moi aussi, faire mon malin et analyser un texte de loi ? On va donc regarder cette loi dans le détail, à partir du texte en ligne sur http://profscontrelahausse.org/ (ou pour être complètement honnête, le début seulement de la loi)

L=paste("Dans la presente loi, a moins",
"que le contexte nindique",
"un sens different, on entend par:",
"association d’etudiants}}:une",
"association ou un regroupement ",
"d’associations de niveau postsecondaire",

etc

"la qualite de l’enseignement, les ",
"services requis de façon a ",
"tenir compte des circonstances particulieres",
"resultant de l’interruption ",
"de la session d’hiver de l’annee 2012 ou ",
"de la session d’ete de l’annee 2012.")
 
L2=substring(toupper(L),1:nchar(L),1:nchar(L))
alphabet=c("A","B","C","D","E","F","G"
,"H","I","J","K","L","M","N","O",
"P","Q","R","S","T","U","V","W","X","Y","Z")
L2=L2[L2%in%alphabet]
ML=t(matrix(L2,nc,nl))

On peut visualier le début du texte ci-dessous (l'idée est de ne garder que les lettres, pour faire simple, mais on peut tenter plus largement en gardant la ponctuation, et les chiffres, ça ne change rien aux algorithmes que nous allons développer)

Commencons pas des trivialités: le texte de cette loi est fondamentalement inégalitaire: ce sont toujours les mêmes qui sortent gagnants, et les mêmes qui sortent perdants,

> sort(table(L))[c(1,2,3,4,5,20,21,22,23)]
L
Y   J   X   H   F   N   A   S   E
7   8  10  24  26 355 386 402 820 
On retrouve ici reproduites les inégalités observées dans presque tous les textes (sauf peut être La Disparition, mais c'est une autre histoire): N, A, S et E en tête, et toujours les mêmes petites lettres qui se trouvent honteusement oubliées...  
Maintenant, si on prend le temps de regarder en détail le contenu de la loi, force est de constater qu'on y voit la présence du Malin constamment ! 
Pour reprendre une idée de Bahye ben Asher ibn Halawa (רבינו בחיי) on peut tenter d'utiliser le code de la ThorahJulien Prévieux avait utilisé cette technique il y a quelques temps pour retrouver toute une terminologie associée aux crashs boursiers dans une page du Capital de Karl Marx, prise au hasard, ou aux paradis fiscaux dans la Richesse des Nations d'Adam Smith. Les amateurs de la théorie de Ramsey prétendront que c'est une supercherie, et ils auront probablement raison. Mais après avoir joué l'autre jour sur la recherche de schémas graphiques dans des images (en l’occurrence des rayures rouges et blanches, afin de repérer Charlie), j'ai eu envie de jouer à chercher des mots dans une matrice de lettres.
Commencons par nous donner un mot, au hasard, e.g.`
MOT="DEMON"
Le but est de voir si ce mot figure, ou pas, dans la matrice de lettre, suivant un schéma simple (c'est l'idée du code de la Thorah). Ma stratégie est ici assez simple: il faut commencer par chercher la lettre la moins présente dans la matrice, qui figure dans le mot. En fait, dans DEMON, je vais me focaliser sur les lettres du centre (en excluant la première et la dernière), i.e. DEMON
lettrelaplusrare=function(mot=MOT){
lettres=substring(mot,2:(nchar(mot)-1),
2:(nchar(mot)-1));
fL=TL[lettres];
return(names(fL)[which.min(fL)])}
rare=lettrelaplusrare()
On trouve ainsi la lettre qui va servir de base dans notre algorithme. A la rigueur, je peux faire une petite fonction qui va colorer une lettre dans la matrice,
dessinlettre=function(ltr="A",clr="yellow"){
which(L==ltr);
CLR=rep(NA,length(ML));
CLR[which(as.vector(ML)==ltr)]=clr;
dessin(ML,CLR,plt=FALSE)}
mais pour le moment, ça ne sert à rien. Dans un second temps, on regarde les lettres qui entourent la lettre la plus rare (c'est pour ça que j'ai exclus les lettres au bord, afin d’être certain d'avoir toujours une lettre avant, et une après),
lettres=substring(MOT,1:(nchar(mot)),
1:(nchar(mot)))
avant=lettres[which(lettres==rare)-1]
apres=lettres[which(lettres==rare)+1]
avant=avant[1]
apres=apres[1]
dessinlettre(rare,"yellow")
dessinlettre(avant,"green")
dessinlettre(apres,"pink")
On a ainsi trois lettres importantes: une relativement rare au centre, en jaune, une avant, en vert, et une après en rose. C'est ce qu'on peut visualiser sur le dessin ci-dessous,
Bon, ensuite on rentre dans l'algorithme, un peu lourd à mon gout: on va étudier tous les triplets possibles de ces trois lettres, en identifiants ceux pour lesquels la lettre au centre est précisément au milieu, coincée entre les deux autres lettres. C'est le principe de ce code de la Thorah. La première étape est de récupérer les coordonnées (dans la matrice) de ces trois lettres,
lrr=which(as.vector(ML)==rare)
lav=which(as.vector(ML)==avant)
lap=which(as.vector(ML)==apres) 
LIGrr=(lrr-1)%%nrow(ML)+1
COLrr=(lrr-1)%/%nrow(ML)+1
LIGav=(lav-1)%%nrow(ML)+1
COLav=(lav-1)%/%nrow(ML)+1
LIGap=(lap-1)%%nrow(ML)+1
COLap=(lap-1)%/%nrow(ML)+1
ptrr=cbind(LIGrr,COLrr)
ptav=cbind(LIGav,COLav)
ptap=cbind(LIGap,COLap)
Ensuite, on balaye l'ensemble des triplets (et ça peut être long),
pointsalignes=function(x,y,z){
((y[1]+2*(x[1]-y[1]))==z[1])&((y[2]+
2*(x[2]-y[2]))==z[2])}
LEQUEL=rep(NA,3)
for(a in 1:length(LIGrr)){
for(b in 1:length(LIGav)){
for(c in 1:length(LIGap)){
if(pointsalignes(ptrr[a,],ptav[b,],ptap[c,])==
TRUE){LEQUEL=cbind(LEQUEL,c(a,b,c))}
}}}
Parmi les triplets, on retient seulement ceux pour lesquels la lettre du milieu est précisément au centre, entre les deux autres. Enfin arrive la dernière étape pénible: regarder suivant la suite logique ce que sont les mots que l'on obtient (une fois placées ces trois lettres, la récurrence est complètement décrite). On retient alors les quintuplets correspondant au mot DEMON,
V=LEQUEL
for(k in 2:ncol(V)){
CLR=rep(NA,length(ML));
i1=ptrr[V[1,k],];
CLR[(i1[2]-1)*ncol(ML)+i1[1]]=clr;
lettres=substring(MOT,1:(nchar(mot)),
1:(nchar(mot)))
milieu=which(lettres==rare)
obtenu=ML[i1[1],i1[2]]
for(u in 1:(milieu-1)){
i2=ptav[V[2,k],]+(u-1)*(ptav[V[2,k],]-ptrr[V[1,k],]);
indice=(i2[2]-1)*ncol(ML)+i2[1]
if((i2[2]<=ncol(ML))&(i2[1]>=1)&(i2[1]<=
nrow(ML))&(i2[1]>=1)){CLR[indice]=clr;
obtenu=paste(ML[i2[1],i2[2]],obtenu,sep="")}}
for(u in (milieu+1):nchar(MOT)){
i3=ptap[V[3,k],]+(u-milieu-1)*(ptap[V[3,k],]-
ptrr[V[1,k],]);
indice=(i3[2]-1)*ncol(ML)+i3[1]
if((i3[2]<=ncol(ML))&(i3[1]>=1)&(i3[1]<=
nrow(ML))&(i3[1]>=1)){CLR[indice]=clr;
obtenu=paste(obtenu,ML[i3[1],i3[2]],sep="")}}
if(obtenu==MOT)
{print(paste(k,"  ",obtenu))}}
Reste à faire tourner cette boucle. Et ca ne manque pas, DEMON apparaît effectivement dans le texte,
On peut aller vérifier, pour ceux qui douteraient,
Coïncidence ? Peut-être.... alors tentons un autre mot, comme MALIN, qui apparait aussi,


ou alors SATAN... ce dernier apparaît huit fois dans le début du texte !
Étonnant non.... en tous les cas, je me suis bien amusé à coder ce petit algorithme, et de voir qu'il fonctionnait aussi bien ! Les amateurs de probabilités pourront se lancer dans des calculs, ou consulter quelques articles qui évoquent cette technique, comme McKay,Bar-Natan, Bar-Hillel & Kalai (2001), et tous les articles évoqués dans les références... Quant au juge qui devra se prononcer sur ma condamnation pour blasphème et pour hérésie, il va de soi que ce billet est une blague. Le diable n'existe pas ! On n'est plus au moyen-age quand même....

Saturday, May 19 2012

Births and week-ends, in France

This week, I have seen on the internet (sorry, I cannot find proper references) the graph produced here on the right: which birthday is most likely ? The fact that I have no further information is important, since I do not know in which country such a graph was obtained. At least, I know it should not be France...

In France, I have already mentioned that there is a strong week-end effect: nowadays, there is 25% less deliveries during week-ends than during the week. Calot (1981) observed already that there were less deliveries on Sundays. This has been confirmed more recently, e.g. in http://www.lepoint.fr/ or http://www.prepabl.fr/, with a significant difference between week days, and week-ends. Here  is the number of birth per day, over 40 years, with in blue the average trend during the week, and in red, during week-ends,

naissance=read.table(
"http://freakonometrics.free.fr/naissanceFR2.txt")
attach(naissance)
date=as.Date(date)
plot(date, nbre,cex=.5)
t2=as.POSIXlt(date)
jour=t2$wday
X=naissance$date
Y=naissance$nbre
J=jour
df=data.frame(X,Y,J)
library(splines)
regs=lm(Y~bs(X,df=20),data=df[jour%in%c(0,6),])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="red",lwd=3)
regs=lm(Y~bs(X,df=20),data=df[jour%in%1:5,])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="blue",lwd=3)

If we look at the evolution of the ratio week-ends over weeks days, we have the following graph

t2=as.POSIXlt(date)
jour=t2$wday
jour=jour[1:(1982*7)]
nbre2=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbre2[1:5+7*(i-1)]=nbre[1:5+7*(i-1)]*taux
nbre2[6:7+7*(i-1)]=nbre[6:7+7*(i-1)]
nbre2[1:7+7*(i-1)]=
mean(nbre[1:7+7*(i-1)])/mean(nbre2[1:7+7*(i-1)])*
nbre2[1:7+7*(i-1)]
}
nbretaux=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbretaux[1:7+7*(i-1)]=taux
}
plot(date[1:length(nbre2)],nbretaux)
X= date[1:length(nbre2)]
Y=nbretaux
library(splines)
reg=lm(Y~bs(X,df=20))
Yp=predict(reg)
lines(X,Yp,col="red",lwd=3)

In the beginning of the 70's, during week-ends, there were 5% less deliveries, but 25% less around 2000. It is then possible to produce the same kind of graphs as the one above, per year of birth. And here, we clearly observe the importance of the week end effect (maybe also because of color choice)

naissance=read.csv(
"http://freakonometrics.free.fr/naissanceFR.csv",
sep=";")
M=as.matrix(naissance[,3:ncol(naissance)])
BIRTH=as.vector(t(M))
YEAR=rep(1968:2005,each=12*31)
MONTH=rep(rep(1:12,each=31),38)
DAY=rep(1:31,12*38)
X=NA
for(y in 1968:2005){
sbase=base[YEAR==y,]
X=c(X,sbase$BIRTH/sum(sbase$BIRTH,
na.rm=TRUE))
}
base=data.frame(YEAR,MONTH,DAY,
BIRTH,BIRTHDAYPROB=X[-1])
 
m1=min(base$BIRTHDAYPROB,na.rm=TRUE)
m2=max(base$BIRTHDAYPROB,na.rm=TRUE)
y=1980
colr=rev(heat.colors(100))
sbase=base[YEAR==y,]
plot(0:1,0:1,col="white",xlim=c(-1,12),
ylim=c(-31,1),axes=FALSE,xlab=
paste("Naissance en",y,sep=" "),ylab="")
for(x in 1:nrow(sbase)){
a=sbase$MONTH[x];b=sbase$DAY[x]
polygon(c(a-.9,a-.9,a-.1,a-.1),-c(b-.9,b-.1,
b-.1,b-.9),col=colr[(sbase$BIRTHDAYPROB[x]-m1)/
(m2-m1)*100],border=NA)
}
text((1:12)-.5,.5,c("J","F","M","A","M","J","J",
"A","S","O","N","D"),cex=.7)
text(-.5,-(1:31)+.5,1:31,cex=.7)

http://freakonometrics.free.fr/ANNIVFRANCE.gif

Friday, May 18 2012

Non transitivity of correlation for random vectors in dimension 3

Dependence in dimension 2 is difficult. But one has to admit that dimension 2 is way more simple than dimension 3 ! I recently rediscovered a nice paper, Langford, Schwertman & Owens (2001), on transitivity of the property of being positively correlated (which inspired the odd title of this post). And more recently, Castro Sotos, Vanhoof, Van Den Noortgate & Onghena (2001) conducted a study which confirmed that there are strong misconceptions of correlation (and I guess, not only because probabilistic reasoning is extremely weak, as mentioned in Stock & Gross (1989)) and association, or correlation (as already stated in Estapa & Bataneor (1996), or Batanero, Estepa, Godino and Green (1996)). My understanding is that is it possible to have almost anything... even counterintuitive results. For instance, if we want to mix independence and comonotonicity (i.e. perfect positive dependence), all the theorems you might think of should probably be incorrect. Consider the following result (based on some old examples I have been using in my courses 5 or 6 years ago, see e.g. here)

"If X and Y are comontonic, and if Y and Z are comonotonic, then X and Z are comonotonic"

Well, this result seems to be intuitive, and probably valid. But it is not. Consider the following triplet,

Projections on bivariate planes of the three dimensional vector are

Here, X and Y are comonotonic, so are Y and Z, but X and Z are independent... Weird, isn't it ? Another one ?

"If X and Y are comontonic, and if Y and Z are independent, then X and Z are independent"

Again, even if it is intuitive, it is not correct... Consider for instance the following 3 dimensional distribution,

Here, X and Y are comonotonic, while Y and Z are independent, but X here and Z are countercomonotonic (perfect negative dependence). It is also possible to consider the following distribution,e

that can be visualized below,

In that case, X and Y are comonotonic, while Y and Z are independent, but X here and Z are comonotonic (perfect positive dependence). So obviously, we should be able to construct any kind of counterexample, on any kind of result we might think as intuitive.

To be honest, the problem with intuition is that is usually comes from the Gaussian case, and from the perception that dependence is related to correlation. Pearson's linear correlation. Consider the case of a 3 dimensional random vector, with correlation matrix

http://freakonometrics.blog.free.fr/public/perso6/CORRMATRICE.gif

Given two pairs of correlations, http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif and http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif, what could we say about http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif ? For instance, the intuition is that if http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif and http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif are positive, then http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif is likely to be positive too (perhaps). The only property (at least the most important) we have on that correlation matrix is that it should be positive-semidefinite. So if we play on eigenvalues, it should be possible to derive inequalities satisfied by  http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif.

Langford, Schwertman & Owens (2001) claim (in Theorem 3) that correlations have to satisfy some property, like

http://freakonometrics.blog.free.fr/public/perso6/kendall1.gif

which is simply the fact that the determinant of the correlation matrix has to be positive, that property was already mentioned in Kendall (1948), as an exercise,

But is that a sufficient and necessary condition ? Since I am extremely lazy, let us run some numerical calculation to visualize possible values for http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif, as function of http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif and http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif. Consider the following code

U=seq(-1,1,by=.1)
V=seq(-1,1,by=.001)
FSUP=function(a,b){
DF=function(c){min(eigen(matrix
(c(1,a,b,a,1,c,b,c,1),3,3))$values)};
V[max(which(Vectorize(DF)(V)>0))]}
FINF=function(a,b){
DF=function(c){min(eigen(matrix(
c(1,a,b,a,1,c,b,c,1),3,3))$values)};
V[min(which(Vectorize(DF)(V)>0))]}
MSUP=outer(U,U,Vectorize(FSUP))
MINF=outer(U,U,Vectorize(FINF))
library(RColorBrewer)
clr=rev(brewer.pal(6, "RdBu"))
U=U[2:20]
MSUP=MSUP[2:20,2:20]
MINF=MINF[2:20,2:20]
persp(U,U,MSUP,col="green",shade=TRUE)
image(U,U,MSUP,breaks=((-3):3)/3,col=clr)
persp(U,U,MINF,col="green",shade=TRUE)
image(U,U,MINF,breaks=((-3):3)/3,col=clr)
Here, we can derive the lower and the upper bound for http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif, as function of http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif and http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif.

In the dark blue area, the bound for the correlation can be really low, while in the dark red, the bound is very high (either the lower bound on the left, or the upper bound on the right). Since it might be hard to read, it is possible to fix for instance http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif, and to derive bonds for http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif, as function of http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif.
V=seq(-1,1,by=.001)
U=seq(-1,1,by=.1)
U=U[2:(length(U)-1)]
V=V[2:(length(V)-1)]
U=c(-.9999,U,.9999)
V=c(-.99999,V,.99999)
FSUP=function(a){
DF=function(c){min(eigen(matrix(
c(1,a,-.7,a,1,c,-.7,c,1),3,3))$values)};
V[max(which(Vectorize(DF)(V)>0))]}
FINF=function(a){
DF=function(c){min(eigen(matrix(
c(1,a,-.7,a,1,c,-.7,c,1),3,3))$values)};
V[min(which(Vectorize(DF)(V)>0))]}
 
VS=Vectorize(FSUP)(U)
VI=Vectorize(FINF)(U)
plot(c(U,U),c(VS,VI),col="white")
polygon(c(U,rev(U)),c(VS,rev(VI)),
col="yellow",border=NA)
lines(U,VS,lwd=2,col="red")
lines(U,VI,lwd=2,col="red")
On the graph below, we have bound for a negative correlation for http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif (on the left, with -0.7) and a positive correlation for http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif (on the right, here +0.7),

We do observe here extremely nice ellipses... Consider the case of a null correlation http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif then the region for possible values for http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif and http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif is the unit circle.

The interpretation is that if http://freakonometrics.blog.free.fr/public/perso6/correl-b.gif is null, and so is http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif then http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif might take any value between -1 and 1 (under the assumption that marginal distribution allow such values, e.g. marginal Gaussian distributions). On the other hand if http://freakonometrics.blog.free.fr/public/perso6/correl-a.gif is either -1 or +1 (perfect negative/positive correlation) then http://freakonometrics.blog.free.fr/public/perso6/correl-c.gif has to be null...

Wednesday, April 25 2012

Short versus long papers, in academic journals

This Monday, during my talk on quantile regressions (at the Montreal R-meeting), we've seen how those nice graphs could be interpreted, with the evolution of the slope of the linear regression, as a function of the probability level. One illustration was on large hurricanes, from Elsner, Kossin & Jagger (2008). The other one was on birthweight, from Abrevaya (2001).

It is also to illustrate that technique to academic publication, e.g. the length of papers, over time. Actually, the data we can extract from Scopus are quite similar to the ones uses on hurricanes. For several journals, it is possible to look at the length of articles. Since Scopus is quite expensive ($60,000 per year for the campus, as far as remember, so I can imagine the penalty I might have to pay for sharing such a dataset)

base=read.table("/home/scopus.csv",
header=TRUE,sep=",")
pages=base$Page.end-base$Page.start
year=base$Year

Again, a first idea can be to look at boxplots, and regression on (nonparametric) quantiles, here for Econometrica,

boxplot(pages~as.factor(year),col="light blue")
Q=function(p=.9) as.vector(by(pages,as.factor(year),
function(x) quantile(x,p)))
u=1:16
points(u,Q(p),pch=19,col="blue")
abline(lm(Q(p)~u,weights=table(year)),lwd=2,col="blue")

Consider now (as in the slides in the previous post) a quantile regression (instead of a regression on quantiles), for instance in the Annals of Probability,

library(quantreg)
u=seq(.05,.95,by=.01)
coefstd=function(u) summary(rq(pages~year,
tau=u))$coefficients[,2]
coefest=function(u) summary(rq(pages~year,
tau=u))$coefficients[,1]
CS=Vectorize(coefstd)(u)
CE=Vectorize(coefest)(u)
k=2
plot(u,CE[k,],ylim=c(min(CE[k,]-2*CS[k,]),
max(CE[k,]+2*CS[k,])))
polygon(c(u,rev(u)),c(CE[k,]+1.96*CS[k,],
rev(CE[k,]-1.96*CS[k,])),
col="light green",border=NA)
lines(u,CE[k,],lwd=2,col="red")
abline(h=0)

We have the following slope, for the year, as a function of the probability level,

The slope is always positive, so size of papers is increasing with time, short and long papers. But the influence of time is much larger for long paper than short one: for short papers (lower decile) every year, the size keeps increasing, with one more page every three years. For long paper (upper decile), it is two more pages every three years.

If we look now at the Annals of Statistics, we have

and for the evolution of the slope of the quantile regression,

Again the impact is positive: papers are longer in 2010 than 15 years ago. But the trend is the reverse: short papers (lower decile) are much longer, almost one more page every year, with long paper increase only by one more page every two years... Initially, I want to run such a study on a much longer term, with quantile regressions and splines to see when there might have been a change, both in lower and upper tails. Unfortunately, as suggested by some colleagues, there might have been some changes in the format of the journal (columns, margins, fonts, etc). That's a shame, because I rediscover nice short papers of 5-10 pages published 20 or 30 years ago. They are nice to read (and also potentially interesting for a post on the blog). 5 pages, that's perfect, but 40 pages, that's way too long. I wonder if I am the only one having this feeling, missing those short but extremely interesting papers....

Tuesday, April 24 2012

Régression médiane et géométrie

Suite à mon exposé d'hier, Pierre me demandait d'argumenter un point que j'avais évoqué oralement sans le justifier: " une régression médiane passe forcément par deux points du nuage ". Par exemple,

x=c(1,2,3)
y=c(3,7,8)
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
library(quantreg)
abline(rq(y~x,tau=.5),col="red")

Essayons de le justifier... de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe "entre" les points (on admettra que la régression médiane sépare l'espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

c'est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c'est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue",lwd=2)
for(i in seq(-2,3,by=.5)){
abline(a=i,b=3.2,col="blue",lty=2)}

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d)(H)
plot(H,D,type="l")

Formellement, l'optimum est ici
> optimize(d,lower=-5,upper=5)
$minimum
[1] -0.2000091
 
$objective
[1] 2.200009

Bref, retenons cette courbe, et notons qu'elle passe par un des points,

Et c'est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l'on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu'il n'y a qu'un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s'arrêter dès qu'on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2,
col="blue",lwd=2) points(x[1],y[1],cex=1.8,lwd=2,col="blue") for(i in seq(-1,5,by=.25)){ abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d2)(H)
plot(H,D,type="l")

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
$minimum
[1] 2.500018
 
$objective
[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
h=optimize(d2,lower=-5,upper=5)$minimum
abline(a=(y[1]-h*x[1]),
b=h,col="purple",lwd=2)

Pour comprendre pourquoi, comparons les deux cas,

  • en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue

  • en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l'autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c'est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu'à l'atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c'est toujours pareil...

Saturday, April 21 2012

Quotas et sondages

Apres @Christian il y a quelques semaines (suite aux premiers billets ou je faisais un parallèle entre l'estimation d'un paramètre sur le simplexe et les sondages électoraux), un ancien collègue de Rennes m'a demande si je pouvais faire un billet sur les sondages par quotas (au lieu de ne parler que de sondages aléatoires, ou de faire comme si les sondages politiques l’étaient). La question de Benoit était sur une phrase dans les transparents en ligne a Lyon (je n'ai pas trouve l'auteur)

"Dans les échantillonnages non probabilistes, comme la méthode des quotas utilisée par les instituts de sondage, on ne peut pas estimer la marge d’erreur due à l’élaboration de l’échantillon et généraliser les conclusions à l’ensemble de la population. En toute rigueur, les instituts de sondage ne peuvent pas calculer des intervalles de confiance. Dans la pratique, on sait que la marge d’erreur dans un échantillonnage par quotas est d’environ 3% pour un échantillon de 1000 individus (comme pour un échantillonnage aléatoire)"

N'ayant pas forcement ☐ les compétences  le temps  le courage (coche la bonne réponse), je vais me contente de rappeler qu'il y a quelques documents en ligne qui évoquent les méthodes de quotas. Par exemple le document de Philippe Cibois, ou celui en ligne a l’université de Pau. Historiquement Jerzy Neyman evoquait la methode des quotas (representative survey) en 1934 dans un article très intéressant. Dans Les sondages, publie a l'issu des journées de statistiques organisées en 1986, un chapitre explique ce qu'est la methode des quotas. 

La référence en question étant

Une référence passionnante semblerait etre l'article de Jean Claude Deville A Theory of Quota Surveys paru en 1991 dans Survey Methodology (mais non accessible en ligne). Si quelqu'un a le temps de regarder tout ça, éventuellement en faisant des simulations, je suis preneur ! Sinon merci a @wulab pour la référence des Simpsons pour l'illustration.

Wednesday, April 18 2012

25% ? sérieusemement, vous y croyez ?

(pour reprendre les statistiques postées par @guybirenbaum l'autre jour)

Réfléchissons deux minutes (ou un peu plus, c'est ce qu'on a fait ce midi avec ). On vous interroge sur un sondage: si vous croyez aux sondages, vous le dites. Mais si vous n'y croyez pas... vous répondez n'importe quoi (sinon c'est que vous y croyez, non ?). Bon essayons de formaliser ça afin de mieux comprendre... Soit la réponse donnée, et la vérité (ce que pense vraiment la personne interrogée). On sait que

Mais ce qu'on veut calculer, c'est , la probabilité qu'une personne y croit vraiment aux sondages. Pour ça, on sait que

Pour le premier terme, une personne qui y croit le dira (sinon c'est qu'elle n'y croit pas), donc . Par contre, une personne qui n'y croit pas peut dire n'importe quoi. Disons que vaut . Dans ce cas

ou encore

Supposons par exemple qu'une petite proportion des personnes qui n'y croient pas disent y croire, disons 10%. Dans ce cas, la vraie probabilité qu'une personne croit au sondage serait plutôt 16%. Si était un peu plus grand (25%), on serait juste tombé sur la proportion des rigolos qui répondent n'importe quoi aux sondages, car personne n'y croit !

Friday, April 13 2012

Ex-æquo aux élections ?

La saison 1 de Mad Men se termine avec l'élection de 1960, Nixon contre Kennedy, qui s'est jouée dans un mouchoir de poche.

L'autre jour, Benoit Rittaud évoquait sur http://images.math.cnrs.fr/ la possibilité que des candidats soient ex-æquo dans une élection, et les aspects légaux (car visiblement la constitution française n'a pas vraiment prévu ce cas). Le plus intéressant aurait été - à mon avis (mais mon avis est fortement biaisé) - de calculer la probabilité qu'un tel évènement survienne. Par exemple, qu'à l'issu du premier tour, le deuxième et le troisième candidat soient ex-æquo, de telle sorte qu'on ne puisse déterminer qui ira au second tour. Justement, un des commentaire affirme "la probabilité qu’un tel événement arrive est environ l’inverse du nombre typique de voix séparant le deuxième homme (ou femme) du troisième". N'ayant réussi à démontrer ce résultat, je me suis dit qu'on pourrait se lancer dans des simulations pour quantifier cette probabilité.

Commençons par supposer qu'il y a quatre candidats (et la possibilité de s'abstenir). On suppose que le taux d'abstention est de 20%, et que les intentions de vote pour chacun des candidats soit de 20%. 

EXE=function(N=100,ns=1000000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("A","B","C","D","E"),size=N,prob=rep(.2,5),
replace=TRUE)
tb=table(M)
Ms=rev(sort(tb[-(names(tb)=="A")]))
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Le soucis est que faire plusieurs centaines de millions de tirages de plus de 40 millions de lois multinomiales pourrait prendre du temps (genre vraiment beaucoup). Une stratégie intuitive peut être de faire des simulations pour une population plus faible, puis d'extrapoler. 

E=function(N) EXE(N)
taille=c(20,50,100,200,500,1000,2000)
proba=Vectorize(E)(taille)
plot(taille,proba,type="b",log="xy")
 
base=data.frame(y=proba,x=taille)
reg=lm(log(y)~I(log(x)),data=base)
s2=summary(reg)$sigma
m=exp(predict(reg,newdata=
data.frame(x=41194689)))*exp(s2^2/2)

L'extrapolation linéaire n'étant peut-être pas parfaite, on peut tenter un ajustement quadratique (comme cela a été fait sur la figure ci-dessus)

reg2=lm(log(y)~I(log(x))+I(log(x)^2),data=base)
s2=summary(reg2)$sigma
m2=exp(predict(reg2,newdata=data.frame(x=41194689)))*exp(s2^2/2)

La différence de prédiction avec les deux modèles est significative comme on le voit sur le graphique. Toutefois, avoir des ex-æquo avec une très très grande population est tellement rare qu'il faut réellement beaucoup de simulations pour avoir un ordre de grandeur valide. L'extrapolation semble être une alternative valide ici (mais je suis ouvert à toute critique, ou discussion).

Maintenant que l'on a une méthodologie (ou presque), utilisons des ordres de grandeurs réalistes. Par exemple en se basant sur l'élection présidentielle de 2002. 41.19 millions d'électeurs, 28.4% d'abstentions, et pour les votes exprimés, un candidat en tête avec 19.9% des votes, un second à 16.8% et le troisième à 16.2% (et les autres plus loin). Les intentions de votes pour les principaux candidats figurent dans la loi multinomiale ci-dessous

EXE=function(N=100,ns=1000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("abs","JC","JMLP","LJ","FB","AL","JPC"
,"NM","autres"),size=N,prob=c(.284,.1988*.716,
.1686*.716,.1618*.716,.0684*.716,.0572*.716,
.0533*.716,.0525*.716,.1714104),replace=TRUE)
tb=table(M);
Ms=rev(sort(tb[-(names(tb)%in%c("abs","autres"))]));
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Si on en croit l'affirmation du commentaire, la probabilité d'avoir eu des ex-æquo (entre le second et le troisième candidat) serait de l'ordre de 

> 1/(41194689*(.1686*.716-.1618*.716))
[1] 4.985823e-06

Je dois avouer que les calculs n'ont rien donné (on verra tout à l'heure ce qui se passe si on change la définition d’ex-æquo). Numériquement, les calculs tournent encore (après 4 jours) sur le serveur... Par contre, si la population ne contenait que 10000 électeurs, la probabilité d'avoir des ex-æquo (entre le 2ème et le 3ème candidat) serait de l'ordre de 3 sur 10000 (avec un intervalle de confiance que je n'ai toujours pas calculé)

> EXE(10000,100000000)
[1] 3.245562e-05

A titre d'illustration, on peut regarder ce qui s'est passé, commune par commune en France (les données sont en ligne sur le site du Ministère de l'Intérieur). L'exercice est biaisé parce qu'il y a de l'hétérogénéité spatiale... mais ça permet toujours d'illustrer... Les nombres de voix obtenues, pour le second et le troisième candidat, par commune est le suivant

> election=read.table(
+ "/Users/UQAM/Documents/Pres2002Tour1.csv",sep=";")
+ election=election[-1,]
> VOIX=election[,seq(14,74,by=4)]
> DEUXIEME=TROISIEME=TOTAL=rep(NA,nrow(VOIX))
> for(i in 1:nrow(VOIX)){
+ u=rev(sort(as.numeric(as.character(
+ election[i,seq(14,74,by=4)]))))
+ TOTAL[i]=sum(u)
+ DEUXIEME[i]=u[2]
+ TROISIEME[i]=u[3]
+ }

Sur 37000 communes,

> nrow(election)
[1] 37513

ou plutôt sur les 24000 communes qui ont eu plus de 5000 votes exprimés

> sum((TOTAL>5000))
[1] 24230

il y a eu 56 cas avec de deuxième et troisième candidats ex-aequo,

> base=data.frame(TOTAL,DEUXIEME,TROISIEME,COMMMUNE=election[,3])
> selection=(TOTAL>5000)
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 56

On a même le détail,

> sousbase[indice,]
TOTAL DEUXIEME TROISIEME                     COMMMUNE
1134   5819      699       699                        Crouy
1507   5386      717       717          Louroux-Bourbonnais
2603   6235     1151      1151                  Saint-Peray
2916   6346      935       935                        Coucy
2970   7529     1413      1413                   Haraucourt
5089   5567      717       717                      Vignats

(etc). On peut alors regarder, en fonction de la taille des communes (ici le nombre de votes exprimés) la probabilité d'avoir un ex-æquo, dans une commune donnée. On va subdivisionner en 50 classes (pour commencer), et compter à chaque fois le nombre d’ex-æquo,

> n=50
> niveau=seq(1/(2*n),1-1/(2*n),by=1/n)
> base$CLASSE=cut(TOTAL,quantile(TOTAL,seq(0,1,by=1/n)),
+ labels=niveau)
> base$EXAEQUO=(base$DEUXIEME==base$TROISIEME)*1
> m=by(base$EXAEQUO, base$CLASSE, sum)

Par exemple pour les 2% des plus grosses villes, on a eu 5 cas de communes avec des ex-aequos 

------------------------------------------------------------
base$CLASSE: 0.99
[1] 5

Si on regarde plus en détails, ce sont les communes avec plus de 7500 votes exprimés,

> quantile(TOTAL,.98)
98%
7427.76
> selection=(TOTAL>quantile(TOTAL,.98))
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 5
> indice=which(sousbase$DEUXIEME==sousbase$TROISIEME)
> sousbase[indice,]
TOTAL DEUXIEME TROISIEME       COMMMUNE CLASSE EXAEQUO
2970   7529     1413      1413     Haraucourt   0.99       1
14359  7496     1392      1392          Renac   0.99       1
25028  8104     1381      1381         Heloup   0.99       1
30014  8087     1340      1340         Barnay   0.99       1
36770  7435      764       764 Pont-sur-Yonne   0.99       1

On peut représenter ces nombres de cas d'ex-æquo et faire une régression de Poisson (car on fait du comptage ici),

> taille=length(TOTAL)/n
> plot(quantile(TOTAL,niveau),m/taille)
> sb=data.frame(niveau,exaequo=as.vector(m))
> library(splines)
> reg=glm(exaequo~bs(niveau),family=poisson,data=sb)
> u=seq(0,1,by=.001)
> lines(quantile(TOTAL,u),predict(reg,
+ newdata=data.frame(niveau=u), + type="response")/taille,lwd=2,col="red")

pour des facilités de lecture, on a en ordonnées des probabilités et en abscisse les tailles moyennes des villes

On peut faire la même chose avec 200 classes,

Bref, pour une commune de taille moyenne (disons moyenne supérieure, avec 5000 votes exprimés), on a 2 chances sur 1000 d'avoir un second et un troisième candidat ex-æquo. Ce qui laisserait entendre que pour 40 millions, la probabilité devrait être faible, voire très faible...

Bon, au delà du problème conceptuel des ex-æquo (ou légal, car je doute qu'il soit possible d'avoir une égalité parfaite, le Ministère de l'Intérieur ferait un recompte et trouverait forcément un résultat différent), il semble que la probabilité semble beaucoup plus faible que celle obtenue par la règle ad-hoc....

Mais peut-être serait-il réaliste de dire que l'élection n'a pas su dégager une majorité (c'est le but d'une élection il me semble) si le nombre de voix séparant deux candidats est trop faible. Par exemple, si moins de 100 voix séparent les candidats, on refait un vote ! Pour revenir sur mon exemple précédant, avec 10000 votants

> EXE100(10000,1000000)
[1] 0.012641

soit une chance sur 100 (je me suis autorisé à faire moins de simulations car l'évènement me semble moins improbable... donc il y a besoin de moins de simulations pour l'observer à l'occasion). Si on va plus loin, et que l'on fait des simulations que l'on projette, comme auparavant, pour estimer la probabilité d'avoir des ex-æquo non pas pour une ville, mais au niveau de la France entière, on est plutôt aux alentours de

> (m=exp(predict(reg,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
1.505422e-12 

pour un modèle linéaire, et

> (m=exp(predict(reg2,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
4.074588e-98  
pour un modèle quadratique (avec la confiance - relativement faible - que l'on peut accorder à une telle projection). Ce qui donne une différence énorme... Si on visualise tout ça,
lines(c(1,100000000),exp(predict(reg,newdata=
data.frame(x=c(1,100000000))))*exp
(s2^2/2),lty=2,col="light blue")

que l'on représente en bleu, en vert un ajustement quadratique

lines(10^(seq(1,8,by=.1)),exp(predict(reg,newdata=
data.frame(x=10^(seq(1,8,by=.1)))))*exp (s2^2/2),lty=2,col="dark green")

alors que le calcul suggéré dans le commentaire donnerait la courbe rouge,

lines(c(1,100000000),1/(c(1,100000000)*
(.1686*.716-.1618*.716)),col="red",lty=2)

Peut-on parler d'évènement improbable ? A titre de comparaison, quand on parle de quelque chose d'improbable, ou qui arrive avec une probabilité infinitésimale, on pense au loto. A l'ancien loto (français), il fallait choisir 6 numéros parmi 49 (sans remise). La probabilité de gagner le gros lot était de 

Pour l'Euromillions, il faut choisir 5 numéros parmi 50 numéros auxquels il faut ajouter les combinaisons possibles avec les Étoiles, soit 2 numéros à choisir parmi 9. La probabilité de gagner le gros lot est alors de

Bref, tomber sur des ex-æquo lors d'une élection nationale, en France, c'est a priori plus rare que gagner le gros lot en jouant une fois au loto... On peut donc dire qu'il ne devrait pas y avoir d’ex-æquo à l'issu du premier tour !

- page 1 of 10