# Monte Carlo techniques to create counterfactuals

In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample $\{y_1,\cdots,y_n\}$, and more generally, in econometrics $\{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}$. But let’s get back to statistics (without covariates) to illustrate. We assume that observations $y_i$ are realizations of an underlying random variable $Y_i$. We assume that $Y_i$ are i.id. random variables, with (unkown) distribution $F_{\theta}$. Consider here some estimator $\widehat{\theta}$ – which is just a function of our sample $\widehat{\theta}=h(y_1,\cdots,y_n)$. So $\widehat{\theta}$ is a real-valued number like . Then, in mathematical statistics, in order to derive properties of the estimator $\widehat{\theta}$, like a confidence interval, we must define $\widehat{\theta}=h(Y_1,\cdots,Y_n)$, so that now, $\widehat{\theta}$ is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, $\widehat{\theta}$ is .

There are two strategies here. In classical statistics, we use probability theorem, to derive properties of $\widehat{\theta}$ (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, $\{y_1,\cdots,y_n\}$, and that’s a pity. But maybe we can create another one $\{y_1^{(1)},\cdots,y_n^{(1)}\}$, as realizations of $F_{\theta}$, and another one $\{y_1^{(2)},\cdots,y_n^{(2)}\}$, anoter one $\{y_1^{(3)},\cdots,y_n^{(3)}\}$, etc. From those counterfactuals, we can now get a collection of estimators, $\widehat{\theta}^{(1)}$,$\widehat{\theta}^{(2)}$, $\widehat{\theta}^{(3)}$, etc. Instead of using mathematical tricks to calculate $\mathbb{E}(\widehat{\theta})$, compute $$\frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}$$That’s what we’ve seen last friday.

I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students,

Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] X=Davis$height[Davis$sex=="M"]

We can visualize its distribution (density and cumulative distribution)

u=seq(155,205,by=.5) par(mfrow=c(1,2)) hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) lines(u,dnorm(u,178,6.5),col="black") Xs=sort(X) n=length(X) p=(1:n)/(n+1) plot(Xs,p,type="s",col="blue") lines(u,pnorm(u,178,6.5),col="black")

Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a $\mathcal{N}(178,6.5^2)$. In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals

hist(X,col=rgb(0,0,1,.3)) lines(density(X),col="blue",lwd=2) Y=rnorm(n,178,6.5) hist(Y,col=rgb(1,0,0,.3)) lines(density(Y),col="red",lwd=2) Ys=sort(Y) plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205)) polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA) lines(Xs,p,type="s",col="blue",lwd=2) lines(Ys,p,type="s",col="red",lwd=2)

We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution $\widehat{F}$ obtained from $\{y_1,\cdots,y_n\}$ (the blue curve), and some conterfactual, $\widehat{F}^{(s)}$ obtained from $\{y_1^{(s)},\cdots,y_n^{(s)}\}$ generated from $F_{\theta_0}$ – where $\theta_0$ is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance.

d=rep(NA,1e5) for(s in 1:1e5){ d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic } ds=density(d) plot(ds,xlab="",ylab="") dks=ks.test(X,"pnorm",178,6.5)$statistic id=which(ds$x&gt;dks) polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA) abline(v=dks,col="red")

If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic $\widehat{d}^{(1)}$, $\widehat{d}^{(2)}$, etc, and compare it with the one observe on our sample $\widehat{d}$. The proportion of samples where the test-statistics exceeded the one observed

mean(d&gt;dks) [1] 0.78248

is the computational version of the $p$-value

ks.test(X,"pnorm",178,6.5)   One-sample Kolmogorov-Smirnov test   data: X D = 0.068182, p-value = 0.8079 alternative hypothesis: two-sided

I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations).

This Monday, I will be giving the first part of the (crash) graduate course on advanced tools for econometrics. It will take place in Rennes, IMAPP room, and I have been told that there will be a visio with Nantes and Angers. Slides for the morning are online, as well as slides for the afternoon.

In the morning, we will talk about smoothing techniques, and in the afternoon, it will be on simulations and bootstrap techniques.

On Thursday, March 9nd, I will give the second lecture of the PhD course on advanced tools for econometrics, on simulation techniques (and bootstrap). Slides are available online.

The first part is this Thurdays, on Nonlinearities in Econometric models.

# Statistical Tests: Asymptotic, Exact, ou based on Simulations?

This morning, in our mathematical statistics course, we’ve been discussing the ‘proportion test‘, i.e. given a sample of Bernoulli trials$\boldsymbol{x}=\{x_1,\cdots,x_n\}$, with $X_i\sim\mathcal{B}(p)$, we want to test

$H_0:p=p_0$against $H_1:p\neq p_0$

