# Freakonometrics

## Tag - Rstats

Monday, June 4 2012

## Longevity and mortality dynamics with R

Following the previous post on life contingencies and actuarial models in life insurance, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The second part of the talk (on Actuarial models with R) will be dedicated to longevity and mortality. A complete set of slides can be downloaded from the blog, but again, only some part will be presented.

As mentioned earlier, the codes are from a book on actuarial science in R, written with Christophe Dutang and Vincent Goulet (so far in French) that should appear, some day... The code used in the slides above can be downloaded from here, and datasets are the following,

```> DEATH <- read.table(
+ "http://freakonometrics.free.fr/Deces-France.txt",
+ "http://freakonometrics.free.fr/Exposures-France.txt",

For additional resources, I will use Rob Hyndman's package on demography, Heather Turner and David Firth's package on generalized nonlinear models (e.g. the slides of the short course Heather gave in Rennes at the UseR! conference in 2009), as well as functions developed by JPMorgan's LifeMetrics (functions are  fully documented in the LifeMetrics Technical Document). All those functions can be obtained using

```> library(demography)
> library(gnm)
> source("http://freakonometrics.free.fr/fitModels.R")```

Friday, January 13 2012

## ACT6420 Méthodes de prévision, introduction à R

Cet après midi, j'assurerais la démonstration, en salle informatique à 14h30, PK S1 590, au sous sol (salle L). On fera une introduction à R, afin de présenter les rudiments du langage. Je mets en ligne quelques lignes de code,

```> X=cars\$speed
> Y=cars\$dist
> n=length(X)
> (mx=mean(X))
[1] 15.4
> sum(X)/n
[1] 15.4
> s=0
> for(i in 1:n){s=s+X[i]}
> print(s/n)
[1] 15.4
>
> my=mean(Y)
> (sx=sd(X))
[1] 5.287644
> sqrt(sum((X-mx)^2)/(n-1))
[1] 5.287644
> sy=sd(Y)
> (r=cor(X,Y))
[1] 0.8068949
> sum((X-mx)*(Y-my))/((n-1)*sx*sy)
[1] 0.8068949
>
> X<=15
[1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[9]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[17]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[25]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[33] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[41] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[49] FALSE FALSE
> sum(X<=15)
[1] 26
> X[X<=15]
[1]  4  4  7  7  8  9 10 10 10 11 11 12 12 12 12 13 13
[18] 13 13 14 14 14 14 15 15 15
> mean(Y[X<=15])
[1] 27.5
> mean(Y[X>15])
[1] 59.75
>
> (b1=r*sy/sx)
[1] 3.932409
> (b0=my-b1*mx)
[1] -17.57909```
On peut aussi faire quelques graphiques afin de visualiser la régression,
```> plot(X,Y)
> plot(X,Y,xlab="Vitesse de la voiture",
+ ylab="Distance de freinage",pch=4,
+ xlim=c(0,30),axes=FALSE)
> axis(1);axis(2)
> abline(a=b0,b=b1,col="red",lty=2)
> Yh=b0+b1*X
> points(X,Yh,pch=19,col="blue")```

Si on continue à créer et manipuler des objets (qui seront utiles pour analyser la régression),
```> (SCT=sum((Y-my)^2))
[1] 32538.98
> (SCR=sum((Y-Yh)^2))
[1] 11353.52
> (SCE=sum((Yh-my)^2))
[1] 21185.46
> SCR+SCE
[1] 32538.98
>
> (s=sqrt(sum((Y-Yh)^2)/(n-2)))
[1] 15.37959
>
> (sb1=s/(sx*sqrt(n-1)))
[1] 0.4155128
> t=(b1-0)/sb1
> 2*(1-pt(t,df=n-2))
[1] 1.489919e-12
>
> lm(Y~X)

Call:
lm(formula = Y ~ X)

Coefficients:
(Intercept)            X
-17.579        3.932

> regression=lm(dist~speed,data=cars)
> names(regression)
[1] "coefficients"  "residuals"     "effects"
[4] "rank"          "fitted.values" "assign"
[7] "qr"            "df.residual"   "xlevels"
[10] "call"          "terms"         "model"
> regression\$coefficient
(Intercept)       speed
-17.579095    3.932409
> summary(regression)

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
Min      1Q  Median      3Q     Max
-29.069  -9.525  -2.272   9.215  43.201

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791     6.7584  -2.601   0.0123 *
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511,	Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

> summary(regression)\$sigma
[1] 15.37959```
On approfondira ces codes cet après midi.

Wednesday, January 11 2012

## Construire une courbe ROC

Juste avant les vacances, Jean-Pierre Liégeois, un jeune lecteur du var, me demandais par courriel, "A partir d'une régression logistique(ou d'une matrice de confusion 2x2), comment programmer en R, un programme qui construit la courbe ROC associée". Avant d'aller plus loin (et de répondre a la question), je vais renvoyer vers un vieux billet sur les matrices de confusion. L’idée est que l'on suppose que l'on dispose d'un prédicteur d'une variable prenant des valeurs 0 et 1 (ou pour reprendre la terminologie classique "positif" et "négatif"), par exemple un modèle logistique. Formellement, pour l'ensemble de nos observations, on a une valeur observée  et (comme je l'expliquais dans un autre billet) et d'un score . Et c'est ce score qu'on va utiliser pour construire la courbe ROC. Ce score sera utilise pour prédire . La règle d'affectation est alors simple: on se fixe un seuil , et

• si , alors   est "positif"
• si , alors   est "négatif"
On peut alors construire une matrice dite de confusion, qui est simplement un table de contingence,
 valeur observée valeur prédite "positif" "négatif" "positif" TP FP "négatif" FN TN
où TP désigne les vrais positifs (true positive), TN  les vrais négatifs (true negative),FP désigne les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
Quid de la mise en œuvre sous R ? Commençons par générer des données, et estimons un modèle de régression.
```set.seed(1)
n=50
X=rnorm(n)
Y=rbinom(n,size=1,prob=
exp(2*X-1)/(1+exp(2*X-1)))
B=data.frame(Y,X)
reg=glm(Y~X,family=binomial,data=B)
S=predict(reg,type="response")```
On a maintenant notre observation (variable prenant les valeurs 0 ou 1) et notre score. On va ensuite pouvoir choisir plusieurs valeurs possible pour le seuil, et visualiser le taux de vrais positifs, en fonction du taux de faux positifs.
```plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP=sum((Ps==1)*(Y==0))/sum(Y==0)
TP=sum((Ps==1)*(Y==1))/sum(Y==1)
points(FP,TP,cex=.5,col="red")
}```
On a alors le graphique suivant,

Si on relit les points, on a alors la courbe ROC,
```FP=TP=rep(NA,101)
plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP[1+s*100]=sum((Ps==1)*(Y==0))/sum(Y==0)
TP[1+s*100]=sum((Ps==1)*(Y==1))/sum(Y==1)
}
lines(c(FP),c(TP),type="s",col="red")```

En fait, le code est assez simple, et il traîne dans différents packages de R, e.g.
```library(ROCR)
pred=prediction(P,Y)
perf=performance(pred,"tpr", "fpr")
plot(perf,colorize = TRUE)```

On peut aussi s'amuser a bootstrapper l’échantillon pour construire des intervalles de confiance, ou ajuster des modèles théoriques,
```library(verification)
roc.plot(Y,P, xlab = "False Positive Rate",
ylab = "True Positive Rate", main = "", CI = TRUE,
n.boot = 100, plot = "both", binormal = TRUE)```

ou encore (toujours avec des bornes de confiance obtenues par bootstrap)
```library(pROC)
PROC=plot.roc(Y,P,main="", percent=TRUE,
ci=TRUE)
SE=ci.se(PROC,specificities=seq(0, 100, 5))
plot(SE, type="shape", col="light blue")```

Tuesday, November 8 2011

## ACT2040: régression binomiale négative

La régression binomiale négative n'apparait pas dans les familles proposées sous R, mais il suffit d'utiliser un autre package pour la faire (voire plusieurs autres si on veut aller plus loin sur les classes de modèles sur les variables de comptage),

```> library(MASS)
> ? glm.nb```

Si on compare avec la régression de Poisson, on voit que les estimations sont très proches,

```>  regp = glm(nbre~ageconducteur+offset(exposition),
+  data=base,family=poisson)
>  regnb = glm.nb(nbre~ageconducteur+offset(exposition),
+  data=base)
>  coefficients(regp)
(Intercept) ageconducteur
-3.206981838  -0.006714151
>  coefficients(regnb)
(Intercept) ageconducteur
-3.214933688  -0.006614091 ```

On peut également comparer les prédictions, voire rajouter une transformation nonlinéaire de l'âge du conducteur (pour noter que changer de loi influence peu l'estimation, car par défaut le lien est un lien logarithmique, comme pour la régression de Poisson).

``` regf = glm(nbre~as.factor(ageconducteur)+offset(exposition),
data=base,family=poisson)
age=sort(unique(base\$ageconducteur))
Y1=predict(regf,newdata=data.frame(ageconducteur=age,
exposition=1),type="response")
Y2=predict(regp,newdata=data.frame(ageconducteur=age,
exposition=1),type="response")
Y3=predict(regnb,newdata=data.frame(ageconducteur=age,
exposition=1),type="response")
plot(age,Y1)
lines(age,Y2,col="red",lwd=3)
lines(age,Y3,col="blue",lwd=3)

