Category Archives: Probability

Un double centenaire : Treatise on probabilities de John Maynard Keynes et Risk, Uncertainty and Profit de Frank Knight

Il y a tout juste cent ans, paraissait deux ouvrages fondamentaux essayant de lier l’économie, le risque, l’incertain et le calcul des probabilités. Et si l’occasion est rare de se plonger dans des textes aussi anciens, y compris pour la majorité des étudiantes et des étudiants, il est agréable de noter que ces ouvrages contiennent de nombreuses réflexions passionnantes, qui éclairent d’un œil (presque) nouveau ces domaines où nous avons parfois l’impression de tout savoir, depuis tout ce temps…

Cambridge, Keynes, et les probabilités

Commençons par Treatise on probabilities, le premier ouvrage publié par un des plus grands économistes du début du XXème siècle, et qui témoigne de l’importance des mathématiques dans le formalisme économique – John Maynard Keynes avait déjà publié des articles économiques depuis une dizaine d’année lorsque le traité est publié. Et étonnamment, cet ouvrage est souvent oublié par les économistes[i].

Mais commençons par le début… John Maynard Keynes est né à Cambridge en juin 1883, son père (John Neville Keynes) était professeur de logique et d’économie à l’université de Cambridge, et sa mère (Florence Ada Brown) deviendra maire de la ville (en 1932, une autre femme – Eva Hartree – avait déjà été élue maire dix ans auparavant[ii]). Au début du XXème siècle (en réalité entre 1880 et 1940, comme le montre MacLeod & Urquiola (2020)), Cambridge était probablement la plus prestigieuse université au monde. En 1903, Bertrand Russell y publie son Principles of Mathematics (dont le titre de la traduction en français – Écrits de logique philosophique – reflète peut-être mieux le contenu que le titre original). L’ouvrage de John Maynard Keynes est d’ailleurs dans cette lignée, car malgré un formalisme mathématique indéniable, le traité est peut-être avant tout un ouvrage de philosophie et de logique. C’est aussi l’époque où Srinivasa Ramanujan est sur le campus, à Cambridge, invité du Tamil Nadu, en Inde, pour travailler avec Godfrey Harold Hardy et John Littlewood, comme le raconte The Man Who Knew Infinity, de Robert Kanigel, qui décrit le Cambridge des années 1910. John Maynard Keynes y côtoiera aussi de nombreux artistes (en tant que membre du groupe, ou association, Bloomsbury) comme Virginia Woolf, ou E. M. Forster

En 1905, alors qu’il vient d’obtenir un Bachelor of Arts en mathématiques, il suit les cours d’Alfred Marshall, qui vont lui faire découvrir l’économie politique, et en 1908 (de retour de deux ans passés à Bombay au Bureau des Affaires Indiennes), il travaille sa thèse sur les probabilités. Il ne s’agit pas formellement d’une thèse de doctorat[iii], mais d’une « Fellowship Dissertation » afin de devenir Fellow du King’s College. Son manuscrit est presque fini lorsque la première guerre mondiale éclate, ce qui empêchera sa publication, et il retravaillera le document pendant encore 2 ans après la fin de la guerre, pour la voir publiée en 1921, sous le titre Treatise of Probabilities. Coïncidence, 1921 est aussi l’année ou Franck Ramsey arrive à Cambridge (il deviendra également Fellow de King’s College en 1924 avant de décéder en 1930, à tout juste 26 ans), dont la proximité intellectuelle avec John Maynard Keynes est importante, même s’ils auront beaucoup de divergence, en particulier sur le concept de probabilité. 1921, c’est enfin un peu moins de 10 ans avant la publication de l’axiomatique des probabilités par Kolmogorov, en 1929, qui placera les probabilités au sein de la théorie de la mesure, et qui proposera le cadre formel encore enseigné aujourd’hui.

Pour comprendre l’orientation prise par Keynes dans son livre, il faut avoir en mémoire les liens entre probabilités et logique. A cette époque, à Cambridge, William Ernest Johnson travaille sur son ouvrage Logic, dont les trois tomes[iv] paraitront en 1921, 1922 et 1924. Harold Jeffreys, à qui on devra de nombreux travaux en statistique Bayésienne, était aussi à Cambridge à cette époque, ainsi que Dorothy Wrinch, ce dernier ayant publié On Some Aspects of the Theory of Probability en 1919. Si la recherche était bouillonnante à cette époque, R. B. Braithwaite, qui signe la réédition de l’ouvrage de Keynes en 1972, note que malgré tout, rares étaient les ouvrages de référence en anglais sur le sujet, le précédant ayant été publié plus de 50 ans auparavant par John Venn (qui aura laissé à la postérité les fameux diagrammes) avec son Logic of Chance. Et si Harold Jeffreys a proposé une axiomatique des probabilités plus solide que le modeste chapitre 2 du Traité, son ouvrage ne sortira qu’en 1939.

Pour Keynes, une « probabilité » n’est pas un fait objectif, naturel, comme le laisse entendre l’interprétation fréquentiste.  Pour Keynes, une probabilité est avant tout affaire de logique, quantifiant un degré de confirmation, ou de croyance rationnelle, d’une proposition. Formellement,

« Let our premises consist of any set of propositions h, and let our conclusion consist of any set of propositions a, then, if a knowledge of h justifies a rational belief in a of degree a,we say that there is a probability-relation of degree a between a and h. »

