Freakonometrics

To content | To menu | To search

Miscellaneous

Entries feed - Comments feed

Sunday, April 1 2012

Sunday evening, stupid games...

This evening, while I was about to wash the dishes, I heard my elders starting a game (call them Him and Her)
Him: "I have picked - in my head - a number, lower than 50. Try to guess..."
Her: "No way, too difficult..."
Him: "You can try five different numbers..."
Her: ".,. um ... No, no way..."
Me: "Wait... each time we suggest a number, you tell us if yours is either above, or below ?"
You can see me coming clearly, can't you ? Using a simple subdivision rule, we have a fast algorithm (and indeed, if I have to choose between washing the dishes and playing with the kids...)
Him: "um.... ok"
Her: "Daddy, are you sure we will win ?"
Me: "Well... I cannot promise that we will win... but I am rather sure [sic] that we will win quite frequently: more gains than losses..." (I guess).
Her: "Great ! I am playing with daddy..."
Him: "um.... wait, is it one of you trick, again ? I don't to play anymore... Do you want to see the books we've chosen at the library ?"
Her: "Sure..."
Me: "What ? no one wants to see if I was right ? that we have indeed more than 50% chances to win..."
Him and her: "No !
The point of that story ? If we listen to kids, science will not go forward, trust me. But I am curious... I want to see if my intuition was correct. Actually, the intuition was based on the fact that
> 2^5
[1] 32 
> 2^6
[1] 64
so in 5 or 6 steps the algorithm of subdivision should converge. I guess... I mean, I do not know for sure, since 50 is not a power of 2, so it might be difficult, each time, to split in two: we have to deal only with integers here...
To be sure, let us substitute my laptop to my son... to pick up numbers, randomly (yes, sometimes I feel like I am Doctor Tenma, 天馬博士). The algorithm is simple: there are bounds, and at each stop I should suggest the middle of the interval. If the middle is not an integer, I suggest either the integer below or the integer above (with equal probabilities).
cutinhalf=function(a,b){
m=(a+b)/2
if(m %% 1 == 0){m=m}
if(m %% 1 != 0){m=sample(c(m-.5,m+.5),size=1)}
return(round(m))}
The following functions runs 10,000 simulations, and tells us how many times, out of 5 numbers suggested, we got the good one.
winning=function(lower=1,upper=50,tries=5,NS=100000){
SIM=rep(NA,NS)
for(simul in 1:NS){
interval=c(lower,upper)
(unknownnumber=sample(lower:upper,size=1))
success=FALSE
for(i in 1:tries){
picknumber=cutinhalf(interval[1],interval[2])
if(picknumber==unknownnumber){success=TRUE}
if(picknumber>unknownnumber){interval[2]=picknumber}
if(picknumber<unknownnumber){interval[1]=picknumber}
#print(c(unknownnumber,picknumber,success,interval))
};SIM[simul]=success};return(mean(SIM))}
It looks like the probability that we got the good number is higher than 60%,
> winning()
[1] 0.61801
Which is not bad. And if the upper limit was not 50, but something else, the probability of winning would have been the following.
VWN=function(n){winning(upper=n)}
V=Vectorize(VWN)(seq(25,100,by=5))
plot(seq(25,100,by=5),V,type="b",col="red",ylim=c(0,1))

Actually, after losing a couple of times, I am rather sure that my son would have to us that we can suggest only four numbers. In that case, the probability would have been close to 30%, as shown on the blue curve below (where four numbers only can be suggested)
Anyway, as intuited, with five possible suggestions, we were quite likely to win frequently. Actually with a probability of almost 2 out of 3...and 1 out of 3 if my son had decided to pick an number between 1 and 100, or only 4 possible suggestions... Those are quite large actually, when we think about it. It reminds me that McGyver story I mentioned a few months ago... Anyway, calculating probabilities is nice, but I still have to wash the dishes...

Thursday, January 19 2012

Tout le monde peut sortir vainqueur d'un pari ?