library(splines)
regps = glm(nbre~bs(ageconducteur,5)+offset(exposition),
data=base,family=poisson)
regnbs = glm.nb(nbre~bs(ageconducteur,5)+offset(exposition),
data=base)
Y2s=predict(regps,newdata=data.frame(ageconducteur=age,
exposition=1),type="response")
Y3s=predict(regnbs,newdata=data.frame(ageconducteur=age,
exposition=1),type="response")
lines(age,Y2s,col="red",lwd=1)
lines(age,Y3s,col="blue",lwd=1)```

Friday, November 4 2011

## Confidence interval for predictions with GLMs

Consider a (simple) Poisson regression . Given a sample where , the goal is to derive a 95% confidence interval for given , where is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below
```> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1\$fit,pch=4,lwd=3)```
i.e.

Let denote the maximum likelihood estimator of . Then

where is Fisher information of (from standard maximum likelihood theory). Recall that

where computation of those values is based on the following calculations

In the case of the log-Poisson regression

Let us get back to our initial problem.
• confidence interval for the linear combination
A first idea to get a confidence interval for is to get a confidence interval for (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that
thus, an approximation for the variance matrix of will be based on , obtained by plugging estimators of the parameters.
Then, since as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
has a normal distribution, centered on , with variance where is the variance of . All those quantities can be easily computed. First, we can get the variance of the estimators
```> i1=sum(predict(reg,type="response"))
> i2=sum(cars\$speed*predict(reg,type="response"))
> i3=sum(cars\$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)```
Hence, if we compare with the output of the regression,
```> summary(reg)\$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05```
Based on those values, it is easy to derive the standard deviation for the linear combination,
```> x=30
> P2
\$fit
1
5.046034