Elle exprime le degré de croyance qu’il est raisonnable d’avoir envers un ensemble de propositions a, à la lumière d’un ensemble de propositions h. Dans cette mesure, donc, la probabilité peut être qualifiée de subjective, car il s’agit d’une croyance.  Mais dans le sens que lui donne Keynes, la probabilité n’est pas subjective puisqu’elle n’est pas soumise à un caprice humain : c’est une relation logique entre deux propositions (ou deux ensembles de propositions) qui se manifestent dans l’esprit des individus. Ce n’est pas une relation entre une déclaration et une réalité, « probability begins and ends with probability ».  Dans le livre de Keynes, une probabilité peut être exprimée par le symbole ‘a/h‘. La connaissance d’une situation dans laquelle ‘a/h = 1′ est certaine. Une situation dans laquelle ‘a/h = 0′ constitue une impossibilité. Dans la majorité des cas, ce chiffre est situé entre 0 et 1. Et, dans la majorité des cas, les nombres n’ont qu’une signification ordinale. Keynes veut dire par là qu’il est souvent impossible de comparer quantitativement deux probabilités et d’affirmer, par exemple, que les chances qu’un événement se produise sont x fois plus élevées que les chances qu’un autre événement se produise.  De plus, il est possible que deux probabilités soient incomparables : « it is not  always possible to say that the degree of our rational belief in one conclusion is either equal to, greater than, or less than the degree of our belief in another ». Et l’on peut penser que, dans la réalité humaine, les probabilités sont dans la plupart des cas non quantifiables et incommensurables, d’autant plus qu’elles sont contraintes par les limites de la raison humaine. Même lorsque les probabilités individuelles sont mesurables numériquement, nous ne pouvons pas aller très loin dans le raisonnement mathématique. Il est d’ailleurs difficile d’éliminer l’intuition et le jugement direct de la considération des probabilités. Cette difficulté restreint la portée des probabilités fréquentistes fondées sur la loi des grands nombres. Même dans le domaine des sciences naturelles, l’intuition et l’analogie jouent un rôle plus important que la manipulation des fréquences statistiques.

En fait, pour Keynes, la plupart des grands penseurs des probabilité, comme Condorcet, Bernoulli, Bentham, Laplace, ou Edgeworth, se sont trompés en pensant que l’on pouvait appliquer les principes issus de l’équiprobabilité aux sciences morales et ainsi quantifier, mesurer et formaliser mathématiquement la réalité sociale.  Keynes rappelle que certains de ces grands penseurs se sont aventuré dans le domaine de l’éthique, et ont a été amené à croire que les degrés de bonté étaient numériquement mesurables, et parfois additifs, alors que la question relève du jugement intuitif.  Keynes dénonce cette « mathematical charlatanry » qui a fait perdre du temps à la statistique théorique, et qui est soutenue par une confiance totale dans l’inférence statistique.  En économie, cette erreur est soutenue par l’illusion de la quantification selon laquelle des choses comme l’utilité peuvent être quantifiées, mesurées, additionnées et soustraites. Cette erreur conduit à une utilisation abusive des statistiques qui, en tant que moyen de décrire quantitativement la réalité, s’est transformée en un instrument de prédiction. Keynes poursuivra cette réflexion par la suite, avec en particulier des éléments de réflexions sur les probabilités dans sa Théorie Générale, avec le Chapitre 12. En particulier, c’est dans ce second ouvrage que l’on trouve la note de bas de page abondamment citée, où il affirme que « très incertain » et « très improbable » sont pour lui équivalent, laissant entendre que les probabilités serraient l’outils approprié pour modéliser les situations incertaines (récemment, Brady (2019) dédie un article à cette fameuse note 1 de la page 148). Pour une discussion plus poussée sur le contenu du Traité, je vais renvoyer à l’article de Christian Robert sur le sujet.

Pour revenir sur l’histoire du Traité, comme nous l’indiquions en introduction, il est passé inaperçu pour bon nombre d’économistes, tout du moins jusque dans les 1970, cinquante ans après sa parution. Lors de la réédition des travaux de Keynes, il a été décidé de publier ce livre non pas dans l’ordre chronologique comme cela avait été fait pour les sept autres monographies, mais comme un huitième titre. Au début des années 1980, les travaux de Robert Skidelsky, suivis par ceux de Bradley Bateman, Anna Carabelli ou Michael O’Donnell, vont inaugurer un ensemble de travaux étudiant les relations entre philosophie et économie chez Keynes ainsi que la continuité de sa vision philosophique. On peut noter cette réflexion de Peter Waley

« The first major effort to construct a theory of imprecise probability was made by Keynes (1921). Keynes aimed to develop an inductive logic, based on a logical interpretation of probability as a “degree of rational belief” […] Since Keynes, a large literature has grown up concerning the mathematics and interpretation of imprecise, epistemic probabilities »

Mais il ne serait pas honnête d’oublier la critique que Franck Ramsey a faite dans Truth and Probability (publié à titre posthume en 1931, mais basée sur une conférence donnée en 1926).  Pour ce dernier, Le calcul des probabilités consiste à établir un ensemble de règles permettant aux degrés de croyance de former un système cohérent. Les relations de probabilité, telles que Keynes les avait décrites, n’existent pas. Ce questionnement sur les probabilités, l’incertain et les croyances vont se retrouver dans un second ouvrage, publié lui aussi en 1921, mais dans un contexte assez différent.

