# The odds of a cluster of airplane accidents

Recently, there have been a lot of airplane accidents.

• July, 17th 2014, Hrabove, Ukraine, Malaysia Airlines, Boeing 777, fatalities 298 (/298)
• July, 23rd 2014, Magong, Taiwan, TransAsia Airways, ATR 72-500, fatalities 47 (/58)
• July, 24th 2014, Aguelhok, Mali, Air Algerie, Mc Donnell Douglas MD-83, fatalities 116 (/116)

It is simple to find a lot of datasets about airplane crashes. For instance on http://ntsb.gov/aviationquery. The dataset is nice, with a lot of information,

> planes=read.table(
+ sep="|",header=TRUE)

for instance the exact location of the crashes,

> library(maps)
> map("world", interior = FALSE)
> points(planes$Longitude,planes$Latitude,
+ pch=19,cex=planes$Total.Fatal.Injuries/50, + col="red") # On Wigner’s law (and the semi-circle) There is something that I love about mathematics: sometimes, you discover – by chance – a law. It has always been there, it might have been well known by some people (specialized in some given field), but you did not know it. And then, you discover it, and you start wondering how comes you never heard about it before… I experienced that feeling this evening, while working on the syallbus for my course on copulas and extreme values. I discovered the so-called Wigner’s Semicircle Law (see e.g. Fan Zhang’s notes, or Fraydoun Rezakhanlou’s notes on that topic). Consider some $n\times n$ random matrice, with $n$ large (say 100) where elements are centered, such as a collection of random variable taking value $\{-1,1\}$ with equal probability. Then, eigenvalues can be visualized below n=100 M=matrix(sample(c(-1,1),size=n*n,replace=TRUE),n,n) E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))

Consider the symmetric matrix obtained from that matrix,

