Category Archives: Randomness

Des modeles prédictifs en assurance

Cet post est aussi en ligne sur https://hal.archives-ouvertes.fr/hal-02350006, il a été coécrit avec Laurence Barry et Ewen Gallic.

Les compagnies d’assurance émettent des contrats qui prévoient des paiements d’indemnités en cas de survenance d’évènements aléatoires (accident, maladie, décès, etc.). En contrepartie, l’assuré doit s’acquitter d’une prime, dont le montant est déterminé ex-ante, avant le début de la période de couverture. Cette prime se décompose en deux termes : une prime pure (destinée à couvrir les pertes anticipées) et un chargement (incluant des commissions à des agents, divers frais, mais aussi couvrant contre le risque de variabilité des pertes). La prime pure est souvent calculée par classe de risque, et une classification est alors nécessaire.

Assurer une population hétérogène, ou l’importance de la classification

Le regroupement des risques selon diverses informations telles l’âge de l’assuré, son état de santé ou encore sa profession constitue ce que l’on appelle la classification des risques. Cette pratique de segmentation se justifie (à des fins d’admissibilité mais aussi de tarification) par la supposition que les risques sont placés dans des groupes relativement homogènes, au sein desquels les probabilités de survenance sont similaires. Pour Schauer (2006), cette « généralisation », qui vise à voir l’individu sous le prisme de sa classe de risque, de généraliser son comportement à partir de quelques variables explicatives, est probablement la raison d’être de l’actuaire : « to be an actuary is to be a specialist  in  generalization,  and  actuaries  engage  in  a form of decision-making that is sometimes called actuarial ». Statistiquement on cherche une méthode de classification aussi « discriminatoire » que possible[1], en gardant en mémoire que la discrimination est interdite, ce qui rend l’exercice périlleux, et souvent critiqué (nous y reviendrons plus loin).

Les assureurs évoquent souvent deux arguments pour justifier une segmentation. Le premier est qu’elle serait rendue économiquement nécessaire par la concurrence ; ne pas classifier conduit à une anti-sélection, les risques importants restant seuls chez les assureurs qui ne segmentent pas. Dans une telle situation, l’équilibre de marché ne serait pas possible puisque les risques faibles seraient chez un concurrent ayant segmenté. Si le facteur de risque était observable, tant par les assurés que les assureurs, il y aurait un phénomène d’auto-sélection, les assurés à risque faible ayant les polices les moins chères. Cette situation constitue un équilibre séparant de Nash. Mais si le facteur de risque n’est pas observable, un équilibre sous-optimal peut être atteint, résultant d’une externalité négative de cette information non-accessible, à la manière de Wilson (1977), tel que décrit dans Cummins et al. (1982) dans le cas des contrats d’assurance-vie. Cela dit, Kleindorfer & Kunreuther (1980) montrent qu’accéder à davantage d’information ne conduit pas nécessairement à une amélioration du bien-être des consommateurs. De plus si la classification n’est pas autorisée, l’équilibre est maintenu, les risques faibles subventionnant les risques élevés.

Le second argument avancé pour justifier une segmentation est que cette dernière (et le fait, par conséquent, d’ajuster les primes au risque) serait juste et équitable. Mais cette vision de l’équité n’a pas toujours été de mise et semble portée par les développements techniques. Ainsi la classification est devenue de plus en plus fine, multipliant les classes de risque et conduisant à des tarifs « personnalisés ». En plus des avancées statistiques, des facteurs économiques pourraient justifier cette sophistication : la concurrence de plus en plus forte sur certaines branches.

Incertitude en assurance

Il y a plusieurs manières de caractériser l’incertitude en assurance. Comme bien souvent quand on fait des prévisions, il convient de distinguer l’incertitude associée à l’estimation des probabilités et l’incertitude réelle sur le résultat (aléa de l’évènement). Pour la seconde notion, Hacking (1975) parle de probabilité structurelle, et c’est celle qui est souvent utilisée pour introduire les concepts de probabilité, par exemple avec des dés ou des jeux de cartes : les probabilités sont connues, seule l’issue du jeu est incertaine. Par exemple je sais que la probabilité d’avoir 6 en lançant un dé est 1/6 (compte tenu de la géométrie du cube).