Chicago, Knight et l’incertain

Le second ouvrage dont on célèbre le centenaire de la publication, Risk, Uncertainty and Profit de Frank Knight a été publié de l’autre côté de l’Atlantique, indépendamment de ce qui se passait en Angleterre. En 1871, Chicago est dévasté par un immense incendie, et sa reconstruction se fera grâce à l’école de Chicago, qui attirera les plus grands architectes, entre 1890 et 1910, avec William Le Baron Jenney, Henry Hobson Richardson, ou Frank Lloyd Wright qui est à Chicago pendant une bonne partie de cette période. Chicago atteint le million d’habitants et devient un moment une ville importante artistiquement, même si on retient souvent qu’elle était la capitale du crime organisé, avec Al Capone par exemple. Et si notre histoire ne commence pas à Chicago, c’est dans cette ville de Franck Knight s’imposera, en dirigeant, avec Jacob Viner, le département d’économie, entre 1920 et 1950. Même si c’est Milton Friedman qui rendit célèbre l’école de Chicago, Franck Knight a marqué par ses enseignement la génération rendue célèbre après-guerre, comme le rappelle Van Overtveldt (2007).

Dans les années 1890, et pendant une quinzaine d’années, les économistes américains essayent de répondre à la question de la justice sociale du partage et de la répartition des richesses. On peut retenir l’article de Frederick Hawley, The Risk Theory of Profits parmi les nombreux articles parus dans le Quarterly Journal of Economics, édité à Harvard, où les économistes se répondent les uns aux autres, pendant plusieurs années. Pour Frederick Hawley, la prise de risque était une composante inévitable de la dynamique de la production, et ceux qui prenaient des risques devaient avoir droit à une récompense, connue sous le nom de « profit ». Selon Frederick Hawley, le profit est le prix payé par la société pour assumer le risque. En particulier, il note que cette prime de risque devrait correspondre à une compensation supérieure à la valeur actuarielle, c’est-à-dire une prime sur le risque calculable. Autrement dit, le profit naît du risque non assurable, et voyant le profit comme un résidu, ce qui a donné le nom « residual theory of profit ». S’en suivra une longue polémique avec John Bates Clark, jusqu’à l’article de Francis Edgeworth en 1904, qui offrira une synthèse du sujet.

Maintenant que le contexte est posé, revenons un instant à Frank Knight. Frank Knight est né dans l’Illinois, en novembre 1885. Il était le fils aîné d’une famille de fermiers, et les contraintes de la vie à la ferme l’ont forcé à mettre ses études de côté un certain temps. Il a quitté sa famille alors qu’il avait plus de 20 ans pour commencer des études collégiales (évangéliques) dans le Tennessee, pour obtenir un Bachelor of Science en sciences naturelles, puis une maîtrise en allemand quelques années plus tard. Il part ensuite à Cornell, pour un doctorat en économie (discipline qu’il découvre à partir de 1913), et défendra sa thèse, intitulée « A Theory of Business Profit », qui avait remporté la deuxième place dans la catégorie générale du concours d’essais économiques Hart, Schaffner & Marx en 1916. A cette époque à Cornell, on retrouve Allyn Abbott Young qui y restera avant d’enseigner à Harvard, puis à la London School of Economics, qui présidera le jury de doctorat de Frank Knight, et qui lui offrira de rester comme stagiaire post-doctoral, jusqu’à son départ pour Harvard en 1917. Frank Knight de son côté rejoindra alors John Bates Clark à l’Université de Chicago, où il retravaillera le manuscrit de sa thèse, qui sera publiée en 1921 sous le titre Risk, uncertainty and profit, comme le raconte Emmet (2020) – qui revient longuement sur les changements entre le manuscrit de la thèse, et l’ouvrage publié. Les superviseurs de Knight – John Bates Clark et Allyn Abbott Young – ont attesté à l’éditeur que peu de choses avaient été modifiées, mais comme pour John Maynard Keynes, ces années de maturation ont permis de clarifier certains concepts.

Le point de départ de livre de Frank Knight, ce sont les travaux de John Clark, et en particulier Distribution of Wealth paru en 1899,

« Professor Clark’s work is taken as the basis and starting point of the present study. Some criticism of it will be found necessary, leading to disagreement on some rather essential points, but with its general aim and with many of its conclusions our own are in accord ».

Tout au long du livre, il revient sur des défauts des approches de Hawley et Clark, en notant par exemple que si les agents économiques connaissaient les distributions de probabilités, ils devraient alors se comporter exactement de la même manière. Même si Clark avait introduit une notion de « valeur actuarielle subjective », l’idée centrale à l’époque était que les probabilités (et donc la valeur actuarielle) était parfaitement connue par tous les agents

« It is this “true” uncertainty, and not risk, which forms the basis of a valid theory of profit and accounts for the divergence be-tween actual and theoretical competition ».