\$se.fit
[1] 0.05747075

\$residual.scale
[1] 1

>
> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084```
And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval
```> segments(30,exp(P2\$fit-1.96*P2\$se.fit),
+ 30,exp(P2\$fit+1.96*P2\$se.fit),col="blue",lwd=3)```
Based on that technique, confidence intervals are not centered on the prediction, but who cares ?

• delta method
Actually, those who like to use "more or less" expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,
```> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)\$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,+ newdata=data.frame(speed=30))
> P1
\$fit
1
155.4048

\$se.fit
1
8.931232

\$residual.scale
[1] 1```
The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.
```> segments(30,P1\$fit-1.96*P1\$se.fit,30,
+ P1\$fit+1.96*P1\$se.fit,col="blue",lwd=3)```

Note that those quantities - obtained with two different approaches - are rather close here

```> exp(P2\$fit-1.96*P2\$se.fit)
1
138.8495
> P1\$fit-1.96*P1\$se.fit
1
137.8996
> exp(P2\$fit+1.96*P2\$se.fit)
1
173.9341
> P1\$fit+1.96*P1\$se.fit
1
172.9101 ```
• bootstrap techniques
And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,

Tuesday, October 25 2011

## Consecutive number and lottery

Recently, I have been reading odd things about strategies to win at the lottery. E.g.

or

I wrote something a long time ago, but maybe it would be better to write another post. First, it is easy to get data on the French lotteries, including draws, number of winners and gains,

```loto=read.table("http://freakonometrics.blog.free.fr/public/data/loto.csv",sep=";",header=TRUE)
balls=loto[,c("boule_1","boule_2","boule_3","boule_4","boule_5","boule_6")]
q=function(x){quantile(x,(0:5)/5)}
sortballs=balls
consec=balls[,-1]
sortconsec=consec
for(i in 1:nrow(balls)){sortballs[i,]=q(balls[i,])
consec[i,]=sortballs[i,2:6]-sortballs[i,1:5]
sortconsec[i,]=sort(consec[i,])}
winner1=loto[,"nombre_de_gagnant_au_rang1"]
gain1=as.numeric(as.character(loto[,"rapport_du_rang1"]))
winner2=loto[,"nombre_de_gagnant_au_rang2"]
gain2=as.numeric(as.character(loto[,"rapport_du_rang2"]))
winner3=loto[,"nombre_de_gagnant_au_rang3"]
gain3=as.numeric(as.character(loto[,"rapport_du_rang3"]))
winner4=loto[,"nombre_de_gagnant_au_rang4"]
gain4=as.numeric(as.character(loto[,"rapport_du_rang4"]))
winner5=loto[,"nombre_de_gagnant_au_rang5"]
gain5=as.numeric(as.character(loto[,"rapport_du_rang5"]))
which1=(sortconsec[,1]==1)
which2=(sortconsec[,2]==1)
which3=(sortconsec[,3]==1)
which4=(sortconsec[,4]==1)
which5=(sortconsec[,5]==1)```

There several ways to defining "winning at the lottery" (2 out of 6, 3 out of 6, 4 out of 6, etc) and to define "having consecutive numbers" (it can be 2 out of 6, 3 out of 6, etc). For instance,

It is also possible to compare the number of winners obtained for medium winners (3 out of 6, so called vainqueur de rang 4) when there were 2 consecutive numbers

```> t.test(winner4[which1==TRUE],winner4[which1==FALSE])