D’un point de vue statistique, la probabilité se mesure quand on peut observer une fréquence, c’est-à-dire une répétition de risques semblables. Les statisticiens ont ainsi défini une notion de probabilité empirique, basée sur la répétition[2]. Si, en lançant mille dés j’obtiens 173 fois la face 6, la probabilité empirique d’avoir 6 est de 17,3%. La loi des grands nombres nous assure que cette fréquence va tendre vers la vraie valeur en répétant l’expérience, et le théorème central limite permet d’en contrôler les fluctuations. C’est la première incertitude dont nous parlions au début de cette section, que nous appellerions l’erreur d’estimation.

On peut enfin mentionner deux notions supplémentaires ; tout d’abord, les probabilités conditionnelles. Cette idée est introduite en assurance par de Moivre, ou de Witt, lorsqu’ils notaient que pour estimer une probabilité de décès, il fallait considérer des personnes de même âge. C’est cette idée que l’on retrouve quand on considère une classification : on veut des risques homogènes, similaires, sans être pour autant identiques. La probabilité que l’on obtient est alors conditionnelle à ce facteur commun qui caractérise la classe observée. Dans notre exemple des dés, cela revient à dire qu’il ne faut pas lancer mille dés, mais mille fois le même dé – ou à défaut des dés semblables.

Enfin, les probabilités subjectives ont été formalisées par Bruno de Finetti et Leonard Savage (ainsi que plus philosophiquement par Frank P. Ramsey) pour comprendre et modéliser la prise de décision. Elles sont relativement populaires en économie de l’incertain, mais difficile à mettre en œuvre dans un contexte de valorisation de contrats d’assurance automobile ou habitation. Il s’agit d’un jugement, qui ne peut être confronté à la réalité, mais envisageable pour l’assurance de risques encore mal connus (McGrayne (2012) évoque ainsi les premiers contrats d’assurance aviation). Une approche bayésienne consiste alors à combiner cette probabilité subjective avec la probabilité comme fréquence observée d’un phénomène : partant d’une croyance a priori, on affine l’estimation par une mise à jour progressive en répétant les expériences. Classiquement, la probabilité d’avoir la face 6 sera une moyenne entre notre croyance (1 chance sur 6) et une probabilité dite historique, obtenue en faisant quelques lancés (3 sur 20 lancers, par exemple). Les poids attribués aux deux dépendant du nombre d’expériences effectuées : on donnera plus de crédit à l’expérience si on fait mille lancés que si on en fait soixante.

Incertitude sur le résultat, ou aléa fondamental

Les probabilités prédictives, utilisées pour calculer la prime d’un contrat d’assurance, sont la première étape d’un problème de classification. Un outil classique pour juger de la pertinence d’un classifieur est la courbe ROC, décrite dans Kuhn (2018)) : on compare la probabilité individuelle (a priori, telle que résultant du modèle de classification) à un seuil, compris entre 0 et 1; si la probabilité est inférieure au seuil, l’estimation est que la personne survit, sinon qu’elle décède.

On compare ensuite cette estimation aux réalisations (ex-post) de survie et de décès. Pour chaque seuil, on peut considérer la matrice classique dite matrice de confusion de théorie de la décision : elle consiste à répartir les observations suivant le résultat observé (en colonne) et l’estimation résultant du modèle en ligne (en fonction de la probabilité estimée pour l’individu et le seuil que l’on s’est fixé). On peut ainsi partager la population entre les classements corrects, et les erreurs (dont les « faux positifs » si la personne a survécu malgré une probabilité estimée de décès supérieure au seuil, et les « faux négatifs » si la personne décède malgré une probabilité estimée inférieure au seuil).

Figure 1: Courbe ROC et classification pour un seuil de probabilité valant 1.5%.