La première idée de Frank Knight est bien sûr la notion d’incertitude. L’ouvrage récent de John Kay et Mervyn King, Radical Uncertainty : Decision-Making Beyond the Numbers, revient sur cette idée centrale dans le livre de Knight. La deuxième idée est la théorie des prix, et plus précisément la théorie de la concurrence parfaite. La différence la plus importante, selon Knight, entre la théorie classique des prix et la concurrence imparfaite est la présence dans cette dernière de l’incertitude, qui émerge d’un changement dynamique, et donc imprévisible. C’est cette imprévisibilité dynamique qui a créé les conditions économiques pour l’action entrepreneuriale. Cette approche nécessitaient un jugement, une faculté humaine qui n’était pas nécessaire dans la concurrence parfaite, où tout était calculable, et calculé, parce que tout était « prévisible ». L’incertitude exigeait un jugement, non seulement sur les états du monde, mais aussi sur la capacité des autres à porter des jugements dans des situations incertaines. L’organisation et la coopération est donc apparue, tant au sein des entreprises que dans l’organisation économique sociale, comme un moyen de coordonner les jugements. Franck Knight a ainsi été le premier fait une distinction entre les problèmes de décision lors que les probabilités sont connues, et lorsqu’elles ne le sont pas, ce qui impose de clarifier (comme avait tenté de le faire Keynes) la notion de probabilité.

Dans Risk, Uncertainty and Profit, Knight a distingué trois types de probabilité différents, qu’il a appelées « probabilité a priori », « probabilité statistique » et « estimations ». Le premier type se situe sur le même plan logique que les propositions mathématiques, l’exemple classique étant la probabilité d’obtenir une face lors d’un lancer de dé (la symétrie du cube permet d’affirmer qu’il y a une chance sur six d’obtenir une des faces). La « probabilité statistique » dépend de l’évaluation empirique de la fréquence, dans des expériences répétées. Dans le langage bayésien contemporain, dans le premier cas, on retrouve la distribution de probabilité a priori. alors que dans le deuxième cas, on retrouverait une distribution a posteriori. Dans le troisième cas, les données existantes ne se prêtent pas à l’analyse statistique, lorsque des expériences répétées sont impossibles. Ce dernier cas est celui qui a le plus intéressé Knight, en tant qu’économiste essayant de comprendre et modéliser le monde des affaires, et la nature du profit dans ce contexte. Knight a identifié la confusion entre « le problème de l’estimation intuitive » et « la logique de la probabilité », qu’elle soit a priori ou statistique

Le livre n’a pas connu un grand succès, seulement deux articles publiés avant 1930 l’ayant cité. Mais en 1933, Lionel Robbins a fait en sorte que la London School of Economics (LSE) réimprime le livre afin que ses étudiants puissent y avoir accès. Et quarante ans plus tard, l’American Economic Association a décerné à Knight la médaille Walker, décernée tous les cinq ans pour honorer l’économiste américain qui a apporté la plus grande contribution à l’économie tout au long de sa carrière.

Arrous, J. (1982). Keynes et les probabilités: Un aspect du ” fondamentalisme ” keynésien. Revue économique, 33(5), 839-861

Bouvier-Patron, P. (1996). La question de la réduction d’incertitude chez F. Knight. L’Actualité économique, 72 (4), 397–415.

Brady, Michael Emmett, How Should the Post Keynesian School Define ‘Uncertainty’? 1 The Only Correct Answer Is to Use Keynes’s Own Definition Given in Footnote 1 on Page 148 of Chapter 12 of the General Theory: Uncertainty Is an Inverse Function of the Weight of the Argument (2016). SSRN:3438090

Combenale, P. (2010). Introduction à Keynes. La Découverte.

Diebold, F.X. Doherty, N.A. & Herring, R.J. (2010). The Known, the Unknown, and the Unknowable in Financial Risk Management: Measurement and Theory Advancing Practice. Pinceton University Press.

Emmet, R.B. (2020). The Writing and Reception of Risk, Uncertainty and Profit. SSRN: 3591596.

Kent, S. (2019). ‘A degree of a different character’: 100 years of the Cambridge PhD. Cambridge University Library Special Collections Blog

Leroy, S.F. & Singell, L.D. (1987). Knight on risk and uncertainty. Journal of Political Economy, XCV, 2, p. 294-406.

MacLeod, W.B. & Urquiola, M. (2021). Why Does the United States Have the Best Research Universities? Incentives, Resources, and Virtuous Circles. Journal of Economic Perspectives, 35:1, 185-206.

Muchlinski, E. (2003). Épistémologie et probabilité chez Keynes. L’Actualité économique, 79 (1-2), 57–70.

Pradier, PC. & Serrano, D.T. (2000) Frank H. Knight le risque comme critique de l’économie politique. Revue de synthèse 121, 79–116

Robert, C.P (2010). An attempt at reading Keynes’ Treatise on Probability. ArXiv: 1003.5544

Schmidt, C. (2003). Que reste-t-il du Treatise on Probability de Keynes? L’Actualité économique, 79 (1-2), 37–55.

Van Overtveldt, J. (2007). The Chicago school : how the university of Chicago assembled the thinkers who revolutionized economics and business, Chicago, Agate

Walley, Peter ( 1991), Statistical Reasoning with Imprecise Probabilities, London : Chapman and Hall.

Wrinch, D. & Jeffreys H. (1919) On Some Aspects of the Theory of Probability, Philosophical Magazine, no 38,‎ 1919, p. 715-731

Zabell, S.L. (1992) Predicting the unpredictable. Synthese 90, 205–232

[i] Cet ouvrage n’est même pas cité dans l’introduction à Keynes de Pascal Combemale, Combenale (2010)