$\frac{M+M'}{2}$

and more precisely, let us look at its eigenvalues,

E=eigen(.5*(M+t(M)))$values Then the distribution of those eigenvalues is the so-called semi-circle distribution hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="", main="",border="white",xlim=c(-1.2,1.2),ylim=c(0,.65)) u=seq(-1,1,by=.01) v=sqrt(1-u^2)*2/pi lines(u,v,col=CL[6],lwd=2) Now, if we consider some $\mathcal{N}(0,1)$ distribution, instead of our binomial one, we got exactly the same M=matrix(rnorm(n*n),n,n) E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))
E=eigen(.5*(M+t(M)))$values hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="", main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65)) u=seq(-1,1,by=.01) v=sqrt(1-u^2)*2/pi lines(u,v,col=CL[6],lwd=2) Actually, it is a very general result, see the second chapter of an Introduction to Random Matrices by Greg Anderson, Alice Guionnet and Ofer Zeitouni, for instance. If entries of the random matrix are independent centred random variables, symmetric, such that higher moments exist, then this property is valid. That’s awesome, isn’t it? Because if the distribution has too heavy tails, then this property is no longer valid. For instance, if we consider a random matrix where entries have a Student distribution, we get something different… M=matrix(rt(n*n,df=2.1),n,n) M=M/sd(M) E=eigen(M)$values
plot(E,xlim=c(-11,11),ylim=c(-11,11))
E=eigen(.5*(M+t(M)))$values hist(E/sqrt(2*n),probability=TRUE,col=CL[4],xlab="",ylab="", main="",border="white",,xlim=c(-1.2,1.2),ylim=c(0,.65)) u=seq(-1,1,by=.01) v=sqrt(1-u^2)*2/pi lines(u,v,col=CL[6],lwd=2) (here, I do normalize by the standard deviation to get something comparable with the previous graph, where variables were centered, with unit variance) and if we consider a distribution with infinite variance, we get M=matrix(rt(n*n,df=1.75),n,n) we get I guess I will get back on that property in my course! # Halloween and candies (a ballot problem) This year, for Halloween, a post on candies (I promise, next year I will write another post on zombies). But I don’t want to focus on the kids problems (last year, we tried to minimize their walking distance to collect as much candies as possible, with part 1 and part 2), I want to discuss my own problems. Because usually, the kids wear their costumes, and they go in the streets, they knock on the doors, while I stay at home. So I’m the one, with a bag full of candies, waiting for kids to knock on our door, and then I give them some candies (if they wear a costume). Consider the following problem. Assume that we start with $r$ red candies, and $b$ black ones, with $r>b$. The thing is that no one like those black candies. What could be the probability that for the $n$ kids that will get candies after knocking at my door (with $n=r+b$ for convenience, but we will also consider the more general case where I have to many candies, $n\leq r+b$, later on), the probability to get a red candy is always larger than the probability to have a black candy ? This is somehow related to the popular ballot problem, proposed (and solved) by Whitworth in 1878, but he wrote it only in the fourth edition of Choice and Chance, in 1886 (this is what the legend told us). In 1887, Joseph Bertrand proposed a similar problem, and Désiré André introduced the reflection principle to solve it. The problem is simple : consider an election between two candidates, A (who receives $m$ votes) and B (who receives $n$ votes). A wins the election ($m>n$). If the ballots are cast one at a time, what is the probability that A will lead throughout the voting? For those who don’t remember the conclusion, the probability is here quite simple, $\mathbb{P}(\boldsymbol{A} > \boldsymbol{B})=\frac{m-n}{m+n}$ Observe that some geometry proofs were given, later on, by Aebly or Mirimanoff, both in 1923, as well as Howard Grossman in the 1950’s (see the discussion on http://academiclogbook.blogspot.ca/…). Actually, http://futilitycloset.com// produced the following geometric proof (with no clear reference), We start at O, where no votes have been cast. Each vote for A moves us one point east and each vote for B moves us one point north until we arrive at E, the final count, (mn). If A is to lead throughout the contest, then our path must steer consistently east of the diagonal line OD, which represents a tie score. Any path that starts by going north, through (0,1), must cut OD on its way to E. If any path does touch OD, let it be at C. The group of such paths can be paired off as p and q, reflections of each other in the line OD that meet at C and continue on a common track to E. This means that the total number of paths that touch OD is twice the number of paths p that start their journey to E by going north. Now, the first segment of any path might be up to m units east or up to units north, so the proportion of paths that start by going north is n/(m + n), and twice this number is 2n/(m + n). The complementary probability — the probability of a path not touching OD — is (m –n)/(m + n). But let’s try to solve our problem. Let $B_k$ and $R_k$ denote the number of black and red candies, respectively after the $k$th kid git his (or her) candy. Yes, one at a time. Here, $B_0=b$ and $R_0=r$. What we want is $\mathbb{P}(\boldsymbol{R}\geq \boldsymbol{B})=\mathbb{P}(\forall k\in\{0,1,\ldots,r+b\} : R_k \geq B_k)$ Using this formulation, we recognize the ballot problem. Almost. Actually, in the original ballot problem (see Bertrand (1887)), we have to compute the probability that one candidate remains strictly ahead the other one throughout the count. With a strict condition, we get the well-known probability (given previously) $\mathbb{P}(\boldsymbol{R}> \boldsymbol{B})= \frac{r-b}{r+b}$ Here, ties are allowed, and we can prove (easily) that $\mathbb{P}(\boldsymbol{R}\geq \boldsymbol{B})= \frac{r+1-b}{r+1}$ (again, there is some nice geometric interpenetration of that result). It is also possible to get numerically that value using the following function, which will generate a trajectory, and return some indicators (with or without ties) > red_black=function(sd){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v1=rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")) + v2=cumsum(rev(vectcandy)=="R")<= cumsum(rev(vectcandy)=="B") + return(list(evol=cbind(rev(cumsum(rev(vectcandy)=="R")), + rev(cumsum(rev(vectcandy)=="B")),v1),list=vectcandy,test=(sum(v1)==0), + ballot=(sum(v2)==0),when=min(which(v1==1))))} (here I compute the ballot-type index, where ties are not allowed, and the candy-type index). If we generate 100,000 scenarios, starting with 50 red and 25 black candies, we get > r=50 > b=25 > M=sapply(1:100000,red_black) ­­ > mean(unlist(M[3,])) [1] 0.50967 which can be compared with the theoretical value > (r+1-b)/(r+1) [1] 0.5098039 We can also get the distribution of the first time we have more black candies than red ones left (given that this event occur) > Z=unlist(M[5,]) > Z=Z[Z<Inf] > hist(Z,breaks=seq(0,80),probability=TRUE,col="light blue", + border=NA,xlab="",main="") There might be some analytically formula that can be derived, but I have to confess that I am becoming extremely lazy, Assume now that this year, kids do not show up at my door (for some reason). Assume that $n\leq r+b$ kids show up. We can see how the probability $\mathbb{P}(\boldsymbol{R}\geq \boldsymbol{B})=\mathbb{P}(\forall k\in\{0,1,\ldots,n\} : R_k \geq B_k)$ will change, with $n$, > r=50 > b=25 > impact_n = function(n){ + red_black=function(sd,nb=n){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v=(rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")))[1:nb] + return(list(list=vectcandy,test=(sum(v)==0),when=min(which(v==1))))} + M=sapply(1:10000,red_black) + return(mean(unlist(M[2,])))} Yes, not only I am too lazy to derive analytic formulas, I am so lazy that I do not try to optimize my code. Here, the evolution of the probability, as a function of $n$ is > V=Vectorize(impact_n)(25:75) > plot(25:75,V) Fun isn’t it? But now, I have to conclude my post, to work a little bit on my make-up : I have learnt so many thinks at the Montreal Zombie Walk a few days ago that kids willing to knock at my door will be scared to death. I guess I will keep all the candies for me this year ! # Generating a Markov chain vs. computing the transition matrix A couple of days ago, we had a quick chat on Karl Broman‘s blog, about snakes and ladders (see http://kbroman.wordpress.com/…) with Karl and Corey (see http://bayesianbiologist.com/….), and the use of Markov Chain. I do believe that this application is truly awesome: the example is understandable by anyone, and computations (almost any kind, from what we’ve tried) are easy to perform. At the same time, some French students asked me specific details regarding some old lectures notes on Markov chains, and on some introductory example I used as a possible motivation: the stepping stone algorithm. In the notes, I just mentioned the idea of this popular generic algorithm (introduced in Sawyer (1976)) and I use simulations to show – visually – how it works. Again, it was just to motivate the course which actually did focus on the theory of Markov Chains. But those student wanted more, like how did I get the transition matrix, for instance. And that is actually not a simple question, from a computational perspective. I mean, I can easily generate this Markov Chain, but writing explicitly the transition, that was another story. Which took me a bit longer. In a very specific case… But let us get back to the roots, and to the stepping stone algorithm. At least, one of them (the one I used in my notes) because it looks like there are several algorithm. We do consider a grid, say $h\times h$, with some colors inside, say $k$ possible colors. Each cell of the grid has a given color. Then, at some stage, we select randomly one cell in the grid, and it will take the color of one of its neighbor (some kind of absorption, or mutation). This is, more or less, what is also detailed in some lecture notes by James Propp (see also e Sato (1983) or Zähle et al. (2005) for more theoretical details about that Markov chain). This is extremely simple to generate (that’s what I did in my notes, with very big grids, and a lot of colors). But what if we want to write the transition matrix ? First of all, we need to define the state space. Basically, we do have $h^2$ cells, each of them has one color, chosen among $k$. Which gives us $k^{h^2}$ possible states…. And that can be large. I mean, if we consider the smallest possible grid (that might be interesting), say $h=3$, and only $2$ colors, then we talk about $2^9=512$possible states. That is large, not huge. But we should keep in mind that we have to compute a transition matrix, that would be a matrix with $512\times 512=262,144$ elements. More generally, we talk about writing down matrices with $k^{2\cdot h^2}$ elements. If we want black and white $4\times 4$ grids, that would mean a matrix with $2^{2\cdot 4^2}=2^{32}=4,294,967,296$ which mean 4 billion elements ! And if we consider an red-green-blue $3\times 3$ grid, we have to explicit a matrix with $387,420,489$ i.e almost 400 million elements. So, let’s face it: we can only work with $3\times 3$ bi-color grids. So let’s try… The good thing is that it can be related to work I’ve been doing recently on binomial recombining trees (binomial being related to bi-color). First of all, our grid will be describes as follows > h=3 > M=matrix(1:(h^2),h,h) > M [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 with two colors > color=c("red","blue") Then, we should look for neighbors, or derive an neighborhood matrix, > d=function(i,j) dist(rbind(c((i-1)%/%h,(i-1)%%h), + c((j-1)%/%h,(j-1)%%h))) > Neighb=matrix(Vectorize(d)(rep(1:(h^2),each=h^2), + rep(1:(h^2),h^2)),h^2,h^2) > trunc(Neighb*100)/100 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 0.00 1.00 2.00 1.00 1.41 2.23 2.00 2.23 2.82 [2,] 1.00 0.00 1.00 1.41 1.00 1.41 2.23 2.00 2.23 [3,] 2.00 1.00 0.00 2.23 1.41 1.00 2.82 2.23 2.00 [4,] 1.00 1.41 2.23 0.00 1.00 2.00 1.00 1.41 2.23 [5,] 1.41 1.00 1.41 1.00 0.00 1.00 1.41 1.00 1.41 [6,] 2.23 1.41 1.00 2.00 1.00 0.00 2.23 1.41 1.00 [7,] 2.00 2.23 2.82 1.00 1.41 2.23 0.00 1.00 2.00 [8,] 2.23 2.00 2.23 1.41 1.00 1.41 1.00 0.00 1.00 [9,] 2.82 2.23 2.00 2.23 1.41 1.00 2.00 1.00 0.00 > Neighb=(Neighb<2)&(Neighb>0) > Neighb [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE [2,] TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE [3,] FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE [4,] TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE FALSE [5,] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE [6,] FALSE TRUE TRUE FALSE TRUE FALSE FALSE TRUE TRUE [7,] FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE FALSE [8,] FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE TRUE [9,] FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE FALSE Now, let us explicit our 512 possible states. > n=h^2 > states=function(x){ + Base.b=rep(0,n) + ndigits=(floor(logb(x,base=length(color)))+1) + for(i in 1:ndigits){ + Base.b[n-i+1]=(x%%length(color)) + x=(x %/% length(color))} + return(Base.b)} > M=Vectorize(states)(1:(length(color)^n-1)) > liststates=data.frame(rbind(rep(0,h^2),t(M))) > head(liststates) X1 X2 X3 X4 X5 X6 X7 X8 X9 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 3 0 0 0 0 0 0 0 1 0 4 0 0 0 0 0 0 0 1 1 5 0 0 0 0 0 0 1 0 0 6 0 0 0 0 0 0 1 0 1 (for the first six, with 0/1 digits instead of colors). For instance, if we look at a specific one, it is possible to plot the grid, using > plotsteps=function(u){ + plot(0:h,0:h,col="white",xlab="",ylab="",axes=FALSE) + for(i in 0:(h^2-1)){ + x=i%/%h + y=i%%h + polygon(x+c(1,.1,.1,1),y+c(1,1,.1,.1), + col=color[as.numeric(u)[i+1] + 1]) + text(x+.45,y+.45,i) + }} Here, > plotsteps(liststates[100,]) Then, given one state, let us see what could happen next, • let us compute all connected states: all states where we can end up in if we change one cell • we have to check, for each connect state which cell did change • we should compute probabilities to reach those 9 states, based on the fact that each of the cell is chosen with the same probability, and the fact that probability to change the color is based on the colors around. • if some states cannot be reached (if a cell is surrounded by elements of the same color, so it cannot change its color), then, we should remove then from the list of reachable (possible) states. The code will be something like the following > listneighbour=function(i){ + start=liststates[i,] + difference2only=function(j) { + w=which(liststates[j,]!=liststates[i,]) + return((length(w)==1))} + possible=which( Vectorize(difference2only)(1:nrow(liststates))==TRUE ) + P=function(j){ + L=liststates[i,which(Neighb[which(liststates[j,]!=liststates[i,]),]==TRUE)] + T=table(as.numeric(L)) + T=T[as.character(0:(length(color)-1))] + T[is.na(T)]=0 + return(as.numeric(T)/sum(T)) + } + probability=Vectorize(P)(possible) + W=NULL + for(j in possible) W=c(W,which(liststates[j,]!=liststates[i,])) + I=1-liststates[i,W]+1 + vp=diag(probability[as.numeric(I),]) + vproba=0*vp + if(sum(vp)!=0) vproba=vp/sum(vp) + return(list( + color=liststates[i,W], + absorb=W, + possible=possible, + probability=probability, + prob=vproba)) + } For instance, if we start from state 100 (here, on the right) > listneighbour(100)$color
X3 X4 X8 X9 X7 X6 X5 X2 X1
100  1  1  1  1  0  0  0  0  0

$absorb [1] 3 4 8 9 7 6 5 2 1$possible
[1]  36  68  98  99 104 108 116 228 356

$probability [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1 0.8 0.6 0.6667 0.3333 0.4 0.5 0.6 0.6667 [2,] 0 0.2 0.4 0.3333 0.6667 0.6 0.5 0.4 0.3333$prob
[1] 0.17964072 0.14371257 0.10778443 0.11976048 0.11976048
[6] 0.10778443 0.08982036 0.07185629 0.05988024

Let us look more specificaly at the 99th state (which appears above as a state that could be reached from the 100th),

> liststates[99,]
X1 X2 X3 X4 X5 X6 X7 X8 X9
99  0  0  1  1  0  0  0  1  0

If we plot it (here on the right, again), we get

> plotsteps(liststates[99,])

Clearly, here, the cell in the upper corner (number 9) changed from blue to red. Now, about the probability… The probability to select cell 9 is 1/9, and given that cell 9 is chosen, the probability to go from blue to red is 2/3 (the cell is surrounded by 2 red cells, and 1 blue cell). The probability to remain blue is then 1/3. Those are the probabilities computed by our function (the table with two rows, one per color). In order to get a better understanding on the meaning of the last line, with some sort of probabilities), let us look at the following (simpler) example.

> liststates[2,]
X1 X2 X3 X4 X5 X6 X7 X8 X9
2  0  0  0  0  0  0  0  0  1

that can be visualized on the right (on the right). Here,

> listneighbour(2)
$color X9 X8 X7 X6 X5 X4 X3 X2 X1 2 1 0 0 0 0 0 0 0 0$absorb
[1] 9 8 7 6 5 4 3 2 1