La courbe ROC est obtenue en faisant varier le seuil. Chaque seuil correspond à un point de la courbe, rapportant graphiquement les taux de faux positifs (en abscisse) et de vrais positifs (en ordonnée), comme sur la Figure 1.

Considérons un groupe de 1000 assurés, où 20 personnes sont décédées l’an passé. Supposons un modèle dans lequel on admet que la population est parfaitement homogène, la probabilité estimée de décès est de 2% pour tout le monde. Dans ce cas pour tout seuil supérieur à 2%, on estimera que la totalité de la population survit : on aura un taux de faux positifs de 0% et un taux de vrais positifs de 0%, d’où un point (0,0) sur le graphe. A l’inverse pour tout seuil inférieur à 2%, on estimera que la totalité de la population décède : on aura un taux de faux positifs de 100% et un taux de vrais positifs de 100%, d’où un point (1,1) sur le graphe. La courbe de ROC de ce modèle uniforme à 2% est donc la diagonale du carré sur la figure 1.

Mais on peut aussi imaginer qu’il existe un peu d’hétérogénéité avec, par exemple, une probabilité de décès de 1% pour une moitié de la population et de 3% pour l’autre moitié, ou encore que le modèle produit des probabilités comprises entre 1% et 3% de façon non dichotomique. Les données simulées pour construire la courbe noire sur la Figure 1 suppose que la population a des probabilités de décès variables, comprises entre 1% et 3%, obtenues par une régression logistique. Comme le montre le tableau de droite, on commet des erreurs, et comme le montre celle de gauche, la nature de celle-ci varie en fonction du seuil choisi, qui modifie les taux de faux positifs et de faux négatifs.

Le cas extrême serait celui où le modèle aurait correctement attribué une probabilité de 100% aux 20 personnes qui sont effectivement décédées. C’est la courbe rouge sur la Figure 1. Ce partage est possible ex-post, une fois réalisation de l’aléa : a posteriori, il y a une certitude de décès pour ceux qui sont effectivement morts. Mais cela n’a cependant pas grande réalité dans l’assurance, à moins d’imaginer que l’actuaire serait un oracle, qui saurait avec certitude qui va mourir, et qui va survivre. La réalité est plutôt celle de la situation intermédiaire entre la courbe rouge et la diagonale, avant d’arriver dans la région hachurée, où le taux d’erreur est faible, mais pas nul : on ne peut pas prédire, avec certitude, qui va décéder. L’assurance n’est possible que si cette borne supérieure n’est pas trop élevée.

Incertitude statistique, données et modèles

Une question fondamentale pour la survie de l’assurance est de savoir où se situe cette borne supérieure : jusqu’où peut-on aller, entre les deux cas extrêmes (population homogène avec une probabilité de 2% pour tous, et une population très discriminée, avec 2% de la population ayant 100% de chances de mourir, et l’autre 0%) ? Et de quoi cette borne dépend-elle ? En particulier, des modèles plus complexes, tels que les réseaux de neurones très profonds permettent-ils vraiment d’améliorer la prévision ? Et l’enrichissement de données, tel qu’on l’observe grâce aux objets connectés et la fusion avec toutes sortes d’informations externes, va-t-il déplacer la borne supérieure vers le haut ?

Si l’apprentissage profond – voir Goodfellow et al. (2018) – permet d’avoir des classifieurs d’images avec un taux d’erreur proche de 0%, il est difficile d’imaginer qu’il sera possible de prévoir, presque un an à l’avance (à la signature du contrat), qui décèdera dans l’année, qui aura la grippe, qui aura un dégât des eaux, etc. Les modèles plus complexes permettent d’améliorer les prévisions, en tenant compte de non-linéarités, d’effets croisés entre les variables tarifaires, mais pas au point de faire disparaître l’aléa. Et tant que l’assurance est envisagée ex-ante (la prime est fixée au début de la période de couverture), il est difficile d’imaginer que rajouter de l’information fera aussi disparaître l’aléa. C’est d’ailleurs le cas pour les tests génétiques qui n’expliquent qu’une (petite) partie du risque de cancer, par exemple. Et rajouter des données revient souvent à rajouter du bruit, ce qui rend le travail d’analyse plus complexe. Cependant, force est de constater que des modèles plus complexes et des données plus riches tendent effectivement à « améliorer » la prévision, en remontant la courbe ROC vers le haut. Mais se pose-t-on les bonnes questions ? Que signifie vraiment une borne très éloignée du cas homogène, sur la diagonale ?