Petit paradoxe du jour (emprunté à http://www.futilitycloset.com/). Avec un collègue, suite à un pari stupide, nous nous sommes engagé à mettre une cravate, et celui qui a payé sa cravate le plus cher devra la donner à l'autre (ce que nous ne savions pas avant de choisir ladite cravate). A priori, nous avons autant de chance l'un que l'autre de gagner. Si je perds, je perds le prix de ma cravate, disons http://freakonometrics.blog.free.fr/public/perso5/cravate01.gif. Alors que si je gagne, je gagne une cravate qui vaut plus cher que la mienne, disons http://freakonometrics.blog.free.fr/public/perso5/corrrec-01.gif, avec http://freakonometrics.blog.free.fr/public/perso5/corrre-2.gif. Aussi, mon espérance de gain est

http://freakonometrics.blog.free.fr/public/perso5/cravate04.gif

aussi, en moyenne, je vais sortir gagnant de ce pari. Sauf que mon collègue peut faire le même raisonnement... Donc en moyenne, on y gagne tous... Étonnant, non ?

Sunday, December 18 2011

Combien de temps profite-t-on de ses grands parents ?

Ce Hier matin, je suis tombé un peu par hasard sur deux graphiques de l'INSEE (en France) avec l'age moyen des mères à l'accouchement, en fonction de rang de naissance de l'enfant, avec tout d'abord 1905-1965,

puis 1960-2000

Ces graphiques sont passionnants en soi - comme en ont témoigné pas mal de followers sur Twitter - mais ils m'ont fait m'interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n'a - malheureusement - pas accès au données complètes sur le site, mais on peut trouver d'autres données intéressantes (en l’occurrence l'age moyen à la naissance).

> agenaissance=read.table("http://freakonometrics.blog.free.fr/
public/data/agenaissance.csv"
,header=TRUE,sep=",") > agenaissance$Age=as.character(agenaissance$AGE) > agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+ + as.numeric(substr(agenaissance$Age,4,4))/10 > plot(agenaissance$ANNEE+.5,agenaissance$AGE, + type="l",lwd=2,col="blue")
Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus,

On peut alors aller en cran plus loin, en se demandant non pas quel était l'âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)
> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)-
+ agenaissance$AGE
> w=(trunc(agenaissance$NAIS.MERE-.5))
> rownames(agenaissance)=agenaissance$ANNEE
> a1=agenaissance[as.character(w),]$NAIS.MERE
> a2=agenaissance[as.character(w+1),]$NAIS.MERE
> p=agenaissance$NAIS.MERE-(w+.5)
> agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2
> agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5-
+ agenaissance$NAIS.GRD.MERE
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE
2000  2000 30.3 30,3     1970.2       1942.87        57.63
2001  2001 30.4 30,4     1971.1       1943.80        57.70
2002  2002 30.4 30,4     1972.1       1944.92        57.58
2003  2003 30.5 30,5     1973.0       1945.95        57.55
2004  2004 30.5 30,5     1974.0       1947.05        57.45
2005  2005 30.6 30,6     1974.9       1948.04        57.46
> plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE,
+ type="l",lwd=2,col="red")
Là encore, on peut visualiser l'âge de la grand-mère maternelle à la naissance