$possible [1] 1 4 6 10 18 34 66 130 258$probability
[,1] [,2] [,3] [,4]  [,5] [,6] [,7] [,8] [,9]
[1,]    1  0.8    1  0.8 0.875    1    1    1    1
[2,]    0  0.2    0  0.2 0.125    0    0    0    0

$prob [1] 0.65573770 0.13114754 0.00000000 0.13114754 0.08196721 [6] 0.00000000 0.00000000 0.00000000 0.00000000 Things are pretty simple here • if we chose cells $\{1,2,3,4,7\}$, then nothing change, since all the neighbors have the same color. So if we want to focus on changes (or say run the algorithm until the first color change, then choosing those cells is a waste of time) • if we chose cells $\{5,6,8\}$, then it could be possible to change the color. And actually, $\{5\}$ is different from $\{6,8\}$ (since it does have much more neighbors) • if we chose cell $\{9\}$, then definitively, the color will change, since all neighbors have the other color here, Now, the probability to select cell $k$ given that there was a color change would be, if $k$ is in $\{9\}$ $\mathbb{P}(k)\propto \frac{3}{3}=1$ while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so $\mathbb{P}(k)\propto \frac{1}{5}$and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so $\mathbb{P}(k)\propto \frac{1}{8}$ And for the other, $\mathbb{P}(k)\propto 0$. So, it comes – since we assume that cells are drawn independently, and with the same probability, if $k$ is in $\{9\}$ $\mathbb{P}(k)= \frac{1 \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{40}{61}$ while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so $\mathbb{P}(k)= \frac{\frac{1}{5} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{8}{61}$ and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so $\mathbb{P}(k)= \frac{\frac{1}{8} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{5}{61}$ Which are exactly the probability computed above. The point is that we compute probabilities given that a color change will actually occur. The good point is that it should faster convergence to some limiting distribution. If any. What about our transition matrix ? Well, using a simply loop, we should get it easily > M=matrix(0,nrow(liststates),nrow(liststates)) + for(i in 1:nrow(liststates)){ + L=listneighbour(i) + if(sum(L$prob)!=0){
+ j=L$possible + M[i,j]=L$prob
+ }
+ if(sum(L$prob)==0){ + j=i + M[i,j]=1 + } + } One can check that this matrix satisfies some properties of transition matrices. For instance, the sum per row is one, > sum(apply(M,1,sum)!=1) [1] 0 Remember that this matrix is big, so I will not print if here. But trust me, it works (it might take a while on an old laptop, but anyone can do it). Now, if we want to visualize some paths of that chain, we can use the following algorithm. First, we need a starting point, that can be chosen randomly, > j=sample(1:nrow(liststates),size=1) or using a given colored grid, say > j=100 Then we plot it, > plotsteps(liststates[j,]) Now, the code within the loop is here > d=rep(0,nrow(liststates)) > d[j]=1 > d=d%*%M > j=sample(1:nrow(M),size=1,prob=d) > plotsteps(liststates[j,]) Here are some examples. And indeed, we end up either with all cells in blue, or all cells in red. Now, do we have to compute that transition matrix to produce those graph (and to generate that Markov chain) ? No. Of course not… At each step, I use a Dirac measure, and use the transition matrix just to get the probability to generate then the next state. Actually, one can write a faster and more intuitive code to generate the same chain… But I should probably keep that for another post… # Playing cards in Vegas? In a previous post, a few weeks ago, I mentioned that I will be in Las Vegas by the end of July. And I took the opportunity to write a post on roulette(s). Since some colleagues told me I should take some time to play poker there, I guess I have to understand how to play poker… so I went back to basics on cards, and shuffling techniques. Now, I have to confess that I have been surprised, while I was looking for mathematical models for shuffling, to find so many deterministic techniques (and results related to algebra, and cycles). On http://mathworld.wolfram.com/ for instance, one can find nice articles on so-called in-shuffle or out-shuffle techniques. There is also a great article, Golomb (1961), but mainly on algebraic properties of permutations by cutting and shuffling, as well as Diaconis, Kantor and Graham’s The Mathematics of Perfect Shuffle And if you look at Monge’s shuffle, you can find a deterministic recursive relationship. As a statistician (or applied probabilist), I should confess that I did not find answer to the question I wanted to ask : how long should we shuffle before getting cards randomly sorted in ours hands ? • Randomness (from a statistician perspective) First, I need to define (as properly as possible) a notion of “cards randomly sorted“. Consider a game with 32 cards. Why 32 ? Mathematicians will tell you that 32 is a great number, since it is a power of 2, so there might be interesting (algebraic) properties when shuffling. From a computational point of view, 32 is smaller than 52, so my random generations will run faster. This is basically why I used 32. 10 would have been better, but not realistic with cards. So, our 32 cards can be seen as a vector, or a list, of 32 items, say $(a_1,a_2,\cdots, a_{31},a_{32})$ In order to assess if my cards are randomly sorted, let us get back to number properties (real valued numbers). If there were 10 cards, the list can be seen as an element of the following set $(a_1,a_2,\cdots, a_{9},a_{10})\in\{0,1,\cdots,8,9\}^{10}$ (or to be more specific, a subset of that set, since numbers have to be different – it has to be a permutation – we cannot have duplicates, we’ll get back to that point in a few seconds). Let us see this list as a decimal number, with 10 digits. More precisely, $u=\frac{a_1}{10}+\frac{a_2}{10^2}+\cdots+\frac{a_9}{10^9}+\frac{a_{10}}{10^{10}}\in[0,1]$ Now, it is natural to say that cards are randomly sorted is this number is uniformly distributed on the unit interval, isn’t it ? (if we use the same shuffle many times, with the same starting point) Well, if we think about it twice, uniform on the unit interval is probably not the proper distribution, since (as mentioned above) all digits have to be different. For instance, the smallest number would be $0.0123456789$ and the largest $0.9876543210$ . But as we will see, it this uniform assumption might not be too strong, actually. And if we want to get back to our initial problem, with 32 cards, we simply have to use a decomposition in the 32-basis. $u=\frac{a_1}{32}+\frac{a_2}{32^2}+\cdots+\frac{a_{31}}{32^{31}}+\frac{a_{32}}{32^{32}}\in[0,1]$ So if we have an algorithm to shuffle cards, we just have to run it several times (with the same starting value) and see when $u$ starts to be uniformly distributed. We start with a Dirac distribution, we have some kind of transition matrix, we expect our limiting distribution to be uniform and we wonder when the limiting distribution is reached… And from a statistical point of view, that should not be that difficult to assess, since we do have several goodness of fit tests that can be used. Actually, it is possible to check if our technique passes the test of a uniform distribution, when digit are randomly generated (without replacement). The code to generate $u$ is > j = 32 > X3 = (0:(j-1))[sample(1:j)] > x3 = sum(j^(-(1:j))*X3) If we run it a few times, and check if the assumption of a uniform distribution is valid (on samples with, say, 500 observations), > P3=NULL > for(i in 1:10000){ + U3=NULL + for(s in 1:500){ + X3 =(0:(j-1))[sample(1:j)] + x3 =sum(j^(-(1:j))*X3) + U3 =c(U3,x3)} + P3 =c(P3,ks.test(U3,punif)$p.value)
+ }

in 95% of the scenarios, the $p$-value exceeds 5%

> mean(P3>.05)
[1] 0.9529

(which is something we should have under the null), More precisely, we can check that the $p$-value is uniformly distributed on the unit interval.

> hist(P3,freq=FALSE)

So assuming that our number is uniform on the unit interval might be a good notion for “cards are randomly sorted“.

What we need now is some shuffling algorithms. Or to be more specific, some feasible shuffling algorithm. I mean here that I just start playing with cards, so it should be some techniques that I should be able to perform, to understand how it works…. So you will have to wait a few weeks before I start talking about the riffle or dovetail shuffle (you know the kind of shuffle in which half of the deck is held in each hand, and then cards are released by the thumbs so that they fall to the table interleaved… like in the movies) !

• Top in at random shuffle, and related (simple) algorithm

My first algorithm is simple: the top-in at random shuffle. We start with the following ordering

    N=1:m

There are $m$ cards, and n denote the place where the card on top will go.

    n=sample(2:m,size=1)
if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])
if(n==m) N=c(N[2:n],N[1])

Then, we repeat that transfer of the card on top several times.

schuffle1=function(m,ns=10){
N=1:m
for(i in 1:ns)
{
n=sample(2:m,size=1)
if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])
if(n==m) N=c(N[2:n],N[1])
}
return(N)}

Now, it is also possible to consider a bottom-in at random shuffle. The idea is the same, the only difference it that you start from the card at the bottom of the deck. But that would be the same as the one before (in terms of time before reaching randomness)

    n=sample(1:(m-1),size=1)
if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])
if(n==1) N=c(N[m],N[1:(m-1)])

Why not mixing ? Randomly. Call it randomly mixed top-bottom in at random shuffle. You start either with the card on top, or at bottom (with identical probability), of the deck and then move the card somewhere,

     card=sample(c("top","bottom"),size=1)