Homogénéité, équité et causalité

Comme nous l’avons vu, la tarification en assurance repose sur une répartition des risques (des contrats) en catégories, au sein desquelles la distribution des pertes peut être estimée, afin de fixer un niveau de primes. La répartition se fait à partir des caractéristiques de l’assuré, et du bien assuré. En retraçant l’histoire de l’assurance, Ewald (1986) montre que les mécanismes de prévoyance se sont mis en place en déplaçant la charge des accidents du travail sur la société : on abandonne l’idée d’une responsabilité individuelle de l’accident en faveur de la solidarité. L’assurance distingue « entre le dommage que subit tel ou tel individu — c’est affaire de chance ou de malchance — et la perte liée au dommage dont l’attribution est, quant à elle, toujours collective et sociale ». Ce principe de solidarité sociale, de mutualisation des risques, fait que le risque (en assurance) est toujours pensé collectivement.

Aujourd’hui, les tarifs sont considérés comme « justes », ou « actuariellement équitables » si chaque prime correspond à la perte attendue (pour ne pas dire « espérée », au sens mathématique) pour chaque assuré. Dans cette perception de l’équité, une hypothèse essentielle est que les classes soient « homogènes ». En effet, dans l’hypothèse inverse, les personnes les moins risquées subventionnent les personnes les plus risquées, ce qui est perçu comme socialement injuste.

On peut décrire cette version de l’équité actuarielle à l’aide de la formule de décomposition de la variance. La variance globale se décompose en effet en deux termes, la variance inter-classes et la variance intra-classes : l’ « équité actuarielle » vise à ce que  les classes de risque soient relativement distinctes les unes des autres, donc une variance inter-classes forte, accompagnée d’une homogénéité des classes, donc une variance intra-classes faible. D’un point de vue statistique, chercher à augmenter l’une est équivalente à faire diminuer l’autre. Cette mécanique n’est pas toujours claire pour des observateurs non avertis ; ainsi dans l’affaire Manhart, un des cas les plus documentés sur la discrimination par le genre en assurance, le juge Stevens affirme : « we focus on fairness to individuals rather than on fairness to classes […] even a true generalization about a class is an insufficient reason for disqualifying an individual to whom the generalization does not apply» (cité dans Anzalone (2016)). Autrement dit, pour la justice, un critère statistique de type « true generalization » ne peut être opposable à un individu.

Une autre critique importante, que l’on retrouve dans la « gender directive », est le lien entre discrimination et causalité. En effet, statistiquement, les actuaires vont chercher des facteurs de classification fortement corrélés avec la sinistralité. Mais il est possible que ces facteurs ne soient qu’un proxy de la vraie variable causale, restée elle inobservée, conduisant à une mauvaise estimation du risque pour certains. Comme le notent Antonio et Charpentier (2017), le genre a ainsi été longtemps utilisé en assurance automobile car très corrélé avec des variables associées au style de conduite et à d’autres variables historiquement non observables (mais qui le sont aujourd’hui grâce aux objets connectés, comme le kilométrage, les heures de conduites, les types de routes utilisés, etc).