A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d'espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.
> Deces <- read.table("http://freakonometrics.free.fr/
Deces-France.txt",header=TRUE)
> Expo  <- read.table("http://freakonometrics.free.fr/
Exposures-France.txt",header=TRUE,skip=2)
> Deces$Age <- as.numeric(as.character(Deces$Age))
> Deces$Age[is.na(Deces$Age)] <- 110
> Expo$Age <- as.numeric(as.character(Expo$Age))
> Expo$Age[is.na(Expo$Age)] <- 110
>  library(forecast)
>  library(demography)
>  YEAR <- unique(Deces$Year);nC=length(YEAR)
>  AGE  <- unique(Deces$Age);nL=length(AGE)
>  MUF  <- matrix(Deces$Female/Expo$Female,nL,nC)
>  POPF <- matrix(Expo$Female,nL,nC)
>  BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE,
+ years=YEAR, type="mortality",
+  label="France", name="Femmes", lambda=1)
> LCF <- lca(BASEF)
> LCFf<-forecast(LCF,h=100)
> A <- LCF$ax
> B <- LCF$bx
> K1 <- LCF$kt
> K2 <- K1[length(K1)]+LCFf$kt.f$mean
> K <- c(K1,K2)
> MU <- matrix(NA,length(A),length(K))
> for(i in 1:length(A)){
+ for(j in 1:length(K)){
+ MU[i,j] <- exp(A[i]+B[i]*K[j]) }}
> esp.vie = function(xentier,T){
+ s <- seq(0,99-xentier-1)
+ MUd <- MU[xentier+1+s,T+s-1898]
+ Pxt <- cumprod(exp(-diag(MUd)))
+ ext <- sum(Pxt)
+ return(ext) }
> EVIE = function(x,T){
+ x1 <- trunc(x)
+ x2 <- x1+1
+ return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) }
> agenaissance$EV=NA
> for(i in 1:100){
+ t <- 2006-i
+ agenaissance$EV[agenaissance$ANNEE==t]=
+ EVIE(x=agenaissance$age.GRD.MERE[
+ agenaissance$ANNEE==t],t) }
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE       EV
2000 30.3 30,3     1970.2       1942.87        57.63 29.13876
2001 30.4 30,4     1971.1       1943.80        57.70 29.17047
2002 30.4 30,4     1972.1       1944.92        57.58 29.39027
2003 30.5 30,5     1973.0       1945.95        57.55 29.52041
2004 30.5 30,5     1974.0       1947.05        57.45 29.72511
2005 30.6 30,6     1974.9       1948.04        57.46 29.80398
Autrement dit, sur la dernière ligne, l'espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d'environ 29.80 ans. On peut alors visualiser non seulement l'âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle,
> plot(agenaissance$ANNEE+.5,agenaissance$EV,
+ type="l",lwd=2,col="purple")

On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s'est stabilisé à une trentaine d'années. On peut aussi continuer, et remonter d'un cran (en refaisant tourner le code avec quelques modifications): on a alors l'âge (moyen) de son arrière grand mère à la naissance,

et la durée de vie (résiduelle) des arrière grand mères

On manque ici de données, mais il semble que l'on profite - en moyenne - environ 5 ans de son arrière grand mère. Maintenant on peut aussi s'interroger sur les limites de cette étude rapide. En particulier, de même qu'il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d'enfants et de petits-enfants a un impact sur la durée de vie résiduelle d'une personne (ou si on peut supposer l'indépendance comme on l'a fait ici).

Saturday, December 17 2011

Le paradoxe de l'examen surprise...

La semaine dernière, Steve Landsburg reprenait (sur http://thebigquestions.com) sous une forme amusante un paradoxe que j'avais évoqué voilà deux ans sur le blog. Je ne résiste pas au plaisir de me replonger dans ce paradoxe...
Un professeur rédigeant son plan de cours veut mentionner un examen surprise, qui compterait pour 10% de la note finale.  Il annonce donc à ses étudiants, au premier cours, qu'un examen surprise (et un seul) se tiendrait à une date (surprise) parmi les 12 séances de cours. Mais comme le fait noter un élève, l'examen ne pourrait avoir lieu lors de la dernière séance (la 12ème), car cela signifierait qu'il n'y a pas eu d'examen les 11 premières séances, et donc l'examen surprise devrait forcément se tenir au 12ème cours (mais alors il ne serait plus surprise du tout). Bref, impossible d'avoir un examen surprise ce jour là. Mais dans ce cas, surenchérit une étudiante, il ne peut non plus avoir lieu à l'avant dernière séance (la 11ème), car là encore cela signifierait qu'il n'y a pas eu d'examen les 10 premières séances, et donc l'examen surprise devrait forcément au 11ème cours (car il ne peut pas avoir lieu au 12ème comme on vient de le voir, et là encore, il ne serait plus surprise du tout). Bref, de manière récurrente, on peut montrer que l'examen surprise a alors lieu forcément le premier jour. Mais là encore, il n'y a pas de surprise. Et donc impossible qu'il y ait un examen surprise si ce dernier a été annoncé. Moralité, pour faire un examen surprise, il ne faut surtout pas évoquer qu'il y en aura un. Mais comme on ne peut pas faire un examen non annoncé dans les modalités de contrôle des connaissances, il est impossible de faire des surprises aux étudiants à l'université... C'est triste, non ?

Tuesday, October 18 2011

Halloween is in two weeks

Looking for scary things before Halloween ? I mean like really scary ? Just go on http://wearethe99percent.tumblr.com/, some stories are amazing, and scary....