Welch Two Sample t-test

data:  winner4[which1 == TRUE] and winner4[which1 == FALSE]
t = -3.2132, df = 4792.491, p-value = 0.001321
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-6864.430 -1662.123
sample estimates:
mean of x mean of y
33887.82  38151.10 ```

With a simple mean comparison test, we have that there is a significant difference between average number of winners when there were at least 2 consecutive numbers out of 6 balls drawn. And actually, the average number of winners was lower when there are consecutive numbers. And if we look at the average gain, we have also a significant difference

```> t.test(gain4[which1==TRUE],gain4[which1==FALSE])

Welch Two Sample t-test

data:  gain4[which1 == TRUE] and gain4[which1 == FALSE]
t = 5.8926, df = 3675.361, p-value = 4.143e-09
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
11.06189 22.09337
sample estimates:
mean of x mean of y
173.9788  157.4012 ```

Here we see that if we play consecutive numbers, on average, the gain is larger. Perhaps it would be better to look at that on graphs.

```which0=which1
WIN=c(mean(winner5[which0==TRUE]),mean(winner5[which0==FALSE]),
mean(winner4[which0==TRUE]),mean(winner4[which0==FALSE]),
mean(winner3[which0==TRUE]),mean(winner3[which0==FALSE]),
mean(winner2[which0==TRUE]),mean(winner2[which0==FALSE]),
mean(winner1[which0==TRUE]),mean(winner1[which0==FALSE]))
MWIN=matrix(WIN,2,5)