Ce lien avec les mécanismes causaux est d’ailleurs relativement profond, et Hacking (1975) y voit une connexion avec la « révolution probabiliste » : on peut assez facilement mettre en évidence des corrélations, mais les causes, si elles existent, nous restent plus opaques. Laplace au début du 19e siècle déclare ainsi que « la probabilité est relative en partie à nos connaissances, en partie à notre ignorance », liant les probabilités à la fois à une vision newtonienne déterministe du monde et à notre incapacité à le connaitre parfaitement. Cette dernière composante fait que l’on ne peut pas annoncer la date exacte du décès d’un individu, mais statistiquement, dans un groupe homogène, on peut prédire le nombre de décès au cours d’une année. Et pour revenir à la relation causale, le tabagisme par exemple ne cause pas forcément une mort prématurée mais fumer sera vu comme dangereux car il augmente la probabilité de décès pendant une période donnée. Ainsi nous montre Hacking (1975), la causalité est pensée aujourd’hui dans un contexte probabiliste, et non plus déterministe.

Abraham, K.S. (1985) Efficiency and fairness in insurance risk classification, Virginia Law Review 71: 403-451.

Anderson, A.W. (1978). A Critique of the Manhart Brief. The Actuary, 12:5.

Antonio, K. & Charpentier, A. (2017). La tarification par genre en assurance, corrélation ou causalité ? Risques, 109.

Anzalone, C.A. (2016). U.S. Supreme Court Cases on Gender and Sexual Equality. Routledge.

Bailey, H., Hutchison, T. & Narber, G. (1975) The regulatory challenge to life insurance classification, Drake Law Review Insurance Law Annual 4: 779-827

Barry L. (2019). Justice ou justesse? L’équité de l’assurance. Working paper #15, chaire PARI.

Charpentier, A. & Denuit, M. (2004). Mathématiques de l’Assurance Non-Vie : Principes Généraux de Théorie du Risque. Economica.

Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

Cummins, J.D., Smith, B.D., Vance, R.N. & VanDerhai, J.L. (1982). Risk Classification in Life Insurance. Kluwer-Nijhoff Publishing.

Ewald F. (1986). L’État providence. Grasset.

Fisher, R. A. (1936).  The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics. 7 (2): 179–188.

Frézal, S. & Barry, L. (2019). Fairness in Uncertainty: Some Limits and Misinterpretations of Actuarial Fairness, Journal of Business Ethics.

Goodfellow, U., Bengio, Y. & Courville, A. (2018) L’apprentissage profond. Massot Edition.

Hacking, I. (1975) The Emergence of Probability. Cambridge University Press

Hoy, M. (1982) Categorizing risks in the insurance industry, The Quarterly Journal of Economics 97: 321-336.

Kleindorfer, P. & Kunreuther, H. (1980) Misinformation and Equilibrium in Insurance Markets, in Economic Analysis of Regulated Markets, Jörg Finsinger Editor, Springer Verlag, 67-90

Kuhn, M. & Johnson, K. (2018). Applied Predictive Modeling. Springer Verlag.

Martin,G.D. (1976) Gender Discrimination in Pension Plans, Journal of Risk and Insurance 43.

McGrayne, S.B. (2012) The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy. Yale University Press.

Ramsey P.F. (1926). Truth and Probability.

Schauer, F. (2006) Profiles, Probabilities, and Stereotypes. Harvard University Press.

Von Mises, R. (1957). Probability, Statistics and Truth. Dover publications.

Wilson, C. (1977). A model of insurance markets with incomplete information. Journal of Economic Theory, 16:2, 167-207.

[1] Au sens statistique du mot, dans le sens introduit par Fisher (1936).

[2] Dans cette approche fréquentiste, et notamment pour Ronald Fisher et Richard von Mises, la probabilité d’un évènement unique (dit « one shot ») n’a pas de sens.

The “probability to win” is hard to estimate…

Real-time computation (or estimation) of the “probability to win” is difficult. We’ve seem that in soccer games, in elections… but actually, as a professor, I see that frequently when I grade my students.

Consider a classical multiple choice exam. After each question, imagine that you try to compute the probability that the student will pass. Consider here the case where we have 50 questions. Students pass when they have 25 correct answers, or more. Just for simulations, I will assume that students just flip a coin at each question… I have n students, and 50 questions

set.seed(1)
n=10
M=matrix(sample(0:1,size=n*50,replace=TRUE),50,n)