if(card=="top"){
n=sample(2:m,size=1)
if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])
if(n==m) N=c(N[2:n],N[1])}
if(card=="bottom"){
n=sample(1:(m-1),size=1)
if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])
if(n==1) N=c(N[m],N[1:(m-1)])}

All those codes can be together (within the same function),

schuffle1=function(m,ns=10,which="top"){
N=1:m
if(which=="top"){
for(i in 1:ns)
{
n=sample(2:m,size=1)
if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])
if(n==m) N=c(N[2:n],N[1])
}}
if(which=="bottom"){
for(i in 1:ns)
{
n=sample(1:(m-1),size=1)
if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])
if(n==1) N=c(N[m],N[1:(m-1)])
}}
if(which=="mixed"){
for(i in 1:ns)
{card=sample(c("top","bottom"),size=1)
if(card=="top"){
n=sample(2:m,size=1)
if(n<m)  N=c(N[2:n],N[1],N[(n+1):m])
if(n==m) N=c(N[2:n],N[1])
}
if(card=="bottom"){
n=sample(1:(m-1),size=1)
if(n>1)  N=c(N[1:(n-1)],N[m],N[n:(m-1)])
if(n==1) N=c(N[m],N[1:(m-1)])
}
}}
return(N)}

But why do we take only one card ? It won’t be more complex to take 2. Or 3. Or more.

• Tops in at random shuffle, and related (mixed) algorithm

Yes, I used tops to say that we would take several cards on top of the deck. Say a random number of cards. And then, the strategy is the same, so the previous code is (slightly) adapted, as follows

     k=sample(1:(m-1),size=1)
n=sample((k+1):m,size=1); if(k==m-1) n=m
if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])
if(n==m) N=c(N[(k+1):n],N[1:k])

The idea is the following, here

As earlier, it is possible to take cards at the bottom of the deck, or, one more time, to use a mixed strategy. The codes would be

     card=sample(c("top","bottom"),size=1)
if(card=="top"){
k=sample(1:(m-1),size=1)
n=sample((k+1):m,size=1); if(k==m-1) n=m
if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])
if(n==m) N=c(N[(k+1):n],N[1:k])}
if(card=="bottom"){
k=sample(2:m,size=1)
n=sample(1:(k-1),size=1); if(k==1) n=1
if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])
if(n==1) N=c(N[k:m],N[n:(k-1)])}

Again, it is possible to have all those codes in the same function,

schuffle2=function(m,ns=10,which="top"){
N=1:m
if(which=="top"){
for(i in 1:ns)
{
k=sample(1:(m-1),size=1)
n=sample((k+1):m,size=1); if(k==m-1) n=m
if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])
if(n==m) N=c(N[(k+1):n],N[1:k])
}}
if(which=="bottom"){
for(i in 1:ns)
{
k=sample(2:m,size=1)
n=sample(1:(k-1),size=1); if(k==1) n=1
if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])
if(n==1) N=c(N[k:m],N[n:(k-1)])
}}
if(which=="mixed"){
for(i in 1:ns)
{card=sample(c("top","bottom"),size=1)
if(card=="top"){
k=sample(1:(m-1),size=1)
n=sample((k+1):m,size=1); if(k==m-1) n=m
if(n<m)  N=c(N[(k+1):n],N[1:k],N[(n+1):m])
if(n==m) N=c(N[(k+1):n],N[1:k])
}
if(card=="bottom"){
k=sample(2:m,size=1)
n=sample(1:(k-1),size=1); if(k==1) n=1
if(n>1)  N=c(N[1:(n-1)],N[k:m],N[n:(k-1)])
if(n==1) N=c(N[k:m],N[n:(k-1)])
}
}}
return(N)}
• How long should we shuffle before having cards randomly sorted ?

With the codes mentioned above, it is possible to run generations of shuffles,

distu=function(k=100,j=32){
U1B=U1T=U1M=U2B=U2T=U2M=U3=NULL
for(s in 1:100){
X1T=(0:(j-1))[schuffle1(j,k,"top")]
X1B=(0:(j-1))[schuffle1(j,k,"bottom")]
X1M=(0:(j-1))[schuffle1(j,k,"mixed")]
X2T=(0:(j-1))[schuffle2(j,k,"top")]
X2B=(0:(j-1))[schuffle2(j,k,"bottom")]
X2M=(0:(j-1))[schuffle2(j,k,"mixed")]
X3 =(0:(j-1))[sample(1:j)]

x1T=sum(j^(-(1:j))*X1T)
x1B=sum(j^(-(1:j))*X1B)
x1M=sum(j^(-(1:j))*X1M)
x2T=sum(j^(-(1:j))*X2T)
x2B=sum(j^(-(1:j))*X2B)
x2M=sum(j^(-(1:j))*X2M)
x3 =sum(j^(-(1:j))*X3)

U1T=c(U1T,x1T)
U1B=c(U1B,x1B)
U1M=c(U1M,x1M)
U2T=c(U2T,x2T)
U2B=c(U2B,x2B)
U2M=c(U2M,x2M)
U3 =c(U3,x3)
}
B=list(U1T=U1T,...)
}

and then, we run tests to see if the samples can be assumed to be uniformly distributed on the unit interval, e.g. for the very first kind first shuffle describe above, it would be

ks.test(B$U1T,punif)$p.value

More precisely, we use the following function, to estimate to proportion of scenarios where the $p$-value exceeds 5%,

PV=function(k){
P1B=P1T=P1M=P2B=P2T=P2M=P3=NULL
for(i in 1:10000){
B=dist(k,j=32)
P1T=c(P1T,ks.test(B$U1T,punif)$p.value)
P1M=c(P1M,ks.test(B$U1M,punif)$p.value)
P2T=c(P2T,ks.test(B$U2T,punif)$p.value)
P2M=c(P2M,ks.test(B$U2M,punif)$p.value)
P3 =c(P3,ks.test(B$U3,punif)$p.value)}
return(list(
p1T=mean(P1T>.05),
p1M=mean(P1M>.05),
p2T=mean(P2T>.05),
p2M=mean(P2M>.05),
p3=mean(P3>.05)))}

If we plot the results, we have

K=1:100
MP=Vectorize(PV)(K)
plot(K,MP[1,],col="red",type="b",ylim=0:1,pch=1)
lines(K,MP[2,],type="b",pch=19,col="red")
lines(K,MP[3,],col="blue",type="b",pch=1)
lines(K,MP[4,],type="b",pch=19,col="blue")
lines(K,MP[5,],type="b",pch=3,col="black")

Here, we look at the proportion of $p$-values that exceed 5%. We can pretend that we have a uniform distribution if that proportion is close to 95%. So basically, we just have to see when we reached for the first time the 95% region.If we zoom in the upper part of the graph, we get

With 32 cards,

• with a top in at random, we have to shuffle about 70 or 80 cards before having a randomly sorted set of cards. Which is large, but which is quite intuitive. One can imagine that it might take a while before getting the cards at bottom much higher in the pack,
• with a randomly mixed top in at random strategy, it is faster, slightly (we do not have that problem with cards at bottom that stay at bottom), since it takes about 60 or 70 cards.
• with a tops in a random, it goes again faster, with about 35 rounds,
• with a randomly mixed tops-bottoms in at random, it takes about 10 to 15 rounds.

Those results were obtained on tests on samples of size 100. The same code ran on a server during the week-end, with samples of size 500. Note that the output is rather close,

Note that those algorithm were mentioned because they were feasible, not only from a computational point of view, but when playing with real cards, in paper. Like with kids. I can actually ask my kids to perform those shuffle techniques next time we play with cards. The good thing is that randomly mixed tops-bottoms in at random shuffle technique: kids can do it 10 times, and cards should be randomly ordered in the deck…

Now, for those willing to see more algorithms, there are the so-called Fisher-Yates (also Knuth) shuffle. But may I keep that for another post ?

# From a random generator to a sample function

This week-end, I wrote a post since I had some trouble to generate a sample random sample with R, to reproduce one obtained by a co-author, with SAS (generated using Fishman and Moore (1982) used in function RANUNI). I was lucky since another contributor for that book, Christrophe Dutang, got the anwer to the last question I asked: is it possible to reproduce the random generator ? Yes, we can. And it is quite simple, if you use the appropriate library and the appropriate function,

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123)
> runif(10)
[1] 0.7503961 0.3209120 0.1783896 0.9060334 0.3571171
[6] 0.2211140 0.7864383 0.3980819 0.1246652 0.1876858

If you check in the previous post this is exactly what SAS gave us (and that I could not reproduce by myself). But that was only one part of my problems, since the goal was actually to reproduce indices for a training subsample for credit scoring issues.

I have to admit that I had never though about it before: how should we write a sample function? If values can be replaced, that is fine, we just have to split the unit interval correctly. Like

> set.seed(95)
> (U=runif(10))
[1] 0.15171155 0.57584087 0.05309844 0.07044485 0.48887914 0.15276707
[7] 0.37405684 0.30006096 0.96997126 0.30373498
> set.seed(95)
> (R=sample(0:99,size=10,replace=TRUE))
[1] 15 57  5  7 48 15 37 30 96 30

Here, we just truncate from the values obtained from the random generator. And that is just fine. But how do we write a code to sample without replacements? I mean, how do you get that :

> (S=sample(0:99,size=10,replace=FALSE))
[1] 15 57  5  6 46 14 35 27 89 92