plot(1:5,MWIN[1,],type="b",col="red",log="y",
ylim=c(1,1000000),xlab="two consecutive numbers",ylab="number of winners (log scale)")
lines(1:5,MWIN[2,],type="b",col="blue",pch=4)```

If we focus on the case where "having consecutive numbers" means two consecutive numbers, we have below the number of winners, with first rank (6 out of 6), then second rank (5 out of 6), etc,

Note that the y-axis is on a log scale, and that draws with consecutive balls are in red, and no consecutive balls are in blue. If we focus on average gains, curves are in opposite order,

But if we consider the case of three consecutive balls, we have, for the number of winners,

or for average gains

Here it starts to get slightly different: there are more "big winners" when there are at least three consecutive numbers. And with four consecutive numbers, it is clearly the opposite

Here we see that there are much more winners with four consecutive numbers (actually, it might be a triplet and a pair). So I have to confess that I am not convinced by the conclusion: actually a lot a people pick consecutive numbers... Actually, if we look at draws where there were the more winners, we can clearly see that a lot of players like consecutive numbers (perhaps not has much as playing birthdays, since most numbers are lower than 31),

```loto[loto\$"nombre_de_gagnant_au_rang1">50,
c("combinaison_gagnante_en_ordre_croissant",
"nombre_de_gagnant_au_rang1","date_de_tirage")]
combinaison_gagnante gagnant_au_rang1
3189      2-4-13-16-28-31              103
3475       1-5-9-10-12-25               64
4018       4-5-7-14-15-17               63
4396    26-27-28-35-36-37               64
4477     7-11-15-27-33-44               53
4546      2-9-12-14-19-24               60
4685      2-8-10-12-14-16               96
date_de_tirage
3189       19930626
3475       19920212
4018       19880504
4396       19840919
4477       19830914
4546       19820519
4685       19790919```

On September 1979, there were 5 even consecutive numbers (ok, it can not be considered strictly as consecutive numbers) and 96 winners with 6 numbers out of 6 ! And if we look at the others, even if they are not strictly consecutive, there is a lot of regularity. So I believe that picking consecutive numbers might not be a great strategy if you want to win a lot of money !

Wednesday, July 13 2011

## De la créativité des gangsters

Pendant mon séjour récent en Nouvelle Angleterre, j'ai survolé le livre de Leonard Mlodinow, the drunkard's walk. Et au hasard de mes lectures, je suis tombé sur la petite histoire suivante

Autrement dit, en utilisant les cinq derniers chiffres d'une quantité économique comme le us treasury balance, on aurait un générateur de nombre aléatoire...
Par contre la suite est un peu plus surprenante,

La loi de Benford (que j'avais pu évoquer ici ou ) parle des premiers chiffres, mais cette fois on parle des last five digits. Donc visiblement l'évoquation n'est pas pertinente ici. Mais qui sait ? Ca reste malgré tout une histoire intéressante. Considérons - histoire de tester cette légende - les deux sources de données suivantes, http://www.treasurydirect.gov/ et http://www.economagic.com/. Sur ce dernier site, un petit travail de mise en forme des données est nécessaire.
`b1=read.table("http://freakonometrics.free.fr/debtus1.txt",   header=TRUE,sep="\t")b2=read.table("http://freakonometrics.free.fr/debtus2.txt",   header=TRUE)X1=as.character(b1\$Dollar.Amount)n1=nchar(X1)Y1=substr(X1,n1-8,n1-3)X1=as.numeric(substr(Y1,1,2))*1000+as.numeric(substr(Y1,4,6))x=X1/100000X2=b2\$DEBTY2=trunc(as.numeric(X2))X2=as.character(Y2)n2=nchar(X2)Y2=substr(X2,n2-4,n2)y=as.numeric(Y2)/100000y=y[y<1]`

Pour rappel, un générateur aléatoire (standard) vérifie deux propriétés importantes

• les nombres doivent être tirés suivant une loi uniforme sur [0,1], i.e. ici, si on divise les nombres à 5 chiffres par 10000,
• les tirages doivent être indépendants entre eux.
La première propriété semble assez naturelle, et correspond à l'histoire racontée dans un commentaire posté ici (expliquant comment un casino avait été au bord de la faillite car une roulette faisait sortir certains chiffres trop souvent, et j'essayais de comprendre comment utiliser l'information qu'un chiffre sort plus souvent). La seconde est probablement encore plus importante.
• Visualisation de la distribution
La première idée est de visualiser la densité de nos séries de chiffres. Pour éviter les problèmes de bord (et comme c'est juste en introduction) on va utiliser un histogramme, et pas une estimation à noyau.
`hist(x,col="red")hist(y,col="blue")`
On obtient pour la première série la courbe rouge, et pour la seconde la courbe bleue,

On note qu'a priori, pour la première série, l'hypothèse d'uniformité n'est peut être pas la plus réaliste...
• Test de Kolmogorov-Smirnov
On peut aussi mettre en œuvre un test de Kolmogorov-Smirnov afin de tester si la loi uniforme est adaptée.:
`> ks.test(x,"punif") 	One-sample Kolmogorov-Smirnov test data:  x D = 0.1047, p-value = 0.01645alternative hypothesis: two-sided  > ks.test(y,"punif") 	One-sample Kolmogorov-Smirnov test data:  y D = 0.0456, p-value = 0.3581alternative hypothesis: two-sided `
On retrouve ici confirmée l'intuition précédante: la loi uniforme est pertinente pour la seconde série, pas la première.
• Les autocorrélations de la série
Travaillons uniquement sur la seconde série. On peut étudier l'autocorrélation de notre série de nombres, ou peut-être un peu plus malin, sur les quantiles gaussiens associés (les autocorrélations étant intéressantes pour les séries gaussiennes),:
`plot(acf(y))plot(acf(qnorm(y)))`
ie. pour la série brute

et pour la série normalisée,

Bref, on pourrait être tenté de valider l'hypothèse d'indépendance entre les tirages.
• Run test (de Wald–Wolfowitz)
L'idée est de comparer une série de chiffres à la médiane, s'ils sont plus grands, on note + (ou A) et sinon - (ou B). On crée alors une série du genre "+++−−++−−++++++−−−" et on compte les séries de + et les séries de -, les runs,
```library(lawstat)runs.test(y,plot=TRUE)