Let X_{i,j} denote the score of student i at question j. Let S_{i,j} denote the cumulated score, i.e. S_{i,j}=X_{i,1}+\cdots+X_{i,j}. At step j, I can get some sort of prediction of the final score, using \hat{T}_{i,j}=50\times S_{i,j}/j. Here is the code

SM=apply(M,2,cumsum)
NB=SM*50/(1:50)

We can actually plot it

plot(NB[,1],type="s",ylim=c(0,50))
abline(h=25,col="blue")
for(i in 2:n) lines(NB[,i],type="s",col="light blue")
lines(NB[,3],type="s",col="red")


But that’s simply the prediction of the final score, at each step. That’s not the computation of the probability to pass !

Let’s try to see how we can do it… If after j questions, the students has 25 correct answer, the probability should be 1 – i.e. if S_{i,j}\geq 25 – since he cannot fail. Another simple case is the following : if after j questions, the number of points he can get with all correct answers until the end is not sufficient, he will fail. That means if S_{i,j}+(50-i+1)< 25 the probability should be 0. Otherwise, to compute the probability to sucess, it is quite straightforward. It is the probability to obtain at least 25-S_{i,j} correct answers, out of 50-j questions, when the probability of success is actually S_{i,j}/j. We recognize the survival probability of a binomial distribution. The code is then simply

PB=NB*NA
for(i in 1:50){
  for(j in 1:n){
    if(SM[i,j]&gt;=25) PB[i,j]=1
    if(SM[i,j]+(50-i+1)&lt;25)   PB[i,j]=0
    if((SM[i,j]&lt;25)&amp;(SM[i,j]+(50-i+1)&gt;=25)) PB[i,j]=1-pbinom(25-SM[i,j],size=(50-i),prob=SM[i,j]/i)
  }}

So if we plot it, we get

plot(PB[,1],type="s",ylim=c(0,1))
abline(h=25,col="red")
for(i in 2:n) lines(PB[,i],type="s",col="light blue")
lines(PB[,3],type="s",col="red")

which is much more volatile than the previous curves we obtained ! So yes, computing the “probability to win” is a complicated exercice ! Don’t blame those who try to find it hard to do !

Of course, things are slightly different if my students don’t flip a coin… this is what we obtain if half of the students are good (2/3 probability to get a question correct) and half is not good (1/3 chance),

If we look at the probability to pass, we usually do not have to wait until the end (the 50 questions) to know who passed and who failed

PS : I guess a less volatile solution can be obtained with a Bayesian approach… if I find some spare time this week, I will try to code it…

(A brief) history of randomness, and simulation techniques

Hearing “there is a 10% chance of rain today” or “the medical test has a positive predictive value of 75%” shows that the probabilities are now everywhere. A probability is a quantity that is difficult to grasp, but essential when trying to theorize and measure chance, or randomness. And if mathematical theory finally came very late, as Hacking (2006) points out, this did not prevent insurance from developing early enough, and from having the first (actuarial) mortality tables even before the “probability of death” or “life expectancy” had a mathematical basis. And in the same way, many techniques were invented to “generate randomness“, before the explosion of the so-called Monte Carlo methods, in parallel with the development of computing (and the fact that a machine could generate chance). Continue reading (A brief) history of randomness, and simulation techniques

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(
+ "cbad3ca6-6b8f-4c98-9ee0-601faAviationData.txt",
+ 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")

Continue reading The odds of a cluster of airplane accidents

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  random matrice, with  large (say 100) where elements are centered, such as a collection of random variable taking value  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,

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  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  red candies, and  black ones, with . The thing is that no one like those black candies. What could be the probability that for the   kids that will get candies after knocking at my door (with  for convenience, but we will also consider the more general case where I have to many candies, , 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  votes) and B (who receives  votes). A wins the election (). 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,

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),

ballot box lattice diagram

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  and  denote the number of black and red candies, respectively after the th kid git his (or her) candy. Yes, one at a time. Here,  and . What we want is

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)

Here, ties are allowed, and we can prove (easily) that