[ii] En France, la première femme élue à un conseil municipal était Joséphine Pencalet, à Douarnenez en mai 1925 (suite à la grève des sardinières), mais son élection fut invalidée par le Conseil d’État, car à l’époque, les femmes ne sont ni électrices ni éligibles. Il faudra attendre 1945, et le droit de vote, pour que les femmes deviennent éligibles, en Frnace.

[iii] Le programme de Doctor of Philosophy de Cambridge n’existera qu’à partir de 1919, comme le raconte Kent (2019)

[iv] Zabell (1992) revient sur les apports de William Ernest Johnson, qui a en particulier introduit un des premiers relâchements de l’hypothèse d’indépendance, en probabilités, en définissant le concept d’échangeabilité.

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

Bonferroni, les anniversaires et les terroristes (ou pas)

Le paradoxe des anniversaires est un exercice de probabilité que l’on découvre assez vite. Mais assez souvent, on l’oublie, alors qu’il est en réalité assez profond. Le problème est le suivant: quel est le nombre de personnes qu’il faut réunir pour avoir (au moins) une chance sur 2 pour que deux personnes aient leur anniversaire le même jour. Si on calcule la probabilité que deux personnes dans un groupe de n fêtent leur anniversaire le même jour, on peut calculer la probabilité que personne ne soit né le même jour, parmi n personnes. Un calcul rapide montre qu’avec 23 personnes, la probabilité qu’au moins deux personnes fête son anniversaire le même jour est 50.73%. Autrement dit, il suffit d’une groupe de 23 personnes (c’était la taille d’une classe d’élèves, du moins à mon époque…).

L’intuition, pour comprendre pourquoi nous avons l’impression que ce nombre est petit est que, naturellement, si on a un événement qui survient avec probabilité p, et qu’on répète l’expérience n fois, la probabilité que rien ne se produise est (1-p)n. Aussi, en première approximation (et j’insiste sur ce point parce que c’est avec ce genre d’argument fallacieux qu’on conclue que la probabilité d’avoir une catastrophe nucléaire dépasse 100%, conduisant à du “certitude statistique“), on en déduit que la probabilité qu’un événement survienne est de l’ordre de np. Sauf qu’il, on compte des paires. Et avec n personnes, on constitue environ n2/2 paires. Et donc la probabilité devient proche de n2p/2… ce qui est bien plus grand que np. C’est grossièrement ce qui explique que notre intuition se sent trompée : avec seulement 23 personnes, on peut constituer plus de 250 paires. La probabilité que toutes les paires correspondent à des dates d’anniversaire différente n’est pas si petit qu’il y parait.

Jeff Ullman parle ainsi de principe de Bonferroni,

“if you look in more places for interesting patterns than your amount of data will support, you are bound to find crap”

Et cette intuition faussée quand on travaille sur des paires se retrouve dans les techniques utilisées par les services de renseigngement, tel que le montre Anand Rajaraman dans Mining Massive DatasetsSupposons qu’il existe un milliard d’individus qui pourraient potentiellement être des terroristes. Que les gens, quels qu’il soient vont à l’hôtel 1 jour sur 100. Un hôtel donné peut accueillir 100 personnes. On suppose qu’il existe 100,000 personnes, de sorte qu’en moyenne ils soient plein tous les jours. Et supposons qu’on regarde sur un historique de 3 ans, i.e. 1000 jours pour avoir des calculs ronds. Supposons que pour traquer les terroristes, on cherche les gens qui ont séjourné, ensemble, dans un hôtel, le même jour, au moins deux fois… On imagine que c’est effectivement rare, et donc suspicieux…

La probabilité que deux personnes données dorment à l’hôtel le même jour est faible, 1/100×1/100, soit 0,0001. Et comme il y a 100,000 hôtel, la probabilité que les deux personnes descendent, en plus, dans le même hôtel, il faut diviser 0,0001 par 100,000. Ce qui donne 10-9 . Mais ça, c’est la probabilité pour qu’ils se croisent une fois ! Si on regarde la probabilité qu’ils se croisent deux fois, en trois ans, on arrive à une probabilité de 10-18 Pour ceux qui ont du mal, ça doit faire environ

0,00000000000000001 %

Et encore. Ici il peut en plus se croiser la seconde fois dans un autre hôtel. La probabilité qu’il se croisent dans le même hôtel que la première fois devient encore plus petite !

Supposons maintenant qu’un organisme de surveillance avisé décider de lister les paires de personnes qui, sur des paires de jours, se sont retrouvé dans le même hôtel. Ça doit être suffisamment rare pour être suspicieux. Rappelons que quand on considère n individus, le nombre de paires possibles est de l’ordre de n2/2. Autrement dit, avec un million de personnes, le nombre de paires est de l’ordre de 5×1017. Et pareil, le nombre de paires de jours est de l’ordre de 5×105. Aussi, le nombre d’événements suspicieux que l’on observe, sur ces trois années est de l’ordre de

(5×1017 ) x (5×10) x10-18 = 250000

Autrement dit, 1/4 de millions de paires de personnes sont des terroristes potentiels ! Soit 1/2 millions d’individus.