My initial idea was to use the following technique. The first value is easy to get: we just split the unit interval into 100 subdivision (as before since for the first value, replacement or not, we don’t care) and see in which interval the random value is. And we remove that value from our sample. Then, we split the unit intervall into 99 subdivision, and see in which interval the random value is. It is the 10th? Fine, then our second value is the 10th from our sample (the first value has been removed). Then we split the unit interval in 98 subdivision, etc. The code I wrote to produce that algorithm was the following,

> mysample1=function(N,unif){
+  n=length(N)
+  size=length(unif)
+  V0=N[trunc(unif[1]*n)+1]
+  N=N[-which(N==V0)]
+  V=V0
+  for(i in 2:length(unif)){
+    V0=N[trunc(unif[i]*(n-i+1))+1]
+    N=N[-which(N==V0)]
+    V=c(V,V0)}
+ return(V)}

Unfortunetely, I could not reproduce the sample obtained with the R function,

> mysample1(0:99,unif=U)
[1] 15 58  5  7 49 17 39 31 97 32
> S
[1] 15 57  5  6 46 14 35 27 89 92

Since Christrophe is an expert on random generators, I did ask him, one more time. And he came up with the following code,

> mysample2=function(N,unif){
+   integerset=1:length(N)
+   result=rep(NA,length(unif))
+   for(i in 1:length(unif)){
+     intchosen=integerset[ceiling(U[i]*(length(N)-i+1))]
+     integerset[intchosen]=length(integerset)
+     integerset=integerset[-length(integerset)]
+     result[i]=intchosen}
+   return(N[result])}

which works just fine !

> mysample2(0:99,unif=U)
[1] 15 57  5  6 46 14 35 27 89 92
> S
[1] 15 57  5  6 46 14 35 27 89 92

So now, not only can we reproduce random numbers obtained with other software, we can also obtain the same samples indices, with or without replacement ! Thanks Christophe !

[May, 15th] Note that this is note the generator used in SAS, actually. In order to reproduce the sample function of SAS, the algorithm is much more simple, by clearly not that efficicient since we generate a random sample of size 100 (if we have 100 observations), and then, we keep the values associated to the indices of the 10 smallest (if we want a sample of size 10). The code could be

> mysample3=function(N,unif,size){
+ q=sort(unif)[size]
+ return(N[U<=q])}

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123) #OK

> U=runif(100)
> mysample3(1:100,U,size=10)
[1] 27 37 47 59 60 71 75 82 84 87

Thanks Jean-Philippe for the idea (which works).

# Reproducibility and randomness

With Stéphane Tufféry, we were working this week on a chapter of a book, entitled Statistical Learning in Actuarial Science. The chapter should be based on R functions, and we wanted to reproduce some outputs he previously obtained with SAS. The good thing is that even complex functions (logistic regression, regression trees, etc) produce the same kind of outputs. But we found a problem that we could not fix: generating identical training subsets of observations… Out of 1,000 lines, we subsample about 600 lines. The problem is that we could not generate the same sets of indexes with R, and SAS. Even using similar random generators… (execpt if we want to extract 1 or 2 lines, no more).

Let us try to explain what’s going on (based on code produced by Stéphane). According to Eubank (2010), there are (at least) two generators of random numbers in SAS,

For instance, for the RAND function, if we generate a Gaussian sample – with Mersenne-Twister algorithm – the code should be

%LET SEED =6;
%LET NREP=10;
DATA TESTRANDOM ;
CALL STREAMINIT(&SEED);
DO REP = 1 TO &NREP;
X = RAND ('NORMAL');
OUTPUT;
END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;

And we get here

 Obs. REP X 1 1 2.10680 2 2 -0.25604 3 3 0.28692 4 4 -0.22806 5 5 1.34569 6 6 0.16341 7 7 -0.27788 8 8 0.02133 9 9 1.24050 10 10 1.01054

If we want a Uniform sample, it should be

%LET SEED =6;
%LET NREP=10;
DATA TESTRANDM ;
CALL STREAMINIT(&SEED);
DO REP = 1 TO &NREP;
X = RAND ('UNIFORM');
OUTPUT;
END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;
 Obs. REP X 1 1 0.66097 2 2 0.48044 3 3 0.87849 4 4 0.19916 5 5 0.04838 6 6 0.19966 7 7 0.81353 8 8 0.53807 9 9 0.01105 10 10 0.53753

On good thing (so far) about the latest, is that Mersenne-Twister has been coded in R, in the RNGkind function

> RNGkind("Mersenne")
> set.seed(6)
> runif(10)
[1] 0.64357277 0.91590888 0.09523258 0.29537280
[5] 0.76993170 0.25589353 0.51789573 0.67784993
[9] 0.14722782 0.70052604

But the output is different, even if we’re supposed to start, here, with the same seed. Now, if we want to make sure about what is done here, let us write our own codes of the Fishman and Moore (1982) algorithm (in order to reproduce the SAS output). The R version of

> a = 397204094      # RANUNI multiplier
> seed = 123         # seed
> n = 10             # sample size
> m = (2^31) - 1     # period
> for (i in (1:n-1)) {
+ seed = (a*seed)%%m
+ u = seed / m
+ print(u)
+ }
[1] 0.7503961
[1] 0.3209121
[1] 0.3453204
[1] 0.2683455
[1] 0.241798
[1] 0.9888181
[1] 0.3943279
[1] 0.9710172
[1] 0.001632214
[1] 0.921537

Let us now run a similar code with SAS,

%LET SEED =123;
%LET NREP=10;
DATA FRANUNI (KEEP = x) ;
seed = &SEED ;
DO REP = 1 TO &NREP;
CALL RANUNI(seed, x);
OUTPUT;
END;
RUN;
PROC PRINT DATA = FRANUNI ;
RUN ;

and we get the following output

 Obs. x 1 0.75040 2 0.32091 3 0.17839 4 0.90603 5 0.35712 6 0.22111 7 0.78644 8 0.39808 9 0.12467 10 0.18769

It looks like here, indeed, we start with the same seed, since the first two numbers generated are similar. But then, it looks like we really have random numbers… If we change the seed, the first two numbers are similar, but that’s all.

We might be missing something trivial here, but we did not see it. So if anyone has a clue about reproducibility issues when generating random samples, with R and SAS, we are interested !

# In three months, I’ll be in Vegas (trying to win against the house)

In fact, I’m going there with my family and some friends, including two probabilists (I mean professionals, I am merely an amateur), with this incredible challenge: will I be able to convince  probabilists to go to play at the Casino?

Actually, I also want to study them carefully, to understand how we should play optimally. For example, I hope I can make them play the roulette. Roulette is simple. With a French (or European) roulette, it is probably the simplest: if I bet on black, I win if one of 18 black numbers is out, and I lose if one of the 18 red numbers – or zero (which is green) – is out. This gives a winning probability of 18/37 i.e. a 48.64% chance. But in Vegas, I think it’s mostly American Roulette that can be found in casinos, in which there is a zero and a double zero (both favorable to the bank). Here, the  probability of winning is 18/38, i.e. 47.36% chance. The two roulettes are

Now, let us discuss a little bit about optimal strategy. For instance, suppose I go to Las Vegas with an initial wealth $s$ (say $100). The goal is to find the strategy which maximizes the probability to leave Las Vegas with $2s$ (here$200). Should I play big, or small ?

Assume that I can bet $x$ (that will be, here, for convenience, a fraction of $s$). With probability $p$, I will get $2x$, and with probability $1-p$, I will get $0$ (and lose my $x$). As mentioned above, $p$ is (a little) smaller than 50%. The casino must win (actually, we will see that this assumption has a very strong impact on the strategy).

Suppose my goal is to double my initial sum, as mentioned in the introduction of this post. Maybe there is an optimum value for $x$, to maximize the probability of doubling my bet. To make it simple, the game ends either because I did not, or because I did, manage to double my initial wealth… Assume further that $x$ is fixed, and that I do not revise my bets. One can use monte carlo simulations, to get an intuitive idea…

> bet=function(s=1,t=2*s,x=s/4,p=.4736,nsim=100000){
+     vp=rep(0,nsim); #vw=s
+     for(i in 1:nsim){
+       w=s;
+       while((w>0)&(w<t)){
+          ux=sample(c(min(x,t-w),-x),size=1,prob=c(p,1-p))
+          w=w+ux
+       }
+       vp[i]=(w>=t)}
+     return(mean(vp))
+ }

If we plot this probability as a function of $x/s$, we have the following

> BET=function(x) bet(x=x)
> vx=1/(1:20)
> px= Vectorize(BET)(vx)
> plot(vx,px,log="x")


Let us see if we can do the maths, and actually compute those probabilities.

For example, if $x = s$, I play everything I have, and I double with probability $p$. That one was simple.  And indeed, on the graph above, the point on the right is probability  $p$ (the red horizontal line).

Assume now that I can bet $x = s / 2$, and I will play, at least, two rounds

• with probability $(1-p) ^ 2$ I will lose both rounds (and the game is over)
• with probability $p ^ 2$, I will win both rounds, and I double my bet (and the game is also over)
• with probability $2p (1-p)$, I will lose once, and double once. Anyway, I will find myself again with my (initial) wealth $s$. So the game will start again….