(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  kids show up. We can see how the probability

will change, with ,

> 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  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 , with some colors inside, say  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  cells, each of them has one color, chosen among . Which gives us  possible states…. And that can be large. I mean, if we consider the smallest possible grid (that might be interesting), say , and only  colors, then we talk about 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  elements. More generally, we talk about writing down matrices with  elements. If we want black and white  grids, that would mean a matrix with  which mean 4 billion elements ! And if we consider an red-green-blue  grid, we have to explicit a matrix with  i.e almost 400 million elements. So, let’s face it: we can only work with  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 https://latex.codecogs.com/gif.latex%20?\{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 https://latex.codecogs.com/gif.latex%20?\{5,6,8\}, then it could be possible to change the color. And actually, https://latex.codecogs.com/gif.latex%20?\{5\} is different from https://latex.codecogs.com/gif.latex%20?\{6,8\} (since it does have much more neighbors)
  • if we chose cell https://latex.codecogs.com/gif.latex%20?\{9\}, then definitively, the color will change, since all neighbors have the other color here,

Now, the probability to select cell  given that there was a color change would be, if  is in https://latex.codecogs.com/gif.latex%20?\{9\}

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)\propto%20\frac{3}{3}=1

while if is in https://latex.codecogs.com/gif.latex%20?\{6,8\}, then there are 4 out 5 neighbors that are red, so

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)\propto%20\frac{1}{5}and if is https://latex.codecogs.com/gif.latex%20?\{5\}, then, only one neighbor has a different color, out of 8, so

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)\propto%20\frac{1}{8}

And for the other, https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)\propto%200. So, it comes – since we assume that cells are drawn independently, and with the same probability, if  is in https://latex.codecogs.com/gif.latex%20?\{9\}

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)=%20\frac{1%20\cdot%20\frac{1}{9}}{\left(1+2\times%20\frac{1}{5}+%20\frac{1}{8}+5\times%200\right)\cdot%20\frac{1}{9}}=\frac{40}{61}

while if is in https://latex.codecogs.com/gif.latex%20?\{6,8\}, then there are 4 out 5 neighbors that are red, so

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)=%20\frac{\frac{1}{5}%20\cdot%20\frac{1}{9}}{\left(1+2\times%20\frac{1}{5}+%20\frac{1}{8}+5\times%200\right)\cdot%20\frac{1}{9}}=\frac{8}{61}