Maintenant, supposons qu’il y ait 10 terroristes. S’ils sont assez stupides pour se retrouver dans les mêmes hôtels, les mêmes jours, ça veut dire que sur les 500000 personnes qui vont être surveillées, 10 sont effectivement des terroristes, et on a 499990 faux positifs… Enfin, “faux positifs”… ils ont forcément quelque chose d’autre à se reprocher, non ?

Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that  has a similar distribution… so there is no reason to test for a Gamma model for  before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if  has a joint Gaussien distribution, then both marginals are Gaussian, but also . So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for . The econometric interpretation is that with a standard Gaussian linear model, if is normally distributed, not only the conditional distribution  is Gaussian but also the non-conditional distribution of .

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of  is also Gaussian

> library(nortest)
> ad.test(Y)

	Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

	Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here  is continuous. What if we consider a finite mixture here, i.e. takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for . But in practice, would we really reject the Gaussian assumption, for ? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value}
 
> plot(2:100,Vectorize(pv)(2:100),type="l")
> abline(h=.05,col="red")

So here, it could be possible to have also a Gaussian distribution, for . As least to accept that assumption, statistically.

In the context of a Poisson regression, it is well know that it’s not possible to have at the same time  that is Poisson distributed (that’s a Poisson regression) and also  that is Poisson distributed. That simply comes from the fact that

while

and because of the conditional Poisson distribution, then

Thus,

So  cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for .

More generally, it is very difficult to have a distribution family for   that is also the distribution of the non-conditional variable . In the context of a finite mixture ( takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference).

Hence, as a keep saying, before running a regression model on with some given family, it is never a good idea to check if the non-conditional distribution  has the same distribution. Because there is no reason, usually, to remain in the same family.

Collectionner les Figurines Lego avec ses Ami(e)s

Dans mon précédant billet, je regardais combien de figurines Lego je devais acheter pour avoir toute la série, avec ma fille. Pour la série de 16, il fallait, en moyenne, acheter 54 sachets, à 3 euros pièces. Dans les commentaires, @Guillaume et @Gaelle me faisaient noter que les échanges, ça pouvait aider.

Bon, tout d’abord j’aime pas qu’on donne des conseils à mes filles pour avoir raison contre moi ! Parce que oui, c’est sûr qu’en s’y mettant à plusieurs, le coût individuel sera plus faible… Ensuite, parce que les maths me semblent plus compliqués…. Cela dit, en simulant des parties, ça tombe assez vite…

Pour commencer, je fais une petite fonction (histoire d’alléger mes codes) pour incrémenter mes comptages de figurines, à chaque fois que j’en achète une

> f=function(mat,i){
+   mat[i]=mat[i]+1
+   return(mat)
+ }

Ensuite, on peut avoir une petite fonction, correspondant à une simulation, qui me renvoie le nombre de figurine de chaque type que j’ai du acheter pour avoir les 16 une fois. Ou plus. En fait, je veux avoir en tout k fois toutes les figurines (ce qui permet de partager parmi k ami(e)s)

> tirage=function(k=1){
+ test=FALSE
+ cpte=0
+ M=rep(0,16)
+ while(!test){
+   no=sample(1:16,size=1)
+   M=f(M,no)
+   cpte=cpte+1
+   if(min(M)==k) test=TRUE
+ }
+ return(list(compte=M,
+             nbre=cpte))
+ }

Maintenant on peut faire des simulations. Par exemple, pour retrouver le fait qu’il faut, en moyenne, acheter 54 sachet pour que ma fille (seule) ait toutes les figurines

> N=unlist(lapply(1:10000,function(x) tirage(1)$nbre))
> mean(N)
[1] 54.1854

Si maintenant on le fait avec mes deux filles,

> N=unlist(lapply(1:10000,function(x) tirage(2)$nbre))
> mean(N)
[1] 82.7961

ce qui est (beaucoup) moins que le double. Si au lieux de faire ça entre mes deux filles, mais que ma fille collectionne les figurines avec une amie à elle (et que les papas partagent les frais), on a en moyenne 41 sachets à acheter, chacun,

> mean(N)/2
[1] 41.39805

et pour trois, on passe à 36 sachets

> N=unlist(lapply(1:10000,function(x) tirage(3)$nbre))
> mean(N)/3
[1] 36.04363

etc

> V=unlist(lapply(1:10,function(k) mean(unlist(lapply(1:10000,function(x) tirage(k)$nbre)))))
> plot(1:10,V/1:10,type="b",pch=19)

On peut aller plus loin, en regardant les quantiles. Avec 4 ami(e)s, si on achète tous 45 sachets, on a 90% de chances que tout le monde ait les 16 figurines

D’ailleurs, si on regarde l’évolution du quantile (à 90%) sur la moyenne, on voit qu’il baisse vite avec la taille…. Avec 2 ami(e)s, il faut 35% de sachets en plus pour s’assurer d’avoir toutes les figurines dans 90% de cas (par rapport au nombre moyen) alors qu’il en faut juste 20% de plus avec 8 ami(e)s.

Maintenant, si quelqu’un a le courage (et le temps) de faire les maths, je suis preneur !

Collectionner les Figurines Lego

 

Tous ceux qui ont des enfants – ou qui jouent encore au Lego – le savent : régulièrement, Lego sort des séries de figurines, vendus en sachets individuels hermétiques, un peu comme les anciennes cartes à jouer. On est à la saison 15, et il y a 16 figurines par séries,

D’ordinaire, les enfants réclament un paquet, quand on va dans un magasin de jeux (genre pour acheter un cadeau pour un anniversaire) mais là, ma fille les a réclamé tous.