To make the story short the probability of doubling my earnings is

$p ^ 2+ 2 p (1-p)\big( p^ 2 + 2 p (1-p)\big( \cdots$

which is

$p ^ 2 \left (1 +2 p (1-p) + [2p (1-p)] ^ 2 + \cdots \right) = \frac {p ^ 2} {1-2p (1-p)}$

Let’s try something more general: I have initial wealth $s$, I can bet $x$ and the goal is to reach $2s$ (or, more generally, say, $t$). Now, the probability to reach $t$ from $s$ betting (always) $x$ is exactly the same as the probability to reach $t/x$ from $s/x$ betting only 1. Let $P_b(a)$ denote the probability to go from $a$ to $b$ betting 1 (let us use generic parameters). We can easily get the following equation

$P_b(a) = p\cdot P_b(a+1) + (1-p) \cdot P_b(a-1)$

Thus, we can write

$p\cdot (P_b(a+1)-P_b(a)) = (1-p)\cdot (P_b(a)-P_b(a-1))$

or equivalently

$(P_b(a+1)-P_b(a)) =\frac{1-p}{p}\cdot (P_b(a)-P_b(a-1))$

$\left(\frac{1-p}{p}\right)^a\cdot (P_b(1)-P_b(0))$

Now, observe that $P_b(0)=0$ (since I cannot have a gain without any money).

Let us write $P_b(a+1)-P_b(0)$ using a domino technique :

$[P_b(a+1)-P_b(a)]+[P_b(a)-P_b(a-1)]+\cdots+[P_b(1)-P_b(0)]$

i.e.

$\left(\frac{1-p}{p}\right)^a P_b(1)+\left(\frac{1-p}{p}\right)^{a-1} P_b(1)+\cdots+ \left(\frac{1-p}{p}\right)^0 P_b(1)$

so this geometric sum can also be written

$\left(1 -\left[\frac{1-p}{p}\right]^{a+1} \right) \left(1 -\left[\frac{1-p}{p}\right] \right)^{-1}$

Finally, we can write

$P_b(a)=\left(1 -\left[\frac{1-p}{p}\right]^{a} \right)\left(1 -\left[\frac{1-p}{p}\right] \right)^{-1}\cdot P_b(1)$

Here, there is still $P_b(1)$ that I have to explicit. The idea is to observe that $P_b(b)=1$, thus

$P_b(a)=\left(1 -\left[\frac{1-p}{p}\right]^{a} \right)\left(1 -\left[\frac{1-p}{p}\right]^{b} \right)^{-1}$

So finally,

$\mathbb{P}(gain)=\left(1 -\left[\frac{1-p}{p}\right]^{s/x} \right)\left(1 -\left[\frac{1-p}{p}\right]^{2s/x} \right)^{-1}$

Nice isn’t it? But to be honest, there is nothing new here. This is actually an old theorem discovered by Christiaan Huygens in 1657, then extended by Jacob Bernoulli in 1680 and finally properly established by Abraham de Moivre in 1711. It is possible to plot this graph, as a function of $x/s$,

> bet2=function(s=1,t=2*s,x=s/4,p=.4736){
+     vp=(1-((1-p)/p)^(s/x))/(1-((1-p)/p)^(t/x))
+     return(vp)
+ }

The graph is the same as the one with monte carlo simulation (hopefully). Observe, looking carefully at the function above, that the probability is decreasing with $p$. Which makes sense… Further, the probability is decreasing with $t$: the more hungry, the less chance of winning I have.

Now, the interesting part is what is plotted on the graphs above: the smaller $x$ (the size of the bets at each round), the less chances to win: if I want to win, it is important not to play being little player ! I must bet everything I have ! Actually, the funny thing is that if the probability of winning was (slightly) larger than 1/2, on the contrary, I should bet as small as possible

So far, there is nothing new. Everything mentioned in this post can be related to a fundamental result of Lester Dubins and Leonard Savage, in “How to Gamble if You Must : Inequalities for Stochastic Processes” (published in 1965), see also Sudderth (1972). Of course, I can try another strategy, a little less reasonable, I think, which is sometimes called Martingale of D’Alembert. I believe more in luck than coincidence, so, when I win, I drop my bet (do not tempt fate) but when I lose, I increase my bet (I must win someday). But let’s keep it for another post, someday…

Again, that’s a theory. I guess we should try, and see how it works. I’ll try to upload pictures on the blog during the road trip, so if by the beginning in August nothing has been posted on the blog, please send a rescue team to save me at the Bellagio…

# Martingale et journalisme scientifique

Être journaliste scientifique ne doit pas être facile. J’imagine qu’il faut être à l’écoute des nouvelles scientifiques, et d’informer, aussi justement que possible. C’est ce que fait avec brio Pierre Barthélémy (aka @PasseurSciences sur Twitter) dans sa chronique hebdomadaire sélection scientifique de la semaine. Il essaye ainsi de parler de sciences dans un journal qui a (trop) souvent confondu science et technologie (technologie étant aussi souvent un mot savant utilisé pour masquer de la publicité pour des appareils technologiques avancés). Et comme l’espace qui lui est imparti est restreint (c’est le moins qu’on puisse dire), la plupart des informations passent par son blog http://passeurdesciences.blog.lemonde.fr/, qui est riche ! Sylvestre Huet fait aussi ce travail pour Libération, avec son blog http://sciences.blogs.liberation.fr/ (là encore, on peut lire bien plus de choses en ligne que que dans la version papier, ce qui est d’autant plus intéressant pour les lecteurs de l’autre coté de l’océan atlantique). Sylvestre Huet arrive à parler géologie, puis démographie, avant d’évoquer dans un article truffé de références les réformes des instances qui régentent la recherche en France. Oui, les deux blogs sont malgré tout très français. Pour une vision plus internationale, on pourra citer l’admirable travail de Pascal Lapointe, par exemple (aka @paslap) sur http://sciencepresse.qc.ca/.

Mais si je commence à mettre les pieds dans le plat, je poserais la question de la légitimité d’un journaliste scientifique, à l’heure où la blogosphère scientifique explose. Par exemple, si je veux apprendre des choses en physique, ou en biologie, je sais que le blog de Tom Roud (que l’on peut suivre sur Twitter @tomroud) http://tomroud.cafe-sciences.org est une source infinie d’information (ok, peut-être pas “infinie” car comme nous tous, il a une vie en dehors de son blog, mais disons qu’il y a de quoi lire). En sciences humaines, il existe des centaines de blogs hébergés sur http://hypotheses.org/ tenus par des universitaires, ou http://www.cafe-sciences.org/ pour des blogs de sciences, en français. Les bloggers scientifiques s’expriment – le plus souvent – en restant dans leur champ d’expertise. Et c’est tant mieux. Le journaliste scientifique, lui, se doit de parler de tous les sujets scientifiques. Et c’est là que l’exercice devient délicat, car le journaliste se doit d’être critique. Compte tenu de la concurrence féroce qui existe dans le monde académique, on nous pousse à faire croire qu’on vient d’écrire l’article qui va révolutionner la science, sinon le monde. Que notre approche est novatrice, et qu’en plus, on vient de montrer ce qui pourrait ce qui pourrait être le sain Graal dans notre communauté. C’est le travail du journaliste de se demander si c’est vrai….

J’en arrive au point de mon billet. Le dernier article  le hasard, martingale boursière ? (en ligne sur http://lemonde.fr/sciences/…) de Pierre Barthélémy m’a un peu agacé. Pas sur le thème abordé, mais l’impression générale qui m’est restée après avoir lu l’article (samedi matin). Pour l’histoire complète, c’est ma femme qui est tombée la première sur l’article, et qui a été étonnée de lire un article pareil dans Le Monde, en tant que mathématicienne (et probabiliste). Et quand je l’ai lu, c’est en tant qu’économiste que je suis resté sans voix. En résumé (trop succinct, j’en convient), on nous résume l’article de physiciens italiens qui “expliquent qu’au grand dam des traders et autres analystes financiers, en quête perpétuelle de justifications rationnelles aux fluctuations des cours boursiers, les marchés demeurent obstinément imprévisibles“. Pour parcourir l’article mentionné, mais non cité – probablement car il s’agit de la version papier ? – l’article est en ligne sur http://arxiv.org/1303.4351  (trouvé via http://improbable.com/…). On retrouve cette idée dans le texte sous la forme suivante “our main result, which is independent of the market considered, is that standard trading strategies and their algorithms, based on the past history of the time series, although have occasionally the chance to be successful inside small temporal windows, on a large temporal scale, perform on average not better than the purely random strategy, which, on the other hand, is also much less volatile. In this respect, for the individual trader, a purely random strategy represents a costless alternative to expensive professional financial consulting, being at the same time also much less risky, if compared to the other trading strategies.” On apprend dans l’article paru dans Le Monde que “pour Pluchino et compagnie, le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre“. Damned, rien que ça ? Pour les amateurs de lyrisme (il doit y en avoir dans les lecteurs du monde, moins parmi les lecteurs plus scientifiques qui aiment les phrases moins alambiquées), je passe les moments où Pierre Barthélémy est clairement plus journaliste que scientifique, comme lorsqu’il nous explique que ces quatre chercheurs “viennent d ‘attaquer à l’acide de l’aléatoire” les modèles financiers. Cela dit, cela a le mérite de mettre les points sur les i: ce sont des physiciens. Pas des économistes.

Je pense qu’il y a des milliers de raisons d’attaquer les traders (si ce sont bien les “spécialistes” visés). Mais qu’on ne les prenne pas pour plus bêtes qu’ils ne le sont (je le dis d’autant plus volontiers que j’en ai formé un paquet, à l’ENSAE ou à Polytechnique, au cours des 10 dernières années).

Avant de revenir sur la petite histoire de martingale financière (tel que le résume le titre de l’article) faisons un court détour par celle brillamment racontée dans Mansuy (en ligne sur http://math.harvard.edu/~ctm/… et traduit en anglais dans http://jehps.net/Mansuy.pdf) sur l’origine de la notion de martingale. D’un point de vue mathématique (disons comme propriété des processus stochastiques), il faut remonter aux travaux de Joseph Bernstein, Paul Lévy, Émile Borel, et surtout Joseph Leo Doob, au milieu au XXième siècle. Cela dit, le mot est plus ancien: il entre dans le dictionnaire de l’Académie Française en 1762. Jouer à la martingale, c’est jouer toujours tout ce que l’on a perdu.  Mais quelques années plus tot, on pouvait déjà trouver le mot sous la plume de l’Abbé Prévost (le jeu qu’il décrit en 1750 comme variante du jeu du pharaon est aussi appelé martingale d’Alembert) Cela dit, l’origine que je préfère est de relier le mot martingale à l’expression provençale “jouga a la martegalo” qui signifierait jouer de manière incompréhensible, absurde, comme l’évoque Frédéric Mistral dans Lou Tresor dòu Felibrige ou dictionnaire provençal-français. On retrouve une origine proche dans le dictionarie of the French and English tongues de Randle Cotgrave, datant de 1611, qui mentionne l’expression “à la martingale” avec le sens absurdly, foolishly, untowardly, grossely, rudely, in the homeliest manner. Il cite même l’usage de l’expression philosopher à la martingale (sans citer Bernard Henri Levy, mais il semble que ce soit l’idée). Cela dit, même dans Lapinot on parle de martingale,

Bref, à partir de ce comportement incompréhensible, voire absurde, des économistes vont définir une notion importante en finance, que l’on appellera efficience (qui est un mot dangereux car il évoque aux oreilles de tous les économiste une notion d’optimalité au sens de Pareto, on pourra relire Malkiel (2003) The Efficient Market Hypothesis and Its Critics sur ce point). L’idée d’utiliser une marche aléatoire pour modéliser le cours d’un actif est ancienne. Par exemple (je ne remontrais pas au XVIème siècle, mais http://e-m-h.org/history.html le suggère) on pourra relire les travaux de Bachelier datant du début du XXième siècle, mais comme le dit la légende, Bachelier a été peu lu, à l’époque (en tous les cas par des économistes). Publié quelques années plus tard, on pourra relire Cowles & Jones (1937), Some A Posteriori Probabilities in Stock Market ou surtout Samuelson (1965) Proof That Properly Anticipated Prices Fluctuate Randomly. Cet article est passionnant, si on prend le temps de le lire “There is no way of making an expected profit by extrapolating past changes in the future price, by chart or any other esoteric devices of magic or mathematics.” Paul Samuelson est d’une modestie remarquable pour un chercheur, “I have not here discussed what the basic probability distributions are supposed to come from. In whose mind are they ex ante? Is there any ex post validation of them? Are they supposed to belong to the market as a whole? And what does that mean? Are they supposed to belong to the “representative individual”, and who is he? Are there some defensible or necessitous compromises of divergent anticipations patterns? Do price quotations somehow produce a Pareto-optimal configuration of ex ante subjective probabilities? This paper has not attempted to pronounce on these interesting questions“. Il peut y avoir des bulles, des processus non Gaussien, un peu tout ce qu’on peut imaginer, le point important étant qu’il est impossible d’utiliser le passé pour prédire le futur, en finance.

Oui, depuis 50 ans, les économistes savent que la marche aléatoire peut être un bon modèle pour décrire le prix des actifs. Ou pour être plus précis (et rendre à Paul Samuelson ce que Paul Samuelson a dit le premier) les martingales. Mais il faudra surtout attendre les travaux d’Eugène Fama, en particulier Fama (1965) The Behavior of Stock Market Prices et Fama (1970) Efficient Capital Markets: A Review of Theory and Empirical Work, qui vont de manière définitive marquer le début de l’utilisation des martingales en finance. Eugène Fama y retient trois notions d’efficience (pour aller plus loin que le fameux “A market in which prices at any time “fully reflect” available information is called “efficient”“)

1. Expected Returns or “Fair Game” Models
2. Submartingale Models
3. Random Walk Model

(en 1965, l’efficience devait être reliée, pour Eugène Fama, à la notion de marche aléatoire, alors que Paul Samuelson utilisait déjà la notion de martingale). Ces modèles disent tous (on est en 1970) qu’il est inutile d’utiliser l’information passée ou présente “to predict the future in a way which makes expected profits greater than they would be under a naive buy-and-hold model“. N’en déplaise à Pierre Barthélémy , depuis presque 50 ans ce principe est énoncé dans les cours d’asset pricing. Bachelier disait la meme chose en 1900, “l’espérance mathématique du spéculateur est nulle” dans Théorie de la spéculation. Et cela sera confirmé quelques années plus tard par les résultats empiriques de Cowles et Jones. Un an plus tard sera publié Black (1971) Implications of the random walk hypothesis for portfolio management, et trois ans plus tard, Hagerman & Richmond (1973) Random Walks, Martingales and the OTC. La littérature économique va se multiplier pendant les années 70, à tel point qu’en 1978, Michael Jensen, alors professeur à Harvard écrivait “I believe there is no other proposition in economics which has more solid empirical evidence supporting it than the Efficient Market Hypothesis” (cité dans http://economist.com/14030296). Plus récemment, à la fin des années 80 était publié le remarquable LeRoy (1989) Efficient Capital Markets and Martingales, qui étoffait la critique initiée dans LeRoy (1976). Efficient capital markets: Comment. LeRoy et Samuelson ont été les premiers à parler de martingale pour modéliser les prix des actifs. On continue (avec les articles qui font référence, et que tous les étudiants qui ont suivi un cours d’asset pricing ont lu) ? L’année suivante était publié Lehman (1990) Fads, Martingales, and Market Efficiency. Les martingales sont des outils incroyablement riches. Ils ne sont pas équivalent à une marche aléatoire: on peut avoir une martingale, sans avoir de marche aléatoire (par exemple avec un processus ARCH, on pourra relire sur le sujet les notes de cours Predictability of Asset Returns, en particulier le premier paragraphe, ou en français, avec des cours de licence, par Francis Diener,  http://math.unice.fr/~diener/…). Pour ceux qui veulent moins de formalisme, une martingale (en simplifant outrageusement) c’est, comme l’explique Nicolas Poupon,

Mais ces articles posent essentiellement les bases de la théorie financière. Qu’en est-il des aspects empiriques ? Andrew Lo et Craig MacKinlay ont publié un livre, en 2001, qui recense plusieurs articles sur le sujet, A Non-Random Walk Down Wall Street. Qui remet en cause l’idée de marche aléatoire, moins celle de marginale. On pourra aussi penser à Beechey Gruen & Vickrey (2000) et leur étude The Efficient Markets Hypothesis: A Survey, en ligne sur la Reserve Bank of Australia. Enfin, Pour une méthodologie plus proche de celle évoquée dans l’article en ligne sur arxiv, on pourra lire Martingales, the Efficient Market Hypothesis, and Spurious Stylized Facts de Joseph McCauley, Kevin Bassler et Gemunu Gunaratne (si on souhaite aussi utiliser des propriétés de mémoire longue). Moralité ? Non, ce n’est pas nouveau que l’on teste l’hypothèse d’efficience et de martingale. Et si je cite des articles académiques, il faut  ajouter que les journalistes économiques connaissent également cette littérature. Il y a 10 ans, Justin Fox posait la même question que nos chercheurs dans un article dans Fortune intitulé Is The Market Rational ? No, say the experts. But neither are you–so don’t go thinking you can outsmart it. On pourra aussi relire Efficiency and beyond (mentionné auparavant) qui revient sur la notion d’efficience des marchés. On y retrouve que les spécialistes savent tout ca: “on such ideas, and on the complex mathematics that described them [i.e martingales], was founded the Wall Street profession of financial engineering“. Enfin, en 2009, Richard Thaler faisait quelques rappels sur la difficulté d’interpréter correctement l’hypothèse d’efficience des marchés dans Markets can be wrong and the price is not always right. Pour conclure sur l’histoire des martingales en finance, en 2010, Fama répondait d’ailleurs de manière très claire à la question dans une interview – interview with Eugene Fama – lorsque John Cassidy lui demandait “the fundamental insight of the efficient market hypothesis [is] that you can’t beat the market ?“, et qu’il répondait “Right—that’s the practical insight. No matter what research gets done, that one always looks good“.

Pour revenir sur l’article de Pierre Barthélémy, je trouve qu’écrire “le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre” est incroyablement méprisant envers des centaines, pour ne pas dire des milliers de chercheurs en mathématiques financières, et en économétrie de la finance. Mais rassure toi Pierre, je vais continuer à lire tes chroniques. C’est juste que ma grand mère est encore abonnée au Monde, et j’aimerais bien qu’elle continue de croire que je peux sauver l’humanité (et non pas que je refuse d’admettre ce que des chercheurs italiens viennent d’établir, même si cela confirmerait – à ses yeux – que je peux etre têtu comme une mule quand je veux).

Maintenant, pour conclure sur l’article mentionné, et sur l’analyse faite, on ne peut pas dire que les conclusions soient révolutionnaires. Disons, pour quiconque ayant lu un peu de littérature économique dans les 50 dernières années. En fait, pour être honnête, je pense qu’on peut dire que Alessandro Pluchino et Andrea Rapisarda se sont amusé (on inclura Biondo et Helbing). Sans plus de prétention. Pour ceux qui ne se souviennent pas, ils avaient eu un IGNobel en 2011, pour leur étude sur organizations would become more efficient if they promoted people at random. J’avais aussi beaucoup aimé leur article sur les aspects computationnels du principe de Peter, en ligne sur http://arxiv.org/0907.0455. Quand je lis leurs études, j’ai l’impression qu’ils s’amusent. Et le terme n’est absolument pas péjoratif dans ma bouche, loin de là. C’est un peu ce que disait Lucile Quillet la semaine passée dans http://etudiant.lefigaro.fr/…  lorsqu’en mentionnant une étude de Baptiste Coulmont, elle écrivait “c’est le jeu auquel c’est amusé Baptiste Coulmont sur son blog“. Je ne sais pas si le terme était voulu, si c’était basé sur le fait qu’un blog, c’est fait pour jouer… mais j’ai toujours revendiqué que je faisais de la recherche pour mon plaisir, pour m’amuser. Et comme me le rappelait Baptiste, le mot schola signifie école, en latin, et son étymologie grecque est le mot σχολή, signifiant loisir.

# 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….

# 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/
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 ! # De la difficulté de générer du hasard J’avais parlé cet été de l’utilisation des chiffres du us treasury balance comme générateur (involontaire) de nombres aléatoires. Cet été, un débat intéressant a secoué le Québec, suite à la plainte déposée par deux joueurs (compulsifs) de l’Extra. L’histoire est raconté en détails sur le site de canoe.ca, ou de lapresse.ca. Le jeu est assez simple. L’histoire un peu moins (c’est pour cela que j’ai préféré renvoyer vers deux sources racontant la même histoire). Pour essayer d’expliquer simplement ce qui s’est passé revenons un peu sur le jeu. Tous les jours, Loto-Québec tire au hasard un nombre de 7 chiffres. Comme on le voit ci-dessus, le but est d’avoir les premiers (ou les derniers) chiffres bons. Quand on achète un ticket, un numéro nous est attribué… et on espère qu’il sera tiré le soir. Mais le soucis est que les joueurs qui intentent un procès n’achetaient pas un mais dix tickets. Et étrangement, à chaque fois, ils avaient un ticket gagnant, et un seul. En fait, quand on achetait dix combinaisons d’extra, il semble que le premier chiffre sans remise parmi 0 à 9, les cinq suivants avec remise parmi 0 à 9 et le dernier sans remise parmi 0 à 9 (pour reprendre la formulation qu’a utilisé Christian Léger de l’UdM pour me l’expliquer). Sur un achat, on ne voit rien, on a des nombres uniforme… mais avec répétition, ça se complique. En fait, c’est toujours quand on étudie la dynamique des générateurs de nombres aléatoires que l’on a des surprises. Par exemple sur les tirages de l’extra (et non pas les numéros des combinaisons des billets, car on ne les trouve pas sur internet, il faudrait acheter les billets pour refaire l’expérience). On peut récupérer l’historique des tirages en ligne, sur le site de loto-quebec.com. > extra=read.table("http://freakonometrics.blog. free.fr/public/data/extra-loto-quebec.csv",sep=";",header=TRUE) > n=1940 > extra=extra[1:(n-1),1:2] > X=extra[,2]/10000000 Si on regarde les 1940 derniers tirages (car avant, il semble qu’il n’y avait pas 7 mais 6 chiffres), on a la série suivante, Et un test d’uniformité semble confirmer l’intuition que les nombres sont uniformément répartis, > ks.test(X,'punif') One-sample Kolmogorov-Smirnov test data: X D = 0.021, p-value = 0.3895 alternative hypothesis: two-sided Mais générer du hasard, ça n’est pas que ça… Et sur cette série, si on fait des tests sur l’indépendance des observations, on est un peu surpris, > library(lawstat) > runs.test(X, plot.it = FALSE, alternative = "two.sided") Runs Test - Two sided data: X Standardized Runs Statistic = -1.8, p-value = 0.07273 > runs.test(X, plot.it = FALSE, alternative = "positive.correlated") Runs Test - Positive Correlated data: X Standardized Runs Statistic = -1.8, p-value = 0.03637 > library(FitAR) > LjungBoxTest(qnorm(X), k=0, lag.max=12, + StartLag=1, SquaredQ=FALSE) m Qm pvalue 1 7.2 0.0075 2 9.2 0.0102 3 10.9 0.0125 4 10.9 0.0276 5 15.3 0.0090 6 15.4 0.0171 7 16.1 0.0244 8 16.3 0.0380 9 16.7 0.0538 10 20.6 0.0240 11 20.8 0.0357 12 20.8 0.0533 # 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 deslast 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://treasurydirect.gov/ et http://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/100000
X2=b2$DEBT Y2=trunc(as.numeric(X2)) X2=as.character(Y2) n2=nchar(X2) Y2=substr(X2,n2-4,n2) y=as.numeric(Y2)/100000 y=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 courbebleue, 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.01645 alternative hypothesis: two-sided > ks.test(y,"punif") One-sample Kolmogorov-Smirnov test data: y D = 0.0456, p-value = 0.3581 alternative 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. # Cursed numbers ? In Lost, Hugo “Hurley” Reyes played the numbers 4, 8, 15, 16, 23 and 42 at the lottery, and ended up winning the$114-million jackpot. And over the ensuing weeks, everyone around him seems to suffer increasingly bad luck: Hurley’s grandfather dies of a heart attack, his brother’s wife walks out, his mother breaks her ankle while the house Hurley bought her goes up in flames, and Hurley himself is falsely arrested.
Anyway, last week (here) 4 numbers (out of 6) appeared at the lottery in LA. As pointed out by Xian (here), the odds were not that small, i.e. it is a 1‰ chance,

Hence, with one lottery per week, the return period is 16 years. Note this percentage is very close to what we did observe on the French lottery (below the statistics in ‰, from here, in a zip file),

> loto=read.table("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")])
> P=rep(NA,nrow(N))
> for(s in 1:nrow(N)){
+ P[s]=sum(N[s,1]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,2]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,3]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,4]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,5]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,6]%in%c(4, 8, 15, 16, 23, 42))
+ }
> table(P)/nrow(N)*1000
P
0          1          2          3          4
435.732113 405.366057 137.271215  19.966722   1.663894