and if is https://latex.codecogs.com/gif.latex%20?\{5\}, then, only one neighbor has a different color, out of 8, so

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(k)=%20\frac{\frac{1}{8}%20\cdot%20\frac{1}{9}}{\left(1+2\times%20\frac{1}{5}+%20\frac{1}{8}+5\times%200\right)\cdot%20\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

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

(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,

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  and the largest  . 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.

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  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  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 -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 -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 https://latex.codecogs.com/gif.latex%20?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 -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 -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  (say $100). The goal is to find the strategy which maximizes the probability to leave Las Vegas with https://latex.codecogs.com/gif.latex%20?2s (here $200). Should I play big, or small ?

Assume that I can bet https://latex.codecogs.com/gif.latex%20?x (that will be, here, for convenience, a fraction of ). With probability https://latex.codecogs.com/gif.latex%20?p, I will get https://latex.codecogs.com/gif.latex%20?2x, and with probability https://latex.codecogs.com/gif.latex%20?1-p, I will get https://latex.codecogs.com/gif.latex%20?0 (and lose my https://latex.codecogs.com/gif.latex%20?x). As mentioned above, https://latex.codecogs.com/gif.latex%20?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 https://latex.codecogs.com/gif.latex%20?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 https://latex.codecogs.com/gif.latex%20?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 , 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 https://latex.codecogs.com/gif.latex%20?x%20=%20s, I play everything I have, and I double with probability https://latex.codecogs.com/gif.latex%20?p. That one was simple.  And indeed, on the graph above, the point on the right is probability  https://latex.codecogs.com/gif.latex%20?p (the red horizontal line).

Assume now that I can bet https://latex.codecogs.com/gif.latex%20?x%20=%20s%20/%202, and I will play, at least, two rounds

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

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

which is

https://latex.codecogs.com/gif.latex%20?p%20^%202%20\left%20(1%20+2%20p%20(1-p)%20+%20[2p%20(1-p)]%20^%202%20+%20\cdots%20\right)%20=%20\frac%20{p%20^%202}%20{1-2p%20(1-p)}

Let’s try something more general: I have initial wealth , I can bet https://latex.codecogs.com/gif.latex%20?x and the goal is to reach  (or, more generally, say, ). Now, the probability to reach  from  betting (always) https://latex.codecogs.com/gif.latex%20?x is exactly the same as the probability to reach  from  betting only 1. Let  denote the probability to go from  to  betting 1 (let us use generic parameters). We can easily get the following equation

https://latex.codecogs.com/gif.latex%20?P_b(a)%20=%20p\cdot%20P_b(a+1)%20+%20(1-p)%20\cdot%20P_b(a-1)

Thus, we can write

https://latex.codecogs.com/gif.latex%20?p\cdot%20(P_b(a+1)-P_b(a))%20=%20(1-p)\cdot%20(P_b(a)-P_b(a-1))

or equivalently

https://latex.codecogs.com/gif.latex%20?(P_b(a+1)-P_b(a))%20=\frac{1-p}{p}\cdot%20(P_b(a)-P_b(a-1))

https://latex.codecogs.com/gif.latex%20?\left(\frac{1-p}{p}\right)^a\cdot%20(P_b(1)-P_b(0))

Now, observe that https://latex.codecogs.com/gif.latex%20?P_b(0)=0 (since I cannot have a gain without any money).

Let us write https://latex.codecogs.com/gif.latex%20?P_b(a+1)-P_b(0) using a domino technique :

https://latex.codecogs.com/gif.latex%20?[P_b(a+1)-P_b(a)]+[P_b(a)-P_b(a-1)]+\cdots+[P_b(1)-P_b(0)]

i.e.

https://latex.codecogs.com/gif.latex%20?\left(\frac{1-p}{p}\right)^a%20P_b(1)+\left(\frac{1-p}{p}\right)^{a-1}%20P_b(1)+\cdots+%20\left(\frac{1-p}{p}\right)^0%20P_b(1)

so this geometric sum can also be written

https://latex.codecogs.com/gif.latex%20?\left(1%20-\left[\frac{1-p}{p}\right]^{a+1}%20\right)%20\left(1%20-\left[\frac{1-p}{p}\right]%20\right)^{-1}

Finally, we can write

https://latex.codecogs.com/gif.latex%20?P_b(a)=\left(1%20-\left[\frac{1-p}{p}\right]^{a}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]%20\right)^{-1}\cdot%20P_b(1)

Here, there is still  that I have to explicit. The idea is to observe that , thus

https://latex.codecogs.com/gif.latex%20?P_b(a)=\left(1%20-\left[\frac{1-p}{p}\right]^{a}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]^{b}%20\right)^{-1}

So finally,

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(gain)=\left(1%20-\left[\frac{1-p}{p}\right]^{s/x}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]^{2s/x}%20\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 ,

> 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 . Which makes sense… Further, the probability is decreasing with : the more hungry, the less chance of winning I have.

Now, the interesting part is what is plotted on the graphs above: the smaller  (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

http://freakonometrics.hypotheses.org/files/2013/02/Capture-d%E2%80%99e%CC%81cran-2013-02-18-a%CC%80-13.27.42.png

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

http://freakonometrics.hypotheses.org/files/2013/02/Capture-d%E2%80%99e%CC%81cran-2013-02-18-a%CC%80-11.34.23.png

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

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

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

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

Welch Two Sample t-test

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

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

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

Welch Two Sample t-test

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

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

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

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

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

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

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

or for average gains

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

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

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

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

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.

http://freakonometrics.hypotheses.org/files/2015/12/extra06.gif

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