Bon, après le moment où je ris dans ma barbe, j’ai quand même fini par me demander combien ça pourrait coûter si on achetait des paquets, de manière à récupérer les 16 figurines.

La bonne stratégie pour faire le calcul me semble d’ordonner les personnage par ordre d’apparition,  et de modéliser le nombre de paquets achetés entre deux premières apparitions, comme indiqué sur le dessin ci-dessous

Si  désigne le nombre total de paquets à acheter, on note  le nombre de paquets à acheter, entre l’apparition de la ()ème nouvelle figurine, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première figurine), on a

On notera que, si  désigne la loi géométrique,
et à partir de là, les calculs sont simples, puisque

soit (par linéarité de l’espérance)

i.e.

ou encore

Le nombre moyen de paquets à acheter est alors
> sum(16/(1:16)) 
[1] 54.09166

Pour ceux qui ont oublié leur calculatrice, on peut se souvenir que

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)
> log(16)*16 
[1] 44.36142

qui diffère du résultat numérique car il manquait la constante d’Euler

https://latex.codecogs.com/gif.latex?gamma%20=%20lim_{n%20<br /> ightarrow%20infty%20}%20left(%201+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20-%20ln(n)%20<br /> ight)

i.e.

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)+gamma

Numériquement, on obtient, en prenant pour https://latex.codecogs.com/gif.latex?gamma 0.57721 (cf Google)

> (log(16)+ 0.57721)*16 
[1] 53.59678

Mais je m’écarte du sujet… à 3€, il faut débourser, en moyenne, 162€ pour avoir les 16 personnages, soit plus de 10€ (en moyenne) par personnage !

> sum(16/(1:16))*3 
[1] 162.275

On est ici au même prix que l’hotel Lego friends (oui, c’est une demande récurrente)

Par simulation, on peut confirmer cet ordre de grandeur (avec ici un million de simulation de la variable),

> rpaquet=function(i) sum(1+rgeom(16,prob=(1:16)/16))
> S=unlist(lapply(1:1e6,rpaquet))
> mean(3*S)
[1] 162.1709

mais surtout, on peut avoir des quantiles, par exemple

> median(S)*3
[1] 150

après avoir déboursé 150 euros, il y a 50% de chances que l’on ait déjà toutes les figurines.

> quantile(S,.9)*3
90% 
237

 

Et il y a 1 chance sur 10 pour qu’on dépasse les 237€. Mine de rien, pour 16 petites figurines, on dépasse là le prix du Falcon Millenium,

Je pense que les enfants ne sont pas prêt que je leur rachète des petites figurines !

Petit exercice de proba

Mardi dernier, on avait fait une série d’exercices de proba, et je voulais reprendre un exercice pour lequel j’avais proposé une vérification sur ordinateur.

Pour résoudre l’exercice, j’avais suggéré la méthode suivante. est tiré parmi , alors que est tiré parmi .  On note assez facilement que le plus petit nombre est , et comme le ppcm de  et  est ,

(je laisse faire les calculs pour montrer qu’effectivement, au delà, on n’est plus dans les ensembles de départ). Bref, seuls  nombres figurent dans les deux ensembles. La probabilité que l’on cherche est alors donnée par

où on somme sur nos  nombres. On a alors

où chacune des probabilités vaut  (on a une chance sur  de tomber sur n’importe quel nombre), et on a une somme de  termes. Donc la probabilité que l’on cherche est

Pour le vérifier, on peut utiliser le petit code (R) suivant

> list_X=seq(2,by=3,length=100)
> list_Y=seq(3,by=4,length=100)
> n=1e8
> x=sample(list_X,size=n,replace=TRUE)
> y=sample(list_Y,size=n,replace=TRUE)
> sum(x==y)/n
[1] 0.00250221

qui confirme le petit calcul que l’on vient de faire.

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

Les anniversaires de vos amis sur Facebook

J’ai découvert avec un peu de retard le joli billet Les anniversaires de vos amis sur Facebook, qui tentait de répondre à la question

Si je possède n amis, quelle est la probabilité qu’il y ait au moins un jour dans l’année où je n’ai pas d’anniversaire à souhaiter ?

Ce problème, on peut aussi l’analyser en posant plutôt la question suivante,

Il y a en tout 365 jours dans une année. Combien d’amis doit-on espérer avoir pour qu’il n’y ait pas un jour sans fêter un anniversaire d’un de ses amis ?

Le problème est différent. Mais la résolution est assez jolie, si on regarde attentivement.

Si je note  le nombre d’amis qu’il me faut pour couvrir tous les anniversaires, et que je commence à regarder, ami après ami le nombre de jours couverts. Supposons que j’ai  jours de l’année couverts, et notons  le nombre de nouveaux amis qu’il va me falloir avant de couvrir  jours de l’année. On peut voir que

Le point intéressant avec cette notation, c’est qu’on connaît la loi de . En fait,  on reconnaît une loi géométrique

Aussi, on peut assez facilement avoir le nombre moyen d’amis qu’il va falloir pour couvrir tous les jours de l’année,

Je laisse les amateurs de calculs essayer de simplifier, mais ça me semble compliqué. Par contre, numériquement on obtient facilement

> sum(365/(1:365))
[1] 2364.646

En fait, la somme, on peut quand même la simplifier un peu, car on reconnaît