Runs Test - Two sided

data:  y
Standardized Runs Statistic = -0.2462, p-value = 0.8055```

Bref, la légende me semble à prendre avec des pincettes (car fonction de la source considérée), même si l'idée est intéressante (si l'on met de côté les aspects d'aléa moral). Et l'analyse sur la loi de Benford ne semble pas valide: les derniers chiffres sur les grands nombres ne se comportent pas du tout comme les premiers.

Friday, December 17 2010

## Non, les probabilités ne sont pas nos amies...

... et non, nous ne vivons pas le monde merveilleux des Bisounours: la vie n'est pas toujours juste. Les probabilités sont un objet assez curieux à manipuler, pas forcément très intuitif. J'avais évoqué ici, et enfin la notion d'espace probabilisé ... En fait, généralement, les calculs de probabilités sont assez troublants et paradoxaux. On va considérer ici des gens devant prendre des situations en présence d'incertitude, et plus particulièrement deux situations: la difficulté de prendre en compte de l'information (et donc de calculer des probabilités conditionnelles) et le cas où les gens doivent agir en prenant des décisions de manière aléatoire (et donc fixent en même les probabilités).

• Les probabilités conditionnelles

Avant hier, j'avais mis en ligne (ici) un billet sur le paradoxe dit de Monty Hall, où le candidat à un jeu télévisé doit choisir une porte parmi trois (afin de gagner une voiture), une fois son choix effectué, le présentateur ouvre une des deux autres portes (une qui n'a pas de lot derrière) et demande au candidat s'il souhaite revenir sur son choix. Et la réponse est qu'effectivement, le candidat a intérêt à ouvrir la porte qui est restée ouverte (et que l'animateur n'a pas choisie). L'idée du raisonnement est de formaliser proprement cette information à l'aide de lois conditionnelles.

Cette importance de l'information se retrouve par exemple dans l'histoire des cocus de Bagdad ou des yeux bleus (ici). Sinon, sur les paradoxes de probabilité, je peux aussi renvoyer vers des histoires de des lancers de pièces (ici), mais aussi sur le paradoxe des enveloppes () qui pourrait faire penser à l'histoire de Monty Hall. Mais mon préféré sur les lois conditionnelles reste le paradoxe de Simpson (ici, repris sous une forme plus géométrique )...

• Les probabilités comme paramètres de décision
Dans tous ces exemples, les probabilités sont connues, et il s'agit simplement d'une mise à jour afin de prendre en compte une nouvelle information. Mais il existe des cas encore plus vicieux, où les probabilités ne peuvent pas être considérées comme exogènes, ou objectives, car elles sont déterminées par les joueurs. Considérons une jeu simple, un peu dans l'esprit du pierre-feuille-ciseau (que j'avais déjà évoqué ici, en rappelant qu'il était fondamental de prendre des décisions de manière aléatoire). Dans le jeu d'aujourd'hui, on a toujours deux joueurs, disons moi (P) et ma fille (F). On cache nos mains dans le dos, et on les sort au même moment, en montant un ou deux doigts,
• si le nombre de doigts sortis est identique (paires (1,1) ou (2,2)), alors je gagne le nombre de doigts sortis, i.e. 2 ou 4
• si le nombre de doigts sortis est différent (paires (1,2) ou (2,1)), alors ma fille gagne le nombre de doigts sortis, i.e. 3
La matrice des paiements est alors la suivante

 P\F 1 2 1 2 -3 2 -3 4

Si on suppose que l'on a une chance sur 4 pour chacune des paires, le jeu est alors a somme nulle, puisque

Sauf que les probabilités sont fixées par les joueurs, en l'occurrence moi et ma fille. Supposons que je sorte 1 doigt avec probabilité , et 2 avec probabilité (avec ). Ma fille sort 1 doigt avec probabilité , et 2 avec probabilité (avec ). La valeur espérée du jeu (de mon point de vue) est

Ma fille cherche qui va minimiser cette valeur, alors que moi, je cherche qui la maximise. Bref, un problème classique de problème minmax de théorie des jeu. L'équilibre (s'il existe) sera et solution de

A p donné, ma fille cherche à trouver

soit

soit encore

en notant que , on peut réécrire . Autrement dit, si , ma fille a toujours intérêt à sortir un doigt (car ) et ce n'est pas un équilibre... pareillement, si , ma fille a toujours intérêt à sortir deux doigts, et là encore ce n'est pas un équilibre. Autrement dit, le seul équilibre possible est obtenu lorsque , i.e l'équilibre correspond à . De manière dual, pour ma part, je cherchais à trouver

i.e.

soit

qui conduit à un équilibre seulement si . Autrement dit, afin que le jeu soit à l'équilibre, il faut que nous jouions la stratégie suivante: sortir un doigt avec probabilité (et pas ) et deux avec probabilité . Et la valeur espérée du jeu est alors

Moralité, ma fille tu as intérêt à travailler tes maths si tu veux battre un jour ton père ce n'est pas vraiment un jeu "juste"... à condition bien sûr que moi et ma fille soyons ensuite capable de mettre en œuvre cette stratégie... Car si je regarde avec mes yeux de statisticiens, je pense qu'il est assez dur de faire la différence entre et ... La preuve sur cette petite simulation,
`> N=100; pv=rep(NA,100000)> for(s in 1:100000){+ X=sample(1:2,size=N,prob=c(5/12,7/12),replace=TRUE)+ P=prop.test(table(X)[1],n=N,p=1/2)+ pv[s]=P\$p.value}> mean(pv>0.05)[1] 0.66776`
autrement dit, en supposant que ma fille et moi arrivions a jouer 100 fois à ce jeu, dans 2/3 des cas je n'arrive pas à faire la différence entre une probabilité valant et une probabilité valant .