A natural test (which can be related to the maximum likelihood ratio test) is  based on the statistic

$T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$

The test function is here

$\psi(\boldsymbol{x})=\boldsymbol{1}(T(\boldsymbol{x})\notin [c_{1,\alpha},c_{2,\alpha}])$To get the bounds of the acceptance region, we need the distribution of $T(\boldsymbol{X})$, under $H_0$. Consider here a numerical application

n=20
p=.5
set.seed(1)
echantillon=sample(0:1,size=n,
prob=c(1-p,p),
replace=TRUE)
• the asymptotic distribution

The first (and standard idea) is to use the central limit theorem, since

$\sqrt n\frac{\hat{p} - p}{\sqrt{p (1-p)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$

So, under $H_0$,

$\sqrt n\frac{\hat{p} - p_0}{\sqrt{p_0 (1-p_0)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$

Then $c_{1,\alpha}=\Phi^{-1}(\alpha/2)$ while $c_{2,\alpha}=\Phi^{-1}(1-\alpha/2)$. The acceptance region is then between the two red lines, below,

T=sqrt(n)*(mean(echantillon)-.5)/
sqrt(mean(echantillon)*
(1-mean(echantillon)))
u=seq(-3,3,by=.01)
v=dnorm(u)
plot(u,v,type="l",lwd=2)
abline(v=qnorm(.025),col="red")
abline(v=qnorm(.975),col="red")
abline(v=T,col="blue")

• the exact distribution

Here we use the fact that

$\sum_{i=1}^n X_i \sim \mathcal{B}(n,p)$

Using transformation of the ‘density’, we can (at least numerically) compute the (exact) distribution of

$T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$

u=seq(-3,3,by=.01)
v=sqrt(.5*(1-.5))*n*dbinom(round(
(sqrt(.5*(1-.5))*u/sqrt(n)+.5)*n),
size=n,prob=.5)/sqrt(n)

Here I used a round value, it guess it would be better with a floor function, but here the graph looks symmetric (which is something I like)

abline(v=sqrt(n)*(qbinom(.025,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
abline(v=sqrt(n)*(qbinom(.975,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
lines(u,v,type="s")

• distribution based on Monte Carlo simulations

Probably more interesting, here we do not use the fact that we might know the distribution of the mean. We just generate random samples, under $H_0$, and then compute $T(\boldsymbol{X})$,

T=rep(NA,1000)
for(i in 1:1000){
x=sample(0:1,size=n,
prob=c(1-.5,.5),
replace=TRUE)
m=mean(x)
T[i]=(m-.5)/sqrt(m*(1-m))*sqrt(n)}
lines(density(T),lwd=2)
abline(v=quantile(T,.025),col="red")
abline(v=quantile(T,.975),col="red")

# Bristish Statisticians and American Gangsters

A few months ago, I did publish a post (in French) following my reading of Leonard Mlodinow’s the Drunkard’s Walk. More precisely, I mentioned a paragraph that I found extremely informative

But it looks like those gangsters were not only stealing money. They were also stealing ideas, here from a British statistician, manely Leonard Henry Caleb Tippett. Leonard Tippett is famous in Extreme Value Theory for his theorem (the so-called Fisher-Tippett theorem, which gives the possible limiting distributions for a normalized version of the maximum from an i.i.d. sequence, see old posts). According to Martin Gardner, Leonard Tippett suggested to use middle numbers (not the last ones) of larger ones to generate (pseudo) random sequences, or more precisely, in 1927, “published a table of 41,600 random numbers, obtained by taking the middle digits of the area of parishes in England

I could not get a copy of the book Random Sampling Numbers by Leonard Tippett (I could only find reviews, e.g. Nair (1938)) but I do believe that this technique should work to generate sequences that do look like sequences of random numbers. Note that several techniques were mentioned in previous posts (in French) published a few years ago.

Now, I should also take some time to apologize because, sometimes, I am the one playing the gangster: I do steal a lot of illustrations on the internet. And I would like to apologize to the authors. On my previous blog, I did try – once – to add a short line at the end of a post, explaining where the illustration was coming from (trying to give credit to the illustrator). Less than 10 days after adding this short line, I received an email from a ‘publisher’, telling me that there were rights attached to the picture, and that I had 24 hours to remove it (if not, their lawyers will see what to do). Of course, I did remove the picture, and the mention. Now, I use pictures, and no mention. And I feel guilty. So I wanted to apologize for stealing others’ work. I am still discussing to hire an illustrator, to illustrate my blog. Work in progress….

# Ruin probability and infinite time

A couple of weeks ago, I had a discussion with a practitioner, working in some financial company, about ruin, and infinite time. And it reminded me a weird result. Well, not a weird result, but a result I found disturbing, at first, when I was a student (that I rediscovered with the eyes of someone dealing with computational issues, seeing here a difficult theoretical question). Consider a simple ruin problem. A player has wealth . Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time , his wealth is where  is associated to the th coin:  is equal to 1 with probability (tails), and -1 with probability  (heads). It is also possible to write

where  can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume  to be given. Let denote the number of heads and  the number of tails. Then , while . Let  denote the number of paths to go from point A (wealth  at time ) to point B (wealth  at time ). Note that this is a Markovian problem, that can be modeled using Markov chains

But here, we will focus on combinatorial results. Hence,

In order to derive probabilities to reach , let  denote the number of paths going from  to . And let denote the number of paths going from  to  that do reach  at some point between  and . Using a simple reflexion property, then if  and  are positive,

Based on those reflexions, two results can be derived (focusing on probability, instead of counting paths). First, we can obtain that

(given that n and x have the same parity). The second result we can obtain is that

Based on those two expressions, if  denotes the first time  become null, given ,

then

This can be computed easily,

> x=10
> p=.55
> ProbN=function(n){
+ pb=0
+ if(abs(n-x) %% 2 == 0)
+ pb=x/n*choose(n,(n+x)/2)*(1-p)^((n+x)/2)*(p)^((n-x)/2)
+ return(pb)}
> plot(Vectorize(ProbN)(1:1000),type="s")

That looks nice… But if we look closer, we can wonder what

would be ? Since we have the distribution of a probabilty measure, we might expect one. But here

> sum(Vectorize(ProbN)(1:1000))
[1] 0.134385

And this is not due to calculation mistakes that we do not get 1 here. Actually, we should write

which might be interpreted as the probability of ruin, starting from , that we denote  from now on. The term on the left can be approximated using monte-carlo simulations

> p=.55
> x=10
> m=1000
> simul=10000
> S=sample(c(-1,1),size=m*simul,replace=TRUE,prob=c(1-p,p))
> MS=matrix(S,simul,m)
> for(k in 2:m) MS[,k]=MS[,k]+MS[,k-1]
> T0=function(vm) which(vm<=(-x))[1]
> MTmin=apply(MS,1,T0)
> mean(is.na(MTmin)==FALSE)
[1] 0.1328

To check the validity of the relationship above, a simple (theoretical) recursive formula can be derived for the term on the right (ruin probability), namely

with a boundary conditions , and . Then is comes that

Note that it might be tricky to check using monte carlo simulation… since we cannot have an infinite number of runs. And we’re dealing precisely with things that do occur when time is infinite. Actually, we can still check convergence, considering an upper limit  for the number of runs, and then letting  go to infinity. Note that an explicit formula can then be derived (using additional border condition )

Using the following code, it is possible to calculate ruin probability, in order to estimate .

> MSmin=apply(MS,1,min)
> mean(MSmin<=(-x))
[1] 0.1328
> (((1-p)/p)^x-((1-p)/p)^m)/(1-((1-p)/p)^m)
[1] 0.1344306

The following graph shows the evolution of ruin probability as a function of initial wealth (with monte carlo simulation, with a fixed horizon – including a confidence interval – versus the analytical expression)

Hence, with stopping times, one should remember that

and that those two terms can be approximated simply using simulations or standard approximations.

# L’espérance est un opérateur linéaire, so what ?

Oui, “l’espérance est un opérateur linéaire”. On n’arrête pas d’insister sur cette propriété dans la plupart des cours de probabilité. En fait, je pense que cette propriété a deux implications importantes. La première est que l’on interprète souvent cette phrase en disant que si  intégrable, alors pour tout  et . La conséquence de cette propriété (qui est juste) est que  dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).

Mais cette linéarité dit davantage. Elle dit que pour toutes variables aléatoires et  (définies sur le même espace), , pour tout  et . Et ça c’est fort. En particulier, on ne dit pas qu’il faut que  et  sont des variables indépendantes. Car cette propriété est valide quelle que soit la dépendance qui pourrait exister entre  et .

Pourtant @Arnaud (qui passe souvent à l’occasion sur le blog) m’a laissé un couple de questions par courriel: “la dépendance modifie-t-elle la moyenne ?[…] sur des exemples simples en R, la moyenne est toujours perturbée”. L’idée était que les espérances étaient approchées numériquement par simulation. Or les méthodes de Monte Carlo reposent sur la loi des grands nombres et sur le théorème central limite. En particulier, si , et que l’on simule les couples , alors

autrement dit

i.e. la vitesse de convergence est affectée par la dépendance, pas la valeur vers laquelle on va tendre (car la moyenne tend vers l’espérance, qui est linéaire). En particulier, la convergence sera d’autant plus lente que  sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.

Par exemple, si on s’amuse à simuler 100 couples  Gaussiens, avec

on a les valeurs suivantes pour la moyenne de ce couple (obtenu en générant 10,000 échantillons de 100 paires), avec en abscisse la corrélation entre les deuxvariables

le trait rouge est la moyenne empirique des moyennes obtenues, et les régions sont délimitées par les quantiles à 5% (et 95%), 10% (et 90%) et 25% (et 75%). En moyenne, les simulations donnent la même chose, peu importe la dépendance. Mais plus la dépendance est forte, plus la variance sera grande, et plus la convergence sera lente (et les simulations auront d’autant plus tendance à s’éloigner – ou se disperser autour – de la vraie valeur). Mais cette réponse est partielle, car on a supposé que la variance était finie. Or cette hypothèse n’est pas nécessaire pour assurer la convergence de la moyenne vers l’espérance (et donc pour utiliser des méthodes de simulations). Que se passerait-il avec des lois de variance infinie ?

Pour autre chose que les vecteurs Gaussiens (par exemple des variables de variance infinie, comme des lois de Pareto) on peut utiliser le code suivant (on continue pour l’instant à utiliser une copule Gaussienne),

library(mnormt)
library(copula)
library(evir)
set.seed(1)
ns=100
R=seq(-.95,.95,by=.05)
M=BINF1=BSUP1=rep(NA,length(R))
for(i in 1:length(R)){
r=R[i]
norm.cop = normalCopula(r, dim = 2)
S=rep(NA,20000)
for(j in 1:20000){
U=rcopula(norm.cop,ns)
X=cbind(qgpd(U[,1], .7 ,1),qgpd(U[,2], .5 ,1))
S[j]=mean(X[,1]+X[,2])}
BINF1[i]=quantile(S,.05)
BSUP1[i]=quantile(S,.95)
M[i]=mean(S)
}
plot(R,M,col="white",ylim=range(c(BINF1,BSUP1)))
polygon(c(R,rev(R)),c(BINF1,rev(BSUP1)),
col="light yellow",border=NA)
lines(R,BINF1,lty=2)
lines(R,BSUP1,lty=2)
lines(R,M,lwd=2,col="red")
esperance=1+1/(1-.7)+1+1/(1-.5)
segments(min(R),esperance,max(R),esperance,col="blue")

Sur cette somme de deux lois de Pareto (comme dans le code ci-dessus), on obtient numériquement

Autrement dit, ici, l’effet dépendance semble avoir un impact plus faible sur l’approximation de la moyenne. On peut aussi essayer de changer la copule, par exemple mettre une copule de type Clayton (avec en abscisse cette fois le tau de Kendall),

(qui présente de la dépendance dans la queue inférieure) ou la copule duale (qui présente de la dépendance dans la queue supérieure)

Bref, la structure de la dépendance ne semble pas trop changer l’estimation de la moyenne (ou plutôt la vitesse de la convergence). On peut aussi se demander s’il y a un effet dimension. Par exemple, si au lieu de simuler des couples, on simule des vecteurs de plus grande taille. Par exemple, un vecteur échangeable en dimension 10, avec des marges normales centrées réduites,

norm.cop =normalCopula(r, dim = 10,dispstr="ex")

Moralité ? effectivement, la dépendance a un impact sur la vitesse de convergence quand on travaille sur des vecteurs aléatoires. Mais de manière assez surprenante, c’est surtout le cas pour la loi normale. Et étrangement, cet impact ne semble pas lié à de la dépendance extrême, ou à la dimension…

# A Million Random Digits: review of reviews

Recently on his blog (here), Robin mentioned an amazing book, called “A Million Random Digits” published by RAND corporation. The book was initially published in 1955, but RAND published a nice (and expensive) second edition.

A great thing is that on Amazon, there are several extremely interesting reviews of the book. E.g.

Didn’t like the ending, February 10, 2009  By Damien Katz

Even though I didn’t really see it coming, the ending was kind of anti-climatic. But overall the book held my attention and I really liked the “10034 56429 234088” part. It’s nice to know I’m not the only one who feels that way.

I found a typo, September 14, 2007  By fanfan

To whom do I write to report typographical errors? I noticed that the first “7” on the third line page 48 should be a “3”. The “7” that’s printed there now isn’t random. Other than that, this is really an excellent book.

Superb and original plot, April 21, 2007  By Herr Tarquin Biskuitfaß

This one has a very unpredictable plot, sublime character development in a style that stubbornly defies any sort of development in its rare and iconoclastic brilliance, and is told remarkably with numbers instead of letters. Take, for example, this passage on page 202, “98783 24838 39793 80954”. I’m speechless. The symmetry is reminiscent of the I Ching, and it approaches a rare spiritual niveau lacking in American literature. It not only reads well, but it looks great too. I have a tattoo of page 214 on my arm, and I’m hoping to get 202 on my belly to celebrate my next birthday. It is an injustice that Rand Corporation has not received the Nobel Prize for Literature, nor even a Pulitzer.

A serious reference work?, October 16, 2006  By BJ

For a supposedly serious reference work the omission of an index is a major impediment. I hope this will be corrected in the next edition.

Not Nearly A Million, September 3, 2006  By Liron

This book does not even come close to delivering on its promise of one million random digits. My expectations were high after reading the first sentence, which contained ten unique digits. However, the author seems to have exhasted his creativity in this initial burst, because the other 99.999% of the book is filler in which those same ten digits are shamelessly reused!  If you are looking for a larger offering of numerals in various bases, I highly recommend “Peter Rabbit’s ABC and 123”.

Wait for the audiobook version, October 19, 2006  By R. Rosini “Newtype”

While the printed version is good, I would have expected the publisher to have an audiobook version as well. A perfect companion for one’s Ipod.

Wait for it…, February 10, 2009  By Cranky Yankee

It started off slow, single digit slow in the beginning but I stuck with it. I eventually learned all about the different numbers, 1,2,3,4,5,6,7,8,9 and 0 and their different combinations.  The author introduced them all a bit too quickly for my taste. I would have been perfectly happy with just 1,2,3,4 and 5 for the first 20,000 digits, but then again, I’m not a famous random-number author, am I?  After a while, patterns emerged and the true nature of the multiverse was revealed to me, and the jokes were kinda funny. I don’t want to spoil anything but you will LOVE the twist ending!  Like 4352204 said to 64231234, “2242 6575 0013 2829!”

Ok, I have to admit I tried to check a few of them (that’s my freaky part). For instance the first one is a fake: the two first numbers – for instance – never show up together (consecutively),

> DIGIT=read.table("
+ http://freakonometrics.blog.free.fr/public/data/digits.txt")
> DIGIT=DIGIT[,2:11]
> k=1
> I=apply(DIGIT[,1:2]==c(10034,56429),1,sum)==2
> for(k in 2:9){
+ I=cbind(I,apply(DIGIT[,k+0:1]==c(10034,56429),1,sum)==2)
+ }
> I0=which(apply(I,1,sum)>0)
> DIGIT[I0,]
[1] V2  V3  V4  V5  V6  V7  V8  V9  V10 V11
<0 rows> (or 0-length row.names)

Nevertheless, I did have some fun reading those reviews. About the book, unfortunately I have to confess I stopped after 99998 appeared (the first time).

# Lottery, and martingales

I recently got a comment on a post I published one year ago, here, about the fact that in September 2009, on the 6th and the 10th, the 6 same numbers came out at the lottery, in Bulgaria (but  I do not understand the question: the author of the comment ask about the order the numbers came out…)
Xi’an published also a post on that topic, there, since last week, the same thing happened in Israel.
All that reminded me a discussion I had with a colleague about another post (here) where I mentioned that I found a strange distribution of numbers in the French lottery (the old one actually). For those who want to check, all historical events are here, in a zip file. My colleague was wondering if I found the martingale to win the lottery…

First, I do not like that term, since martingale is something different from a mathematical point of view… Second, let us look if it would have been possible to make some money… (free lunch ?)

> loto=read.table("D:\\loto.csv",dec=",",header=TRUE,sep=";")
> ntirage=nrow(loto)
> loto=loto[51:ntirage,]
> ntirage=nrow(loto)
>   N=as.matrix(loto[,c("boule_1","boule_2","boule_3","boule_4","boule_5","boule_6")])
> n=as.vector(N)
> length(n)
[1] 28848
> (TN=table(n))
n
1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
607 576 571 618 579 598 608 582 588 590 562 577 577 580 591 630 558 567 594 608
21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40
578 562 579 583 574 589 602 572 550 598 604 582 545 646 597 618 599 636 609 588
41  42  43  44  45  46  47  48  49
576 589 577 585 618 596 560 571 604

So, it might look nice, but we have to compare that distribution with the one we should have with “independent” draws. It is not possible to look at a discrete uniform distribution: the six numbers are not independent. Each day, the 49 balls are back in the urn, but within a day, we do not have independent draws (it is a sample without replacement of balls). Hence, with 4808 lottery draws, each number cannot be obtained more than 4808 times. So, let us use monte carlo techniques to  look at the theoretical distribution,

> M=matrix(NA,49,1000)
> for(s in 1:1000){
+ B=NA
+ for(i in 1:ntirage){B=c(B,sample(1:49,size=6,replace=FALSE))}
+ B=B[-1]
+ M[,s]=sort(table(B))
+ }
> q50=function(x){quantile(x,.5)}
> Q50=apply(M,1,q50)
> lines(1:49,Q50,col="red",lwd=2)
> q10=function(x){quantile(x,.1)}
> Q10=apply(M,1,q10)
> q90=function(x){quantile(x,.9)}
> Q90=apply(M,1,q90)
> polygon(c(1:49,49:1),c(Q10,rev(Q90)),col="light blue",border=NA)
> lines(1:49,Q10,col="red",lty=2)
> lines(1:49,Q90,col="red",lty=2)
> lines(1:49,Q50,col="red",lwd=2)
> points(1:49,sort(TN),pch=19,type="b")

Looking at the graph, it looks like some numbers appeared too frequently, especially the ones that did not appear frequently (bottom left). So, since I have removed the last 50 draws, let us see if we could have used that information, somehow…

> nb=names(sort(TN))
> loto=loto[1:50,]
> N=as.matrix(loto[,c("boule_1","boule_2","boule_3","boule_4","boule_5","boule_6")])
> n=as.vector(N)
> TN=table(n)
> TN[nb]
> barplot(TN[nb])

Unfortunately, numbers that came out too frequently over 4800 draws did not appear that frequently of the last 50. Playing top number might not have been a great strategy.

(numbers that came out frequently are on the right, while those we did not see much are on the left)… What about worst numbers: if I had decided to play the 6 that did not come out very frequently (we’ve seen earlier that they should have appeared even less, actually), would it have been interesting ? As we can see, our top 2 numbers were numbers that did not appear frequently earlier (29 and 47 appears respectively 10 and 11 times over 50 draws)….
Over 50 draws of 6 balls, the expected frequency of 6 given number is around 36.7,..

> S=rep(NA,10000)
> for(s in 1:10000){
+ B=NA
+ for(i in 1:50){B=c(B,sample(1:49,size=6,replace=FALSE))}
+ B=B[-1]
+ S[s]=sum(B%in%(1:6))
+ }
> mean(S)
[1] 36.7694

But here for the top 6, we have

> z=TN[nb]
> sum(rev(z)[1:6])
[1] 29

i.e. the top 6 appeared 29 times over 50 draw of 6 balls (which looks low) and for the worst 6, it is a bit higher,

>  sum(z[1:6])
[1] 38

If we look at the theoretical density of the frequency of 6 given number, we have

i.e. our worst 6 is a nice average (in green) while top 6 did not appear frequently this time (here in blue) ! So we could not have used that information….
Anyway, if some of you are interesting using statistics to get a free lunch, with the nouveau loto, I did not see any strange pattern (data can be downloaded here).

I am terribly sorry, but I cannot help anyone winning at the French Lottery….

# "Générer" du hasard, partie 5

La dernière méthode que  je présenterais est est celle des générateurs congruentiels linéaires. On retrouve l’idée de linéarité que l’on avait pas exemple dans le premier générateur proposé (ici), la notion de congruence est liée à des problèmes d’arithmétique et d’écriture de nombres dans les bases. On retiendra une fonction de récurrence de la forme suivante

$X_{n+1} = (a cdot X_n + c) mod m$

Pour que l’algorithme marche bien, on essaye de partir d’une valeur initiale correspondant à un nombre premier.
Prenons un premier exemple simple (inspiré par wikipedia). Si on prend des valeurs un peu au hasard, par exemple on obtient

• avec X0 = 10, la suite : 10, 10, 10, 10, 10, …
• avec X0 = 11, la suite : 11, 35, 123, 19, 235, 3, 91, 243, 203, 227, 59, 211, 171, 195, 27, 179, 139, 163, 251, 147, 107, 131, …
• avec X0 = 12, la suite : 12, 60, 236, 28, 204, 252, 172, 220, 140, 188, 108, 156, 76, 124, 44, 92, 12, 60, 236, 28, 204, 252, 172, …

Bref, on notera déjà qu’il est effectivement intéressant de partir d’un nombre premier, car on boucle alors beaucoup moins vite. Ah oui, et pour se ramener sur l’intervalle [0,1], on divise par m. On a alors intérêt à avoir m assez grand.
Prenons un second exemple simple (inspiré du générateur proposé par Robert Sedgewick). On prend , et on part de 0.

[1]         0        1 31415822 40519862 31536702 84962343
[7]  59428603 54128063 38284723  4802584 19281465 53057765
[13] 47900065 67928365 55662665 20022965 84329266 25717386
[19] 95163907 67972648 42464008 74270568 69856328 96165288
[25] 74221448 24728808  5671369 13328949 75902129 98182909

La série a l’allure suivante,
Et si on regarde l’histogramme et la fonction d’autocorrélation sur les 10 000 premières valeurs, on a gardé l’uniformité (ou tout du moins on ne doit pas en être trop loin), et on n’a plus d’autocorrélation.

Parmi les générateurs ultra-classiques, mais mauvais, on retiendra le RANDU, datant des années 60, implémenté alors sur les IBM, $X_{n+1} equiv (65539 times X_n) mod 2^{31}$ Mais cet algorithme présente de très gros défauts numériques. Je n’aborderais pas ici les problèmes de ce générateur (je l’avais fait dans le chapitre sur les simulations dans le tome 2 du livre avec Michel Denuit).
Bref, voilà ce qu’on peut raconter rapidement sur les méthodes de construction de nombres aléatoires, ou comme on dit généralement pseudo-aléatoires. Pour quelques réflexions sur l’utilisation de méthodes déterministes pour générer du hasard, je renvoie au très beau petit livre d’Ivar Ekeland intitulé Chaos. Et histoire d’insister un peu lourdement, ces billets (comme la plupart de ceux présent sur ce blog) n’ont pas grande prétention, si ce n’est présenter des concepts parfois compliqués “avec les mains” (et des dessins). Pour aller plus loin, on trouvera de très bons résumés de ces techniques sur le net, ici ou . Je laisse les plus curieux continuer. Sinon le chapitre 7 du Monte Carlo de Fishman propose une analyse poussée de ces méthodes. Et je voulais aussi noter que ce que j’ai proposé dans ces billets, c’est simplement de tirer des nombres aléatoires uniformément répartis sur [0,1]. Pour générer autre chose que de l’uniformité, je renvoie à la bible de Luc Devroye, non-uniform random variate generation, épuisé chez Springer, mais en ligne ici.

# "Générer" du hasard, partie 4

Bon l’idée de récurrence les dernières fois (ici et ) n’était pas forcément mauvaise, mais telle qu’elle était présentée, ça ne marchait pas très bien… Une autre méthode peut être celle proposée par John von Neumann, parfois appelée middle square. On prend un nombre à 4 chiffres, on l’élève au carré, et on prend les 4 chiffres du “milieu“.  Par exemple à partir de 1234, dont le carré est 1522756, puis on prend le carré de 2275, c’est à dire 5175625… etc. On prend alors comme nombre cette suite de nombres divisés par 10000 (pour revenir entre 0 et 1). Sur cette exemple on note que l’on peut assez rapidement boucler, et retomber sur la même suite (ceci reste le problème central de tous les algorithmes récursifs).

Bref, on n’y est toujours pas… Pourtant, de loi, la suite des nombres peut  sembler aléatoires (en oubliant le fait que l’on s’est restreint à 4 chiffres),

[1]  0.2756 0.5536 0.7296 0.1616 0.1456 0.9936 0.4096 0.7216
[9]  0.7065 0.4225 0.5062 0.3844 0.6336 0.4896 0.7081 0.4056
[17] 0.1136 0.0496 0.6016 0.2256 0.9536 0.5296 0.7616 0.3456
[25] 0.3936 0.2096 0.3216 0.2656 0.4336 0.8008 0.8064 0.8096
[33] 0.5216 0.6656 0.2336 0.6896 0.4816 0.3856 0.8736 0.7696

Mais si on regarde cet exemple de plus prêt, on note que le chiffre 6 apparaît très fréquemment à la fin, par exemple…. Même si cette méthode est souvent présentée (par exemple dans la page de wikipedia consacrée aux méthodes de simulations, ici), elle marche beaucoup moins bien que les méthodes que j’avais présentées auparavant (ici et ) et que l’on cite assez peu.

# "Générer" du hasard, partie 3

Pour poursuivre l’idée développée dans le billet précédant (ici), prenons une autre fonction de récurrence.

Toujours dans ces notes de cours (donné il y a un peu plus de 10 ans par Christian Gouriéroux), il y avait cette tent function que l’on retrouve sur le graphique de gauche. Si je fais tourner l’algorithme à partir de π-3, on obtient, avec une simple fonction R

[1]  0.14159265 0.28318531 0.56637061 0.86725877 0.26548246 0.53096491
[7]  0.93807017 0.12385966 0.24771932 0.49543864 0.99087728 0.01824545
[13] 0.03649090 0.07298179 0.14596358 0.29192717 0.58385434 0.83229132
[19] 0.33541736 0.67083471 0.65833057 0.68333886 0.63332228 0.73335543
[25] 0.53328913 0.93342173 0.13315654 0.26631308 0.53262615 0.93474770

Bref, ça marche bien là aussi (à condition là aussi d’avoir une écriture très longue de π). On peut là aussi montrer simplement que l’on génère une suite de nombre uniformément distribués sur [0,1], car cette loi est invariante par la fonction de récurrence (tout comme celle proposée par Whittle). Mais on n’a toujours pas d’indépendance….

# “Générer” du hasard, partie 2

Une première piste (qui n’est pas de moi mais qui a été proposé par Peter Whittle en 1963, et que j’ai retrouvé dans des notes prises lors d’un cours de Christian Gouriéroux) est une idée simple de récurrence. On prend un nombre u entre 0 et 1 (si possible irrationnel, sinon ça n’a pas grand intérêt), et on considère la suite. Par exemple si on prend comme valeur initiale π-3, on obtient, avec une simple fonction R

[1]  0.14159265 0.28318531 0.56637061 0.13274123 0.26548246 0.53096491
[7]  0.06192983 0.12385966 0.24771932 0.49543864 0.99087728 0.98175455
[13] 0.96350910 0.92701821 0.85403642 0.70807283 0.41614566 0.83229132
[19] 0.66458264 0.32916529 0.65833057 0.31666114 0.63332228 0.26664457

Bref, ça marche bien (à condition numériquement d’avoir une écriture avec un nombre infini de décimales).  Le plus simple est de passer en base 2 (on peut trouver π en base 2 sur le net, par exemple ici). Bref, on peut obtenir la série suivante,
Plus généralement sur le caractère aléatoire des décimales de π, je renvoie aux publications de David H. Bailey (ici).
Bref, cette méthode pourrait marcher car on génère effectivement des nombres uniformément sur l’intervalle [0,1]… sauf que les tirages ne sont pas du tout indépendants. On observe un beau processus autorégressif AR(1), avec une autocorrélation  en 2-h.
Bref, en l’état cette méthode n’est pas vraiment géniale pour générer du “hasard“… à suivre donc.

# “Générer” du hasard, partie 1

J’ai déjà mis plusieurs billets sur mon blog afin d’évoquer les méthodes de Monte Carlo (ici ou là), car de mon point de vue, ce sont des techniques fondamentales en économétrie par exemple, si on travaille sur des petits échantillons. Mais pour l’instant, j’avais toujours éludé le problème délicat du générateur… J’avais toujours admis que l’on avait à notre disposition un générateur de nombres aléatoires sur [0,1], qui génère à chaque appel des nombres indépendants, supposés uniformément distribués sur [0,1].

Je voulais prendre un peu de temps pour présenter quelques pistes de méthodes simples de méthodes permettant de générer des suites de nombre donnant l’illusion du hasard….
Une méthode plus ancienne (que j’ai connue, et oui) est celle des tables de nombres aléatoires, que l’on retrouvait par exemple dans les anciennes éditions du livre de Philippe Tassi ou de Gilbert Saporta. Je laisse les plus curieux traîner dans les vieux livres pour voir comment lire ce genre de table…. Une autre solution est celle que j’évoquais en introduction de mon (vieux) cours de séries temporelles à Dauphine (ici), à savoir l’idée de Slutsky d’utiliser les tirages du loto pour modéliser des prix (et donc inventer ainsi les processus MA, moving average), dans The Summation of Random Causes as the Source of Cyclical Processes. On peut trouver sur  le net les tirages du loto en France depuis des dizaines d’années (ici par exemple) et je suppose qu’en cherchant un peu, on peut trouver des choses similaires dans beaucoup d’autre pays.

Pour ma part, pour revenir sur une loi uniforme sur [0,1], je vais dans un premier temps diviser par 50.

Dans un second temps, on va concaténer deux boules consécutives (en base 50), puis se ramener sur [0,1].
Aussi, les boules 13 et 37 donnent comme nombre 0.2548.

On a ainsi un peu plus de possibilités sur les valeurs atteintes. Graphiquement, on obtient les résultats suivants,
Bref, on a à notre disposition des gens qui ont tiré des boules dans une urne. Alors rigoureusement, il faut faire attention car les tirages sont fait sans remise pour les 7 premières boules, puis on remet les boules pour le tirage suivant… Bref, les tirages ne sont pas rigoureusement indépendant les uns des autres (on ne peut pas avoir 3 fois consécutives le nombre “13”, alors qu’en théorie cela devrait être possible).
Bref, on est un peu limité avec ce genre d’outils, car ce qui m’intéresse c’est que mon ordinateur me génère autant de nombres aléatoires que je le souhaite.

## "sendo l'intento mio scrivere cosa utile a chi la intende…"

Search OpenEdition Search

You will be redirected to OpenEdition Search