Et on sait que

où  est la constante d’Euler. Et effectivement, avec cette approximation, on obtient

> (log(365)+ 0.57721)*365
[1] 2364.144

Ce nombre d’amis est à rapprocher du calcul fait sur https://blogdemaths.wordpress.com/, qui obtenait

on a une chance sur deux pour qu’il y ait au moins un jour sans anniversaire à souhaiter, il faut avoir environ 2285 amis

Calculer un quantile est un peu plus compliqué pour moi (je ne peux calculer que des moyennes), mais on voit que les ordres de grandeurs sont assez proches.

On Hoeffding’s identity

In 1940, Wassily Hoeffding published Masstabinvariante Korrelationstheorie, which was an impressive paper. For those (like me) who unfortunately barely speak German, an English translation could be found in The Collected Works of Wassily Hoeffding, published a few years ago. As I keep saying in my courses about copulas, almost everything was in that paper, by Wassily Hoeffding. For instance, we can see the following graph, of a cumulative distribution function,

What is the difference with a copula? A copula (in dimension 2) is the cumulative distribution function of a random pair with uniform on , as defined by Abe Sklar

But Wassily Hoeffding considered a random pair with uniform on . But everything else is the same. He can even derive the level curves of the density of the Gaussian copula,

> library(mnormt)
> r=.6
> dc=function(u,v) return(
+ as.numeric(dmnorm(cbind(qnorm(u),qnorm(v)),varcov=
+ matrix(c(1,r,r,1),2,2))/dnorm(qnorm(u))/dnorm(qnorm(v))))
> n=500
> vectu=seq(1/n,1-1/n,length=n-1)
> matdc=outer(vectu,vectu,dc)
> contour(vectu,vectu,matdc,levels=
+ c(.325,.944,1.212,1.250,1.290,1.656,3.85),lwd=2)

 

But another interesting point is that there is the so-called Hoeffding’s equality

which is interesting, and quite important, actually, to understand that the covariance (or the correlation) can be seen as some ‘distance‘ to the independence. More precisely, observe that

where  would be the joint cumulative distribution function of some independent variables, with the same marginal distributions.

Of course, it is not exactly a distance, since it can be negative. But still. Now, the thing is that the proof is not trivial. But it is using interesting identities. For instance, in 1885, Franklin wrote a nice paper, Proof of a Theorem of Tchebycheff’s on Definite Integrals, in the American Journal of Mathematics. To get some heuristics about the identity, consider some (finite) sequences  and , then one can prove that

And there is a continuous version of that identity. Consider two bounded functions  and , on some interval,  then

is equal to

In 1979, in Monotone Regression and Covariance Structure, Gerald Shea gave a more probabilistic interpretation of that results, using the fact that

and using a different measure. More precisely, assume now that  functions  and  are integrable, with respect to some measure , on some set . Then

is equal to

In the case where  is a probability measure of , i.e. , this equality is the one used by Wassily Hoeffding, in 1940. The interpretation in terms of random variable is simple that

(with standard assuptions of existence of those quantitites) where  and  are two independent vectors, with identical distribution, . Actually, this relationship can also be found in Some Concepts of Dependence, by E. L. Lehmann, published in 1966. Oh, and by the way, the connection with Chebyshev inequality (claimed in the title of seminal paper by Franklin) come from the fact that if  and  are monotonic, then the left part of the identity is positive, and thus,

But let’s get back to Hoeffiding’s result. How do we get it from that lemma. The idea is to write

as

i.e.

We can then intervert the integral and the expectation, use the fact that

and then, and some integral calculus can be used to rewrite that expression as

So we get here Hoeffding’s identity. Actually, as mentioned by Ben Derrett about the equality above, it can be observed (see http://math.stackexchange.com/105713) that2\text{cov}(X,Y)=2\big(\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y]\big)can also be written

where again,  and  are two independent vectors, with identical distribution, . The later can be writen

Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are  distributions.

  • The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins , and with correlation . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

That was the simple part.

  • The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair , we would like to consider the pair , so that the marginal distributions are . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

  • The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where  and  with independent random variables, given , and that both variables are exponentially distributed, with parameter . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value))
Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value
f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y))
vx=seq(-3,3,length=151)
vy=seq(-3,3,length=151)
z=matrix(NA,length(vx),length(vy))
for(i in 1:length(vx)){
for(j in 1:length(vy)){
z[i,j]=f(vx[i],vy[j])}}
set.seed(1)
Theta=rgamma(1500,a,b)
Z=cbind(rexp(1500,Theta),rexp(1500,Theta))
X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2]))

There is a small technical problem, but no big deal.

Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian.

37% chance

I don’t know if you ever realized, before, but it is quite common to have 37% chance that something happened (or actually “not happened” if we want to be more rigorous). For instance, consider a  grid, and draw  points randomly (and uniformely). Then, around 37% cells are empty. Or if you consider a cell, on that grid, there is 37% chance, that the cell is empty. You can look, on the animation below,

Actually, it is quite simple to prove this result. This come from the fact that

And  is a common probability. To be more specific, we can write it

which is  when  is a Poisson distribution with parameter  (or with mean ). In our case,  is the number of points in a given cell, and since there are   points for   cells (it is a  grid), then . And the Poisson distribution arise here as a limit of a binomial distribution. More precisely, here . We have here the law of small number. Nice isn’t it?

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!