But what about the full sequence…? Imagine that in France, at the official lottery, the exact sequence played by Hugo appears. What a coincidence.The probability that the sequence appears, assuming that there are 48 possible numbers in the lottery, is

i.e. the expected number of draws we need before seeing that sequence for the first time is almost a billion.
Now if we look at all official lotteries around the world, say 100 per week, what is the probability to see Hurley’s sequence shows up – at least once – in 25 years (assuming that after 25 years, no one will remember Lost and those cursed numbers) ? It looks like it is a 1% chance…
So let us wait and see…

# Mandelbrot, fractals and counterexamples in applied probability

Benoît Mandelbrot died yesterday. Like most of the blogs dealing with applied mathematics, it looks like I have to mention this event. Unfortunately, I don’t know much about fractals…

The first time I heard about Mandelbrot and chaos was when I have been working on fractional time series (see eghere). Murad Taqqu gave a very interesting short course in Paris, and I have been using it in two papers (actually one more should appear soon in Climate Change).

The second time was in Québec (city), five years ago, when Roger Nelsen gave a talk on copulas with fractal suppport (here). By that time, we were finishing our paper with Alessandro (in mathematical finance, and limiting theorems). I remember adding a Remark in the paper (that can be found here), since using that kind of copulas was a nice way to show that, without sufficient regularity conditions, the limit we were looking for had no sense.

A few months after, with Johan, in a paper on pitfalls on lower tail dependence for Archimedean copulas (here), we used again this fractal construction (here applied to Archimedean copulas) to find a nice counterexample to a (false) theorem on regular variation. Again, the goal was to understand how the dependence structure of  given  and  changed when  goes to 0. In the case of Archimedean copulas, the copula of the conditional pair is still Archimedean (with another generator, except for Clayton copula). The graph below show how  changes, as  decreases… Actually, I draw  since Archimedean generators are not unique.

I have also decided to plot .

Here, we see that there is no way of talk about a possible limit for the conditional copula because of a fractal behavior in 0 of the generator (even if my fractal are not as nice as the one you can find on the internet….). So thanks Benoît for giving us a nice toy to build interesting counterexamples !