# Articles for the Probability and Statistics Project

Here are some articles for the project for the graduate crash course on probability and statistics,

For those willing to work on datasets, consider success per school in a national exam, over time (here is the file)
 base=read.table("http://freakonometrics.free.fr/brevet_rennes.csv", header=TRUE,sep=";",dec=",")

Some datasets will be uploaded soon

# Probability and Statistics

Next week, I will start a short course on probability and statistics. The slides of the course are now online. There will be more information soon about the exam and the projects.

# The odds of a cluster of airplane accidents

Recently, there have been a lot of airplane accidents.

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

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

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

for instance the exact location of the crashes,

> library(maps)
> map("world", interior = FALSE)
> points(planes$Longitude,planes$Latitude,
+ pch=19,cex=planes$Total.Fatal.Injuries/50, + col="red") # Generating functions Today, I wanted to publish a post on generating functions, based on discussions I had with Jean-Francois while having our coffee after lunch a couple of times already. The other reason is that I publish my post while my student just finished their Probability exam (and there were a few questions on generating functions). • A short introduction (back on a specific exercise) In the Probability exam, I included an exercise we’ve seen in class, last week. The question is the following (question 16 in the form – in French). Let $F(x)=1-e^{-x}/3$ for $x\geq 0$ and $F(x)=0$ for $x<0$ be the cumulative distribution function of some random variable $X$, i.e. $F(x)=\mathbb{P}(X\leq x)$. What is the moment generating function of $X$, i.e. $M(t)=\mathbb{E}(e^{tX})$ ? Consider some $t\in\mathbb{R}$ (we’ll see later on if some additional constraint are necessary). The tricky part of this exercice appears extremely fast, actually: how could you write $\mathbb{E}(g(X))$ ? I mean, in any probability textbook, the standard answer is • if $X$ is discrete, $\mathbb{E}(g(X))=\sum_x g(x)\cdot \mathbb{P}(X=x)$ • if $X$ is (absolutely) continuous, $\mathbb{E}(g(X))=\int g(x)\cdot f(x)dx$ where $f(\cdot)$ is the density of $X$. Here, $X$ is clearly not a discrete variable. But is it (absolutely) continuous. My (strong) belief is that you need to plot that distribution function to see how it looks like, $x\mapsto F(X)$, for all $x\in\mathbb{R}$ (following recent discussions with Philippe Reka, I will try to post more hand-made graphs) Ooops. It looks like we have a discontinuity in 0. So we have to be a bit carefull here : $X$ is neither continuous nor discrete. Let us use the double projection formula, $\mathbb{E}(Y)=\mathbb{E}(\mathbb{E}(Y\vert Z))$ which can also be writen, if $Z\in\{A,B\}$, $\mathbb{E}(Y)=\mathbb{P}(Z=A)\cdot \mathbb{E}(Y\vert Z=A)+\mathbb{P}(Z=B)\cdot \mathbb{E}(Y\vert Z=B)$ This is simply the idea of saying that the overall average is a barycenter of the average per subgroup. Here,$Y=g(X)$ and let $A=\{X=0\}$ while $B=\{X>0\}$ (note that $\mathbb{P}(X\in \{A\cup B\})=1$). Thus, $\mathbb{E}(g(X))=\mathbb{P}(X=0)\cdot \mathbb{E}(g(X)\vert X=0)+\mathbb{P}(X>0)\cdot \mathbb{E}(g(X)\vert X>0)$ Let us consider the three different components. $\mathbb{P}(X=0)=F(0)=1-1/3=2/3$ $\mathbb{P}(X>0)=1-\mathbb{P}(X=0)=1/3$ and $\mathbb{E}(g(X)\vert X=0)=g(0)$ (since it is is a real-valued constant), and here $g(0)=e^{0}=1$. So finally, we should compute $\mathbb{E}(g(X)\vert X>0)$. Observe that $X$ given $X>0$ is a (absolutely) continuous random variable, with a density. To get it, observe that for all $x>0$, $\overline{F}_\star(x)=\mathbb{P}(X>x\vert X>0)=\frac{\mathbb{P}(X>x)}{\mathbb{P}( X>0)}=\frac{e^{-x}/3}{e^{-0}/3}=e^{-x}$ and $f_\star(x)=e^{-x}$, i.e. $X$ given $X>0$ is an exponential distribution. Hence, $X$ is a mixture between an exponential variable and a Dirac mass in $0$. This was actually the tricky part of the question since it is not obvious when we see (only) the formula above. From now on, it is just high-school level computations, $\mathbb{E}(g(X)\vert X>0)=\int_0^\infty g(x) f_\star (x)dx=\int_0^\infty e^{(t-1)\ x} dx=\frac{1}{1-t}$ if $t\leq 1$ (for the first time, we see that the function is not defined everywhere). If we put all the expressions together, $M(t)=\frac{2}{3}\cdot 1 + \frac{1}{3}\cdot \frac{1}{1-t}=\frac{3-2t}{3-3t}$ • Monte Carlo computations If we are lazy (and trust me, I am extremely lazy), it is possible to use Monte Carlo simulations to compute that function, > F=function(x) ifelse(x<0,0,1-exp(-x)/3) > Finv=function(u) uniroot(function(x) F(x)-u,c(-1e-9,1e4))$root

or (to avoid the problem of the discontinuity)

> Finv=function(u) ifelse(3*u>1,0,uniroot(function(x)
+ F(x)-u,c(-1e-9,1e4))$root)) Here, the inverse is simple to get, so we can faster the code using > Finv=function(u) ifelse(3*u>1,0,-log(3*u)) Then, we use > rF=function(n) Vectorize(Finv)(runif(n)) > M=function(t,n=10000) mean(exp(t*rF(n))) > Mtheo=function(t) (3-2*t)/(3-3*t) > u=seq(-2,1 ,by=.1) > v=Vectorize(M)(u) > plot(u,v,type="b",col='blue') > lines(u,Mtheo(u),col="red") The problem with Monte Carlo simulations is that they should be used only if they are valid. If mean, I can compute > set.seed(1) > M(3) [1] 5748134 Finite sum can always be computed, numerically. Even if here, $\mathbb{E}(e^{3X})$ does not exist (or to be more precise, is not finite). It is like the average of a Cauhy sample… I can always compute it, even if the expected value does not exists… > set.seed(1) > mean(rcauchy(1000000)) [1] 0.006069028 This is related to questions I tried to ask a few years ago in a paper, where I wanted to test if $X\in L_1$ (or not). Almost all the tests I know are actually based on that assumption… But this is not the point here. My point is that those generating functions are interesting, when then exist. And perhaps working with characteristic function is a better idea. • Generating functions Now, to get back on the begining of last course, generating functions are interesting for a lot of reasons. But first of all, let us define those function properly. The moment generating function $M_X(t)=\mathbb{E}(e^{tX})$ exists if it is finite on a neighbourhood of $0$ (there is an $a>0$ such that for all $t\in[-a,+a]$, $M_X(t)<\infty$). In that case, there exists some (open) interval $(a,b)\in\overline{R}$ such that for all $t\in(a,b)$, $M_X(t)<\infty$, called the convergence strip of the moment generating function. This function is said to be moment generating, since if $M_X(\cdot)$ exists (as defined in the previous paragraph), then all moments exist, for all $k\in \mathbb{N}\backslash\{0\}$, $\mathbb{E}\left(\vert X\vert^k\right)<\infty$. This is basically due to the fact that, for all $k\in \mathbb{N}\backslash\{0\}$$x^k\exp(-\vert t\vert x)\rightarrow 0$ as $x\rightarrow\infty$, so, for all $x$ large enough, $x^k \leq \exp(\vert t\vert x)$. And before, it is always possible to use a multiplicative constant, $\left\{\begin{array}{l}\text{if }x>x_0, \ \vert x\vert^k \leq \exp(\vert t x\vert)\\\text{if }x\leq x_0, \ \vert x\vert^k \leq K \cdot \exp(\vert t x\vert)\end{array}\right.$ for some $K$. Thus, $\mathbb{E}\left\vert X\vert^k\right) \leq \mathbb{E}\left( e^{\vert tX \vert}\right)\leq \mathbb{E}\left( e^{- tX}\right)+\mathbb{E}\left( e^{ tX}\right)\infty$ if $t$ is small enough (namely $[-t,+t]$ belongs to the convergence strip). Now, if we use Taylor’s expansion, $M_X(t)=\mathbb{E}\left( e^{ tX}\right)=\mathbb{E}\left(\sum_{k=0}^\infty \frac{(tX)^k}{k!}\right)=\sum_{k=0}^\infty\frac{t^k}{k!}\mathbb{E}[X^k]$ and $\frac{\partial^k M_X(t)}{\partial t^k} =\mathbb{E}\left( X^k e^{tX} \right)$ If we look at the value of the derivative of that function at point 0, then $\left. \frac{\partial^k M_X(t)}{\partial t^k}\right\vert_{0} =\mathbb{E}\left( X^k \right)$ As we’ve seen last week in class, it is possible to define a moment generating function in higher dimension, for some random vector $\boldsymbol{X}=(X_1,\cdots,X_d)$, $M_{\boldsymbol{X}}(\boldsymbol{t})=\mathbb{E}\left(e^{\boldsymbol{t}^{{\sffamily T}}\boldsymbol{X}}\right)=\mathbb{E}\left(e^{t_1X_1+\cdots+t_dX_d}\right)$ for some $\boldsymbol{t}\in\mathbb{R}^d$. It is again a moment generating function since crossed derivatives (taken a point $\boldsymbol{0}$) are cross-moments. For instance, $\left. \frac{\partial^2 M_X(\boldsymbol t)}{\partial t_i \partial t_j}\right\vert_{\boldsymbol{0}} =\mathbb{E}\left( X_iX_j \right)$Some, moment generating functions are interesting if you want to derive moments of a given distribution. Another interesting feature is that this moment generating function (under certain conditions) fully characterize the distribution of the random variable, in the sense that if for some $h>0$, $M_X(t)=M_Y(t)$ for all $t\in(-h,+h)$, then $X\overset{\mathcal{L}}{=}Y$. • From moment generating functions to characteristic functions The problem with the moment generating function is that the function is defined (only) on some neighborhood of $0$, and we should be careful. The other problem is that it does exist only for distribution in $L_\infty$. Which might be a strong assumption. Thus, an interesting idea is to consider $\mathbb{E}\left( e^{tX} \right)$ not on the real line, but on the imaginary line. Thus, let $\phi_X(t)=\mathbb{E}\left( e^{i tX} \right)$ for some $t\in\mathbb{R}$. Actually, not some, but all $t\in\mathbb{R}$, since $\vert\phi_X(t) \vert=\vert\int e^{i tX} dF_X(x) \vert\leq \int \vert e^{i tX} \vert dF_X(x) =\int dF_X(x)=1$ so the characteristic function always exists. Paul Lévy proved in 1925 that the characteristic function completely characterizes the distribution. Now, if we look at it quickly, it looks like we did not change a lot of things here, and we should be able to write $\phi_X(t)=M_X(i t)$ If we want to do things properly, let us look at Gut (2005) for instance. Assume that $M_X(\cdot)$ is defined on some interval $(-a,+a)$. It is then possible to define a function $\Gamma_X:\mathbb{C}\mapsto \mathbb{C}$ (this time, it is no longer a real-valued function) as $\Gamma_X(z)=\mathbb{E}(e^{z X})$ which is well defined on some strip $\{z\in\mathbb{C}, \vert\text{Re}(z)\vert\< . $\phi_X(\cdot)$and $M_X(\cdot)$ are then restriction of that function respectively on the imaginary line, and the real line. That function $\Gamma_X(\cdot)$ is clearly holomorphic, and thus, the value it takes on such a strip is fully determined by the values it takes on the real interval $(-a,+a)$. Thus, the moment generating function will completely characterize the distribution. But it has to be defined on some neighbourhood of $0$. Which is not trivial actually… I mean, I nonlife insurance, we see a lot a Pareto distributions. • Fast Fourier Transform Recall Euler’s formula, $e^{it}=\cos(t)+i \ \sin(t)$ Thus, we should not be surprised to see Fourier’s transform. From this formula, we can write $\phi_X(t)=\mathbb{E}\left( e^{i tX} \right) =\mathbb{E}\left( \cos[tX] \right)+i\ \mathbb{E}\left( \sin[tX] \right)$ Using some results in Fourier analysis, we can prove that probability function satisfies (if the random variable has a Dirac mass in x) $\mathbb{P}(X=x)=\lim_{T\rightarrow \infty} \frac{1}{2T}\int_{-T}^{+T}e^{-itx}\phi_X(t)dt$ which can also be written, $f_X(x)=\frac{1}{2\pi}\int_{-\pi}^{+\pi}e^{-itx}\phi_X(t)dt$ And a similar relationship can be obtained if the distribution is absolutely continuous at point $x$, $f_X(x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}e^{-itx}\phi_X(t)dt$ Actually, since we work with real-valued random variables, the complex area was just a detour, and we can prove that actually, $f_X(x)=\frac{1}{2\pi}\int_{-\pi}^{+\pi} \text{Re}\left(e^{-itx}\phi_X(t)\right)dt$ It is then possible to get the cumulative distribution function using Gil-Peleaz’s inversion formula, obtained in 1951, $F_X(x)=\frac{1}{2} +\frac{1}{2\pi}\int_{0}^{\infty} \frac{e^{itx}\phi_X(-t)+e^{-itx}\phi_X(t)}{it}dt$ Nice isn’t it. It means, anyone working on financial markets know those formulas, used to price options (see Carr & Madan (1999) for instance). And the good thing is that any mathematical or statistical software can be used to compute those formulas. • Characteristic function and actuarial science Now, what is the interest of all that in actuarial science ? Characteristic functions are interesting when we deal with sums of independent random variables, since the characteristic function of the sum is simple the product of the characteristic functions. They are also interesting when dealing with compound sums1. Consider the problem of computing the 99.5% quantile of the compound sum of Gamma random variable, i.e. $S=\sum_{n=1}^N X_i$ where $X_i\sim\mathcal{G}(\alpha,\beta)$ are i.i.d. and $N\sim\mathcal{P}(\lambda)$. The strategy is to discretize the loss amounts, > n <- 2^20; > p <- diff(pgamma(0:n-.5,alpha,beta)) Then, the code to compute $\tilde f(s)=\mathbb{P}(S\in[s\pm1/2])$, we use > f <- Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n To compute the 99.5% quantile, we just use > sum(cumsum(f)<.995) That’s extremely simple, isn’it. Want me to do it for real ? Consider the following losses amounts > set.seed(1) > X <- rexp(200,rate=1/100) > print(X[1:5]) [1] 75.51818 118.16428 14.57067 13.97953 43.60686 Let us fit a gamma distribution. We can use > fitdistr(X,"gamma") shape rate 1.309020256 0.013090411 (0.117430137) (0.001419982) or > f <- function(x) log(x)-digamma(x)-log(mean(X))+mean(log(X)) > alpha <- uniroot(f,c(1e-8,1e8))$root
> beta <- alpha/mean(X)
> alpha
[1] 1.308995
> beta
[1] 0.01309016

Whatever, we have the parameters of our  Gamma distribution for individual losses. And assume that the mean of the Poisson counting variable is

> lambda <- 100

Again, it is possible to use monte carlo simulations, if we can easily generate a compound sum. We can use the following generic code: first we need functions to generate the two kinds of variables of interest,

> rN.P <- function(n) rpois(n,lambda)
> rX.G <- function(n) rgamma(n,alpha,beta)

then, we can use (see here for a discussion on possible codes)

> rcpd4 <- function(n,rN=rN.P,rX=rX.G){
+ return(sapply(rN(n), function(x) sum(rX(x))))}

If we generate one million variables, we can get an estimator for the quantile,

> set.seed(1)
> quantile(rcpd4(1e6),.995)
99.5%
13651.64

Another idea is to remember a proporty of the Gamma distribution: a sum of independent Gamma distributions is still Gamma (with additional assumptions on the parameters, but here we consider identical Gamma distributions). Thus, it is possible to compute the cumulative distribution function of the compound sum,

> F <- function(x,lambda=100,nmax=1000) {n <- 0:nmax
+ sum(pgamma(x,n*alpha,beta)*dpois(n,lambda))}

(or at least a approximation). If we invert that function, we get our quantile

> uniroot(function(x) F(x)-.995,c(1e-8,1e8))$root [1] 13654.43 Which is consistent with our monte carlo computation. Now, we can also use fast Fourier transform here, > n <- 2^20; lambda <- 100 > p <- diff(pgamma(0:n-.5,alpha,beta)) > f <- Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n > sum(cumsum(f)<.995) [1] 13654 Now, if it is simple, is it efficient ? Let us compare for instance computation time to get those three outputs, > system.time(quantile(rcpd4(1e5),.995)) user system elapsed 2.453 0.106 2.611 > system.time(uniroot(function(x) F(x)-.995,c(1e-8,1e8))$root)
user      system     elapsed
0.041       0.012       0.361
> system.time(sum(cumsum(Re(fft(exp(lambda*(fft(p)-1)),inverse=TRUE))/n)<.995))
user      system     elapsed
0.527       0.020       0.560

Computations here are comparable with the (numerical) inversion of the cumulative distribution function. Except that here, we were lucky: if the distribution is not Gamma but log normal, the second algorithm cannot be used.

1. This numerical example is taken from the first chapter of Computational Actuarial Science with R, to appear in a few months.

# Couples de variables aléatoires

Vendredi, suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur le thème 13 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

# Proba, intra 1

Un rapide billet pour partager le sujet de l’examen intra de la semaine passée avec des éléments de correction (incluant des statistiques de réponse, comme pour la session passée). Toutes les remarques sur mes corrections sont les bienvenues…

# Variables aléatoires continues

Suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur les thèmes 7 et 8 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

Des éléments de correction de l’intra 1 seront bientôt mis en ligne.

# Loi de Poisson

Suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur le thème 6 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

Pour rappels, l’intra du 27 septembre portera sur les thèmes 1-6, c’est à dire sur ce qui a été abordé dans les feuilles mises en ligne.

# Monty Hall (oh no, not again)

Quite frequently, someone on the internet discovers the Monty Hall paradox, and become so enthusiastic that it becomes urgent to publish an article – or a post – about it. The latest example can be http://www.bbc.co.uk/news/magazine-24045598. I won’t blame them, I did the same a few years ago (see http://freakonometrics.hypotheses.org/776, or http://freakonometrics.hypotheses.org/775, in French).

My point today is that the Monty Hall paradox raise an important question, about information. How comes that something to sounds like non-informative can actually be extremely informative. I will not get back on the blue eyes paradox (see http://freakonometrics.hypotheses.org/1963, in French) or the exam paradox (see http://freakonometrics.hypotheses.org/2328, in French one more time), which are related to information, but not with a probabilistic approach. I will stay close to Monty Hall’s paradox today.

This morning, in my probability class, we were looking at a simple exercise (I say simple because it is only the second course of the session). The problem was the following

Consider an urn $A$, with 15 blue balls, and 10 red balls, and an urn $B$, with 10 blue balls, and 15 red balls. We select randomly one urn (with probability 50% for each urn).
We draw a ball, which turns out to be blue, and we put it back in the urn, Now, we draw a (second) ball. What is the probability that this (second) ball is blue?

Ready? Your first thought should be that since we put back the ball, after the first draw, it does not change the probabilities, right? So, why did we say that? It is necessary? (about the last question, yes, when something is mentioned in an exercise, we should use it).

Let’s forget about this second ball story, as an introduction to this problem. What was, actually, the probability for the first ball to be blue? Trivially, it was

$\mathbb{P}(\text{blue}\vert A)\cdot\mathbb{P}( A)+\mathbb{P}(\text{blue}\vert B)\cdot\mathbb{P}( B)$i.e.

$\frac{15}{25}\cdot\frac{1}{2}+\frac{10}{25}\cdot\frac{1}{2}=\frac{1}{2}$

Let us run a code to get that, using simulations:

> n=1000000
> set.seed(1)

First, let us draw the urn, randomly

> urn=sample(1:2,size=n,replace=TRUE)

Then, let us draw the first, and the second ball,

> urns=matrix(c(15,10,10,15),2,2)
> colnames(urns)=c("blue","red")
> sample.urn=(urns[urn,])
> prob.urn=sample.urn/apply(sample.urn,1,sum)
> u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
> u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]

The probability that the first ball was blue is here

> sum(u1=="blue")/n
[1] 0.499953

and for the second one

> sum(u2=="blue")/n
[1] 0.499221

So, indeed, the probability to have a blue ball is 50%. Now, what was the question? Given that the first ball was blue, what it the probability that the second one is blue? Here, on our simulations, it is

> sum(u2[u1=="blue"]=="blue")/sum(u1=="blue")
[1] 0.5194088

Which is close to 52%.And if you run more simulations, you get

> f=function(seed){
+ set.seed(seed)
+ urns=matrix(c(15,10,10,15),2,2)
+ colnames(urns)=c("blue","red")
+ sample.urn=(urns[urn,])
+ prob.urn=sample.urn/apply(sample.urn,1,sum)
+ u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
+ u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]
+ return(sum(u2[u1=="blue"]=="blue")/
+ sum(u1=="blue"))
+ }
> Vectorize(f)(1:20)
[1] 0.5194088 0.5200931 0.5203338 0.5192104 0.5196960 0.5206121 0.5195453
[8] 0.5184580 0.5203755 0.5200154 0.5196557 0.5179276 0.5188652 0.5204724
[15] 0.5197437 0.5209244 0.5205770 0.5208725 0.5206228 0.5190711

The probability is always close to 52%, and is (significantly) different from 50%.

Still not convinced that we have some information here that should be used? Imagine that in the first urn, we add 1 blue ball, and 24 red balls; and the opposite in the second one. In that case, if we say that the first ball was blue, it means that it is very likely that the urn chosen was the second one. Let’s look at by it running some simulations

> set.seed(1)
> urns=matrix(c(1,24,24,1),2,2)
> colnames(urns)=c("blue","red")
> sample.urn=(urns[urn,])
> prob.urn=sample.urn/apply(sample.urn,1,sum)
> u1=c("blue","red")[1+(runif(n)<prob.urn[,1])]
> u2=c("blue","red")[1+(runif(n)<prob.urn[,1])]

As before, the probability that the second ball is blue is 50% (because of the symmetry actually)

> sum(u2=="blue")/n
[1] 0.500362

But if I tell you that the first one was blue, the probability that the second one is blue becomes

> sum(u2[u1=="blue"]=="blue")/sum(u1=="blue")
[1] 0.9236433

So even if – somehow – we do not change much by replacing the ball in its urn, we do have here some information, since it was mentioned that the ball was blue. And we should use it. Again, the important point is that the sentence was not “we draw a ball and we put it back”, but “we draw a blue ball, and we put it back”. Now, it we do the maths, everything become simple, and clear (as usual).

The question is here to compute

$\mathbb{P}(B_2=\text{blue}\vert B_1=\text{blue})$

and according to Bayes formula, it is

$\frac{\mathbb{P}(B_2=\text{blue}\cap B_1=\text{blue})}{\mathbb{P}( B_1=\text{blue})}$

Now, to compute those two probabilities, we have to condition on the urn,

$\mathbb{P}({\color{red}\bullet})=\mathbb{P}({\color{red}\bullet}\vert A)\cdot\mathbb{P}( A)+\mathbb{P}({\color{red}\bullet}\vert B)\cdot\mathbb{P}( B)$

Given the urn, since we replace the ball,

$\mathbb{P}(B_2=\text{blue}\cap B_1=\text{blue}\vert\text{urn})=\mathbb{P}(B_2=\text{blue}\vert\text{urn})\cdot\mathbb{P}(B_1=\text{blue}\vert\text{urn})$

i.e.

$\mathbb{P}(B_2=\text{blue}\cap B_1=\text{blue}\vert\text{urn})=\mathbb{P}(B_1=\text{blue}\vert\text{urn})^2$

So if we substitute numerical probabilities to get a blue ball in the previous formula, we get

$\frac{\displaystyle{\frac{1}{2}\left(\frac{15}{25}\right)^2+\frac{1}{2}\left(\frac{10}{25}\right)^2}}{\displaystyle{\frac{1}{2}\left(\frac{15}{25}\right)+\frac{1}{2}\left(\frac{10}{25}\right)}}$

which not the same as

$\frac{1}{2}\left(\frac{15}{25}\right)+\frac{1}{2}\left(\frac{10}{25}\right)$

Here, we get

> {(15/25)^2+(10/25)^2}/((15/25)+(10/25))
[1] 0.52

which confirms our empirical 52%, and note that in the second case (where there was only 1 blue ball in one urn, and 24 in the second one)

> {(24/25)^2+(1/25)^2}/((24/25)+(1/25))
[1] 0.9232

which again is close to the empirical 92.3% we got.

I strongly believe that the mis-intuition we might have is close to the one we can observe in Monty Hall paradox. And unless you write things properly, it is difficult to conclude anything….

PS [48  hours later] thanks @mikeandallie for the animated version of my post

# Generating a Markov chain vs. computing the transition matrix

A couple of days ago, we had a quick chat on Karl Broman‘s blog, about snakes and ladders (see http://kbroman.wordpress.com/…) with Karl and Corey (see http://bayesianbiologist.com/….), and the use of Markov Chain. I do believe that this application is truly awesome: the example is understandable by anyone, and computations (almost any kind, from what we’ve tried) are easy to perform. At the same time, some French students asked me specific details regarding some old lectures notes on Markov chains, and on some introductory example I used as a possible motivation: the stepping stone algorithm. In the notes, I just mentioned the idea of this popular generic algorithm (introduced in Sawyer (1976)) and I use simulations to show – visually – how it works. Again, it was just to motivate the course which actually did focus on the theory of Markov Chains. But those student wanted more, like how did I get the transition matrix, for instance. And that is actually not a simple question, from a computational perspective. I mean, I can easily generate this Markov Chain, but writing explicitly the transition, that was another story. Which took me a bit longer. In a very specific case…

But let us get back to the roots, and to the stepping stone algorithm. At least, one of them (the one I used in my notes) because it looks like there are several algorithm. We do consider a grid, say $h\times h$, with some colors inside, say $k$ possible colors. Each cell of the grid has a given color. Then, at some stage, we select randomly one cell in the grid, and it will take the color of one of its neighbor (some kind of absorption, or mutation). This is, more or less, what is also detailed in some lecture notes by James Propp (see also e Sato (1983) or Zähle et al. (2005) for more theoretical details about that Markov chain). This is extremely simple to generate (that’s what I did in my notes, with very big grids, and a lot of colors). But what if we want to write the transition matrix ?

First of all, we need to define the state space. Basically, we do have $h^2$ cells, each of them has one color, chosen among $k$. Which gives us $k^{h^2}$ possible states…. And that can be large. I mean, if we consider the smallest possible grid (that might be interesting), say $h=3$, and only $2$ colors, then we talk about $2^9=512$possible states. That is large, not huge. But we should keep in mind that we have to compute a transition matrix, that would be a matrix with $512\times 512=262,144$ elements. More generally, we talk about writing down matrices with $k^{2\cdot h^2}$ elements. If we want black and white $4\times 4$ grids, that would mean a matrix with $2^{2\cdot 4^2}=2^{32}=4,294,967,296$ which mean 4 billion elements ! And if we consider an red-green-blue $3\times 3$ grid, we have to explicit a matrix with $387,420,489$ i.e almost 400 million elements. So, let’s face it: we can only work with $3\times 3$ bi-color grids.

So let’s try… The good thing is that it can be related to work I’ve been doing recently on binomial recombining trees (binomial being related to bi-color). First of all, our grid will be describes as follows

> h=3
> M=matrix(1:(h^2),h,h)
> M
[,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

with two colors

> color=c("red","blue")

Then, we should look for neighbors, or derive an neighborhood matrix,

> d=function(i,j) dist(rbind(c((i-1)%/%h,(i-1)%%h),
+                            c((j-1)%/%h,(j-1)%%h)))
> Neighb=matrix(Vectorize(d)(rep(1:(h^2),each=h^2),
+                            rep(1:(h^2),h^2)),h^2,h^2)
> trunc(Neighb*100)/100
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 1.00 2.00 1.00 1.41 2.23 2.00 2.23 2.82
[2,] 1.00 0.00 1.00 1.41 1.00 1.41 2.23 2.00 2.23
[3,] 2.00 1.00 0.00 2.23 1.41 1.00 2.82 2.23 2.00
[4,] 1.00 1.41 2.23 0.00 1.00 2.00 1.00 1.41 2.23
[5,] 1.41 1.00 1.41 1.00 0.00 1.00 1.41 1.00 1.41
[6,] 2.23 1.41 1.00 2.00 1.00 0.00 2.23 1.41 1.00
[7,] 2.00 2.23 2.82 1.00 1.41 2.23 0.00 1.00 2.00
[8,] 2.23 2.00 2.23 1.41 1.00 1.41 1.00 0.00 1.00
[9,] 2.82 2.23 2.00 2.23 1.41 1.00 2.00 1.00 0.00
> Neighb=(Neighb<2)&(Neighb>0)
> Neighb
[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
[1,] FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE
[2,]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE
[3,] FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE
[4,]  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
[5,]  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
[6,] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE
[7,] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE
[8,] FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
[9,] FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE

Now, let us explicit our 512 possible states.

> n=h^2
> states=function(x){
+   Base.b=rep(0,n)
+   ndigits=(floor(logb(x,base=length(color)))+1)
+   for(i in 1:ndigits){
+     Base.b[n-i+1]=(x%%length(color))
+     x=(x %/% length(color))}
+   return(Base.b)}
> M=Vectorize(states)(1:(length(color)^n-1))
> liststates=data.frame(rbind(rep(0,h^2),t(M)))
X1 X2 X3 X4 X5 X6 X7 X8 X9
1  0  0  0  0  0  0  0  0  0
2  0  0  0  0  0  0  0  0  1
3  0  0  0  0  0  0  0  1  0
4  0  0  0  0  0  0  0  1  1
5  0  0  0  0  0  0  1  0  0
6  0  0  0  0  0  0  1  0  1

(for the first six, with 0/1 digits instead of colors). For instance, if we look at a specific one, it is possible to plot the grid, using

> plotsteps=function(u){
+   plot(0:h,0:h,col="white",xlab="",ylab="",axes=FALSE)
+   for(i in 0:(h^2-1)){
+   x=i%/%h
+   y=i%%h
+   polygon(x+c(1,.1,.1,1),y+c(1,1,.1,.1),
+   col=color[as.numeric(u)[i+1] + 1])
+   text(x+.45,y+.45,i)
+   }}

Here,

> plotsteps(liststates[100,])

Then, given one state, let us see what could happen next,

• let us compute all connected states: all states where we can end up in if we change one cell
• we have to check, for each connect state which cell did change
• we should compute probabilities to reach those 9 states, based on the fact that each of the cell is chosen with the same probability, and the fact that probability to change the color is based on the colors around.
• if some states cannot be reached (if a cell is surrounded by elements of the same color, so it cannot change its color), then, we should remove then from the list of reachable (possible) states.

The code will be something like the following

> listneighbour=function(i){
+   start=liststates[i,]
+   difference2only=function(j) {
+     w=which(liststates[j,]!=liststates[i,])
+     return((length(w)==1))}
+   possible=which( Vectorize(difference2only)(1:nrow(liststates))==TRUE )
+   P=function(j){
+     L=liststates[i,which(Neighb[which(liststates[j,]!=liststates[i,]),]==TRUE)]
+     T=table(as.numeric(L))
+     T=T[as.character(0:(length(color)-1))]
+     T[is.na(T)]=0
+     return(as.numeric(T)/sum(T))
+   }
+   probability=Vectorize(P)(possible)
+   W=NULL
+   for(j in possible) W=c(W,which(liststates[j,]!=liststates[i,]))
+   I=1-liststates[i,W]+1
+   vp=diag(probability[as.numeric(I),])
+   vproba=0*vp
+   if(sum(vp)!=0) vproba=vp/sum(vp)
+   return(list(
+     color=liststates[i,W],
+     absorb=W,
+     possible=possible,
+     probability=probability,
+     prob=vproba))
+ }

For instance, if we start from state 100 (here, on the right)

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

$possible [1] 36 68 98 99 104 108 116 228 356$probability
[,1] [,2] [,3]   [,4]   [,5] [,6] [,7] [,8]   [,9]
[1,]    1  0.8  0.6 0.6667 0.3333  0.4  0.5  0.6 0.6667
[2,]    0  0.2  0.4 0.3333 0.6667  0.6  0.5  0.4 0.3333

$prob [1] 0.17964072 0.14371257 0.10778443 0.11976048 0.11976048 [6] 0.10778443 0.08982036 0.07185629 0.05988024 Let us look more specificaly at the 99th state (which appears above as a state that could be reached from the 100th), > liststates[99,] X1 X2 X3 X4 X5 X6 X7 X8 X9 99 0 0 1 1 0 0 0 1 0 If we plot it (here on the right, again), we get > plotsteps(liststates[99,]) Clearly, here, the cell in the upper corner (number 9) changed from blue to red. Now, about the probability… The probability to select cell 9 is 1/9, and given that cell 9 is chosen, the probability to go from blue to red is 2/3 (the cell is surrounded by 2 red cells, and 1 blue cell). The probability to remain blue is then 1/3. Those are the probabilities computed by our function (the table with two rows, one per color). In order to get a better understanding on the meaning of the last line, with some sort of probabilities), let us look at the following (simpler) example. > liststates[2,] X1 X2 X3 X4 X5 X6 X7 X8 X9 2 0 0 0 0 0 0 0 0 1 that can be visualized on the right (on the right). Here, > listneighbour(2)$color
X9 X8 X7 X6 X5 X4 X3 X2 X1
2  1  0  0  0  0  0  0  0  0

$absorb [1] 9 8 7 6 5 4 3 2 1$possible
[1]   1   4   6  10  18  34  66 130 258

$probability [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1 0.8 1 0.8 0.875 1 1 1 1 [2,] 0 0.2 0 0.2 0.125 0 0 0 0$prob
[1] 0.65573770 0.13114754 0.00000000 0.13114754 0.08196721
[6] 0.00000000 0.00000000 0.00000000 0.00000000

Things are pretty simple here

• if we chose cells $\{1,2,3,4,7\}$, then nothing change, since all the neighbors have the same color. So if we want to focus on changes (or say run the algorithm until the first color change, then choosing those cells is a waste of time)
• if we chose cells $\{5,6,8\}$, then it could be possible to change the color. And actually, $\{5\}$ is different from $\{6,8\}$ (since it does have much more neighbors)
• if we chose cell $\{9\}$, then definitively, the color will change, since all neighbors have the other color here,

Now, the probability to select cell $k$ given that there was a color change would be, if $k$ is in $\{9\}$

$\mathbb{P}(k)\propto \frac{3}{3}=1$

while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so

$\mathbb{P}(k)\propto \frac{1}{5}$and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so

$\mathbb{P}(k)\propto \frac{1}{8}$

And for the other, $\mathbb{P}(k)\propto 0$. So, it comes – since we assume that cells are drawn independently, and with the same probability, if $k$ is in $\{9\}$

$\mathbb{P}(k)= \frac{1 \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{40}{61}$

while if $k$ is in $\{6,8\}$, then there are 4 out 5 neighbors that are red, so

$\mathbb{P}(k)= \frac{\frac{1}{5} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{8}{61}$

and if $k$ is $\{5\}$, then, only one neighbor has a different color, out of 8, so

$\mathbb{P}(k)= \frac{\frac{1}{8} \cdot \frac{1}{9}}{\left(1+2\times \frac{1}{5}+ \frac{1}{8}+5\times 0\right)\cdot \frac{1}{9}}=\frac{5}{61}$

Which are exactly the probability computed above. The point is that we compute probabilities given that a color change will actually occur. The good point is that it should faster convergence to some limiting distribution. If any.

What about our transition matrix ? Well, using a simply loop, we should get it easily

> M=matrix(0,nrow(liststates),nrow(liststates))
+ for(i in 1:nrow(liststates)){
+ L=listneighbour(i)
+ if(sum(L$prob)!=0){ + j=L$possible
+ M[i,j]=L$prob + } + if(sum(L$prob)==0){
+ j=i
+ M[i,j]=1
+ }
+ }

One can check that this matrix satisfies some properties of transition matrices. For instance, the sum per row is one,

> sum(apply(M,1,sum)!=1)
[1]  0

Remember that this matrix is big, so I will not print if here. But trust me, it works (it might take a while on an old laptop, but anyone can do it). Now, if we want to visualize some paths of that chain, we can use the following algorithm. First, we need a starting point, that can be chosen randomly,

> j=sample(1:nrow(liststates),size=1)

or using a given colored grid, say

> j=100

Then we plot it,

> plotsteps(liststates[j,])

Now, the code within the loop is here

> d=rep(0,nrow(liststates))
> d[j]=1
> d=d%*%M
> j=sample(1:nrow(M),size=1,prob=d)
> plotsteps(liststates[j,])

Here are some examples. And indeed, we end up either with all cells in blue, or all cells in red.

Now, do we have to compute that transition matrix to produce those graph (and to generate that Markov chain) ? No. Of course not… At each step, I use a Dirac measure, and use the transition matrix just to get the probability to generate then the next state. Actually, one can write a faster and more intuitive code to generate the same chain… But I should probably keep that for another post…

# From Simpson’s paradox to pies

Today, I wanted to publish a post on economics, and decision theory. And probability too… Those who do follow my blog should know that I am a big fan of Simpson’s paradox. I also love to mention it in my
econometric classes. It does raise important questions, that I do relate to multicolinearity, and interepretations of regression models, with multiple (negatively correlated) explanatory variables. This paradox has amazing pedogological virtues. I did mention it several times on this blog (I should probably mention that I discovered this paradox via Marco Scarsini, who did learn me a lot of things, in decision theory and in probability). For those who do not know this paradox, here is an example that Marco gave in one of his talk, a few years ago. Consider the following statistics, when healthy people entered in some hospital

 hospital total survivors deaths survival rate hospital A 600 590 10 98% hospital B 900 870 30 97%

while, when sick people entered in the same hospitals

 hospital total survivors deaths survival rate hospital A 400 210 190 53% hospital B 100 30 70 30%

Somehow, whatever your health situation, you should choose hospital A. Now, if we agregate

 hospital total survivors deaths survival rate hospital A 1000 800 200 80% hospital B 1000 900 100 90%

i.e. without any doubts, people should choose hospital B.

and

while

With symbolic notations, one can have at the same time

and

with also

as shown on the graph below

There should be connection between Simpson’s paradox and the ecological fallacy (which is an issue I recently discovered and that I found extremely interesting, related again to difficulties of interpreting
regressions). But that’s another story. My point today is that Colin Blyth did mention another nice paradox, that is related, this time, to stochastic orderings. The idea is the following. Consider the three spinners drawn below (imagine some arrows in those circles)

• spinner A: no matter where the arrow stops, the gain is 3,
• spinner B: 56% chances to gain 2, 22% chances to gain 4, and 22% chances to gain 6,
• spinner C: 51% chances to gain 1, 49% chances to gain 5.

Instead of spinners, it is also possible to consider three different lotteries,

You play against a friend, you pick a spinner, while the friend picks another. Everyone flick his arrow, the highest number wins (no matter the difference). Let us compute the odds. First case, A against B, from
A’s perspective

 B-2 B-4 B-6 A-3 56% +1 win 22% -1 lose 22% -3 lose

In that case, A has 56% chance of beating B. Second case, A against C, from A’s perspective,

 C-1 C-5 A-3 51% +1 win 49% -2 lose
In that case, A has 51% chance of beating C. Third (an final) case, B against C, from B’s perspective. Assuming independence between the spinners, joint probabilities can easily be computed,
 C-1 C-5 B-2 28.56% +1 win 27.44% -3 lose B-4 11.22% +3 win 10.78% -1 lose B-6 11.22% +5 win 10.78% +1 win
In that case, B has 61.78% chance of beating C. So, if we try to summarize,
• A is the best choice, since it beats both with – always – more than 50% chance,
• C is the worst choice, since it is beaten by both with – always – more than 50% chance,
Now, assume that you play not against one friend, but two friends. An everyone picks a different spinner. Let
us compute the odds, one more time. First case, A against B and C, from A’s perspective
 B-2 C-1 B-2 C-5 B-4 C-1 B-4 C-5 B-6 C-1 B-6 C-5 A-3 28.56% +1 win 27.44% -2 lose 11.22% -1 lose 10.78% -1 lose 11.22% -3 lose 10.78% -3 lose
In that case, A has 28.56% chance of beating B and C. Second case, B against A and C, from B’s perspective,
 A-3 C-1 A-3 C-5 B-2 28.56% -1 lose 27.44% -2 lose B-4 11.22% +1 win 10.78% -1 lose B-6 11.22% +3 win 10.78% +1 win
In that case, B has 33.22% chance of beating A and B.Third (an final) case, C against A, from C’s perspective,
 A-3 B-2 A-3 B-4 A-3 B-6 C-1 28.56% -2 lose 11.22% -3 lose 11.22% -5 lose C-5 27.44% +2 win 10.78% +1 win 10.78% -1 lose

In that case, C has 38.22% chance of beating A and B. So, if we try to summarize, this time

• C is the best choice, since has (strictly) more than 1/3 chances to win, which the highest probability
• A is the worst choice, since has (strictly) less than 1/3 chances to win, which the lowest probability

Odd isn’t it ? Now, is there an interpretation of that paradox ? Yes, Martin Gardner, in his paper on induction and probability, mentioned the case of drug testing. The value we had with the spinner is the health level, rated from 1 to 6. Thus, taking drug A, you always get an average health level of 3. With drug C, on the other hand, you get either very sick (level 1) or very well (level 5). Consider now a doctor who wants to maximize the patient’s chance of being well. If only pills A and C are available, then the doctor should choose A. This is what we’ve seen in the first part. Assume that now a company delivers a third pill, called drug B. Then the doctor should find C more interesting…. Odd, isn’t it ?

Colin Blyth gave a more amusing application. Assume that you like to go to the restaurant, and you like get a dessert there. Dessert A – the apple pie – is the average one, with a standard level, that you rank 3 (on a scale from 1 to 6). Dessert C – the cheese cake – can either be awfull (ranked 1) or delicious (ranked 5). You’d better go for the apple pie if you want to maximize the probability of not being disappointed (i.e. maximizing your “best chance” according to Colin Blyth, but I guess it can be interpreted as regret minimization too). Now assume that dessert B – the blueberry pie – is available (with ranks given by the spinner). Then you should go for the cheese cake. I let you imagine the discussion that you can have, then, with your favorite waitress

– Hi Mr Freakonometrics, do you want a piece of apple pie ? (yes, actually she also comes frequently on my blog, and knows me from my pseudo…)

– Probably. But actually, I was wondering if you did have your blueberry pie today ?

– Yes, in fact we do….

– Great, in that case, I’ll go for the cheese cake.

She’ll probably think that I am freak… so I hope she’ll come and read my post, to understand that, actually, it does make a lot of sense to go for what was supposed to be my worst case.

# Pills, half pills and probabilities

Yesterday, I was uploading some old posts to complete the migration (I get back to my old posts, one by one, to check links of pictures, reformating R codes, etc). And I re-discovered a post published amost 2 years ago, on nuns and Hell’s Angels in an airplaine.

It reminded me an old probability problem (that might be known as one on Feymann’s problems): suppose that you have a prescription to take half pills for 6 days. Unfortunately the pharmacist was a bit lazy (or just wanted to help me to write a mathematical problem), and he gives 3 (full) pills in a small box. Day 1, you take a pill, break it in two parts, eat one, and return the other half in the box. Day 2, you draw randomly ‘something’ from the box, i.e. either half a pill, or a pill. If it’s a half one, then you eat it. If it is a fill one, you break it in two, eat one half, and return the other half in the box. Etc.On Day 6, if my story was well explained, you should know that there can only be one half pill. So far, so good. But what about Day 5 ? There were either two half pills, or one full pill. But what was the probability that there was a fill pill in the box on Day 5 ?

Nice problem, isn’t it ?

The good thing is that it can be modeled as a Markovian model. Assume that we do have $n$ pills. After $2n$ days, the box will be empty. Consider the pair $(h,c)$ denoting the number of half pills, and complete pills. $c$ can take all values, from 0 to $n$, and $h$ will be positive, with $c+h\leq n$. Thus, the number of states – possible pairs from Day 1 till Day $2n$ – will be $1+2+\cdots+(n+1)$, i.e. $(n+1)(n+2)/2$. More precisely, define those states in a dataframe,

> n=3
> COMPLETE=HALF=NULL
> for(i in n:0){
+ HALF=c(0:(n-i),HALF)
+ COMPLETE=c(rep(i,length(0:(n-i))),COMPLETE)
+ }
> k=length(COMPLETE)
> state=data.frame(s=1:k,nc=rev(COMPLETE),nh=rev(HALF))
> state
s nc nh
1   1  3  0
2   2  2  1
3   3  2  0
4   4  1  2
5   5  1  1
6   6  1  0
7   7  0  3
8   8  0  2
9   9  0  1
10 10  0  0

Now, we can play to derive the transition matrix of the Markov chain.

> attach(state)
> P=matrix(0,k,k)
> for(i in 1:k){
+ C=state$nc[i] + H=state$nh[i]
+ if((C>0)&(H>0)){
+ P[i,state[(nc==C-1)&(nh==H+1),"s"]]= C/(C+H)
+ P[i,state[(nc==C)&(nh==H-1),"s"]]= H/(C+H)}
+ if((C>0)&(H==0)){
+ P[i,state[(nc==C-1)&(nh==H+1),"s"]]=1}
+ if((C==0)&(H>0)){
+ P[i,state[(nc==C)&(nh==H-1),"s"]]=1}
+ if((C==0)&(H==0)){
+ P[i,state[(nc==C)&(nh==H),"s"]]=1}
+ }

We do have a transition matrix (or a probability matrix) since all elements are positive, and the sum per line is 1,

> apply(P,1,sum)
[1] 1 1 1 1 1 1 1 1 1 1

Here, the transition matrix is the following

> P
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    1 0.00 0.00 0.00  0.0 0.00  0.0    0     0
[2,]    0    0 0.33 0.66 0.00  0.0 0.00  0.0    0     0
[3,]    0    0 0.00 0.00 1.00  0.0 0.00  0.0    0     0
[4,]    0    0 0.00 0.00 0.66  0.0 0.33  0.0    0     0
[5,]    0    0 0.00 0.00 0.00  0.5 0.00  0.5    0     0
[6,]    0    0 0.00 0.00 0.00  0.0 0.00  0.0    1     0
[7,]    0    0 0.00 0.00 0.00  0.0 0.00  1.0    0     0
[8,]    0    0 0.00 0.00 0.00  0.0 0.00  0.0    1     0
[9,]    0    0 0.00 0.00 0.00  0.0 0.00  0.0    0     1
[10,]   0    0 0.00 0.00 0.00  0.0 0.00  0.0    0     1

In order to get our probability, let us start from state 1 – or $(c,h)=(n,0)$ – with probability 1, and let us look at the distribution at different periods,

> dist=c(1,rep(0,k-1))
> MatDist=matrix(NA,2*n+1,k)
> MatDist[1,]=dist
> for(i in 1:(2*n)){dist=as.vector(t(dist)%*%P)
+ MatDist[i+1,]=dist
+ }

(one can check that after $2n$ days, the box is empty). The probability is given in row $2n-1$, and we just have to check which column corresponds to the pair $(c,h)=(1,0)$,

> vs=state[which(MatDist[2*n-1,]>0),]
> proba=MatDist[2*n-1,vs[vs$nc==1,"s"]] > proba [1] 0.3888889 Here the probability of having a full pair on Day 5 is 38.89%. Actually, it is possible to study the evolution of this probability as a function of $n$, > computeproba=function(n=3){ + COMPLETE=HALF=NULL + for(i in n:0){ + HALF=c(0:(n-i),HALF) + COMPLETE=c(rep(i,length(0:(n-i))),COMPLETE) + } + k=length(COMPLETE) + state=data.frame(s=1:k,nc=rev(COMPLETE),nh=rev(HALF)) + P=matrix(0,k,k) + for(i in 1:k){ + C=state$nc[i]
+ H=state$nh[i] + if((C>0)&(H>0)){ + P[i,state[(state$nc==C-1)&(state$nh==H+1),"s"]]= C/(C+H) + P[i,state[(state$nc==C)&(state$nh==H-1),"s"]]= H/(C+H)} + if((C>0)&(H==0)){ + P[i,state[(state$nc==C-1)&(state$nh==H+1),"s"]]=1} + if((C==0)&(H>0)){ + P[i,state[(state$nc==C)&(state$nh==H-1),"s"]]=1} + if((C==0)&(H==0)){ + P[i,state[(state$nc==C)&(state$nh==H),"s"]]=1} + } + dist=c(1,rep(0,k-1)) + MatDist=matrix(NA,2*n+1,k) + MatDist[1,]=dist + for(i in 1:(2*n)){dist=as.vector(t(dist)%*%P) + MatDist[i+1,]=dist + } + vs=state[which(MatDist[2*n-1,]>0),] + proba=MatDist[2*n-1,vs[vs$nc==1,"s"]]
+ return(proba)
+ }

If we plot the probability as a function of $n$, we get

> P=Vectorize(computeproba)(2:40)
> plot(2:40,P,ylim=c(0,.5))

One can observe that the probability is decreasing. But slowly, extremely slowly. With a log scale on the y-axis, we have

> plot(2:40,P,ylim=c(0,.5),log="y")

If we look for ‘high’ values, we can get

> computeproba(100)
[1] 0.14218

I do not know if this limit goes to 0 as $n$ goes to infinity. Actually, since we do have to compute a matrix with $[(n+1)(n+2)/2]^2$  entries i.e. roughly $n^4/4$, $n$ cannot be that large… Too bad. If anyone knows how this probability behaves as a function of $n$, when $n$ is large, I’d be glad to know…

# UEFA, is that it ?

Following my previous post, a few more things. As mentioned by Frédéric, it is – indeed – possible to compute the probability of all pairs. More precisely, all pairs are not as likely to occur: some teams can play against (almost) eveyone, while others cannot. From the previous table, it is possible to compute probability that the last team plays against team 1. Or team 2 (numbers are from the  xls file mentioned previously). To make it simple

> table(M[,2*n])/length(M[,2*n])*100

1        2        3        5        7       10       11
11.82500 12.61212 12.61212 13.25279 19.31173 18.70767 11.67856

Here, the last team (as I did rank them) has 11.8% chances to play against team 1, and 19.3% to play against team 7. If we compute all the probabilities, we obtain

> S
1     2     3     5     7    10    11    13
4   0.00 14.16 14.16  0.00 22.22 21.25 13.05 15.13
6  12.52 13.19 13.19 14.11 20.13  0.00 12.35 14.47
8  18.78  0.00 19.54 21.50  0.00  0.00 18.39 21.76
9  18.78 19.54  0.00 21.50  0.00  0.00 18.39 21.76
12 14.68 15.54 15.54 16.56  0.00 23.19 14.47  0.00
14 11.64 12.37 12.37 13.05 18.96 18.25  0.00 13.34
15 11.77 12.55 12.55  0.00 19.36 18.59 11.64 13.50
16 11.82 12.61 12.61 13.25 19.31 18.70 11.67  0.00

that can be visualized below

White areas cannot be reached, while red ones are more likely. Here, we compute probability that home team (given on the x-axis) plays against some visitor team (on the y-axis). The fact that those probabilities are not uniform seems odd. But I guess it comes from those constraints…

Another weird point: it is possible to reach a deadlock. At least with the technique I have been using. So far, I did not count them. But we can, simply the following code

> U=c(4,6,8,9,12,14,15,16)
> a1=U[1]
> b1=U[2]
> c1=U[3]
> d1=U[4]
> e1=U[5]
> f1=U[6]
> g1=U[7]
> h1=U[8]
> a2=b2=c2=d2=e2=f2=g2=h2=NA
> posa2=(1:n)%notin%c(LISTEIMPOSSIBLE[,a1])
> if(length(posa2)==0){na=na+1}
> for(a2 in posa2){
+ posb2=(1:n)%notin%c(LISTEIMPOSSIBLE[,b1],a2)
+ if(length(posb2)==0){na=na+1}
+ for(b2 in posb2){
+ posc2=(1:n)%notin%c(LISTEIMPOSSIBLE[,c1],a2,b2)
+ if(length(posc2)==0){na=na+1}
+ for(c2 in posc2){
+ posd2=(1:n)%notin%c(LISTEIMPOSSIBLE[,d1],
+ a2,b2,c2)
+ if(length(posd2)==0){na=na+1}
+ for(d2 in posd2){
+ pose2=(1:n)%notin%c(LISTEIMPOSSIBLE[,e1],
+ a2,b2,c2,d2)
+ if(length(pose2)==0){na=na+1}
+ for(e2 in pose2){
+ posf2=(1:n)%notin%c(LISTEIMPOSSIBLE[,f1],
+ a2,b2,c2,d2,e2)
+ if(length(posf2)==0){na=na+1}
+ for(f2 in posf2){
+ posg2=(1:n)%notin%c(LISTEIMPOSSIBLE[,g1],
+ a2,b2,c2,d2,e2,f2)
+ if(length(posg2)==0){na=na+1}
+ for(g2 in posg2){
+ posh2=(1:n)%notin%c(LISTEIMPOSSIBLE[,h1],
+ a2,b2,c2,d2,e2,f2,g2)
+ if(length(posh2)==0){na=na+1}
+ for(h2 in posh2){
+ s=s+1
+ V=c(a1,a2,b1,b2,c1,c2,d1,d2,e1,e2,f1,f2,g1,g2,h1,h2)
+ }}}}}}}}

On the initial ordering of home team, the number of deadlocks was

> na
[1] 657

The probability of obtaining a deadlock is then

> 657/(657+5463)
[1] 0.1073529

(657 scenarios ended in a dead end, while 5463 ended well). The worst case was obtained when we considered

 [1]    6    4   16   14   12   15    8    9

In that case, the probability of obtaining a deadlock was

> 4047/(4047+5463)
[1] 0.4255521

Here, it clearly depends on the ordering. So if we draw – randomly – the order of the home teams, i.e.

> Urandom=sample(U,size=8)

the distribution of the probablity of having a deadlock is

All those computations were based on my understanding of the drawings. But Kristof (aka @ciebiera), on his blog krzysztofciebiera.blogspot.ca/… obtained different results. For instance, based on my previous computations, the probability to obtain identical pairs was 0.018349% (1 chance out of 5463), but Kristof obtained – based on the UEFA procedure (as he called it) – a probability of 0.0181337%. Which is not _ strictly – the same, but both computations yield relatively close results…

# UEFA, what were the odds ?

Ok, I was supposed to take a break, but Frédéric, professor in Tours, came back to me this morning with a tickling question. He asked me what were the odds that the Champions League draw produces exactly the same pairings from the practice draw, and the official one (see e.g. dailymail.co.uk/…).

To be honest, I don’t know much about soccer, so here is what happened, with the practice draw (on the left, on December 19th) and the official one (on the right, on December 20th),

Clearly, the pairs are identical, but not the order. Actually, at first, I was suprised that even which team plays at home first, was iddentical. But (it seams that) teams that play at home first are the ones that ended second after the previous stage of the competition.

And to be more specific about those draws, those pairs were obtained using real urns, real balls, so it is pure randomness (again, as far as I understood). But with very specific rules. For instance, two teams from the same country cannot play together (or one against the other) at this stage. Or teams that ended first after the previous turn can only play with (or against) teams that ended second. Actually, Frederic sent me an xls file, with a possibility matrix.

Let us find all possible pairs, regardless which team plays at home first (again, we do not care here since the order is defined by the rule mentioned above). Doing the maths might have been a bit complicated, with all those contraints. With a small code, it is possible to list all possible pairs, for those eight games. Let us import our possibility matrix,

 > n=16
+ "http://freakonometrics.blog.free.fr/public/data/uefa.csv",
> LISTEIMPOSSIBLE=matrix(
+ (rep(1:n,n))*(uefa[1:n,2:(n+1)]=="NON"),n,n)

I can fix the first team (in my list, the fourth one is the first team that ended second). Then, I look at all possible second one (that will play with the first one),

 > a1=1
> "%notin%" <- function(x, table){x[match(x, table, nomatch = 0) == 0]}
> posa2=((a1+1):n)%notin%LISTEIMPOSSIBLE[,a1]

Then, consider the second team that ended second (the sixth one in my list). And look at all possible fourth team (that will play this second game), i.e exluding the one that were already drawn, and those that are not possible,

 > b1=6
> posb2=(1:n)%notin%c(LISTEIMPOSSIBLE[,b1],a2)

Etc. So, given the list of home teams,

 > a1=4
> b1=6
> c1=8
> d1=9
> e1=12
> f1=14
> g1=15
> h1=16

consider the following loops,

 > posa2=(1:n)%notin%c(LISTEIMPOSSIBLE[,a1])
> for(a2 in posa2){
+ posb2=(1:n)%notin%c(LISTEIMPOSSIBLE[,b1],a2)
+ for(b2 in posb2){
+ posc2=(1:n)%notin%c(LISTEIMPOSSIBLE[,c1],a2,b2)
+ for(c2 in posc2){
+ posd2=(1:n)%notin%c(LISTEIMPOSSIBLE[,d1],a2,b2,c2)
+ for(d2 in posd2){
+ pose2=(1:n)%notin%c(LISTEIMPOSSIBLE[,e1],a2,b2,c2,d2)
+ for(e2 in pose2){
+ posf2=(1:n)%notin%c(LISTEIMPOSSIBLE[,f1],a2,b2,c2,d2,e2)
+ for(f2 in posf2){
+ posg2=(1:n)%notin%c(LISTEIMPOSSIBLE[,g1],a2,b2,c2,d2,e2,f2)
+ for(g2 in posg2){
+ posh2=(1:n)%notin%c(LISTEIMPOSSIBLE[,h1],a2,b2,c2,d2,e2,f2,g2)
+ for(h2 in posh2){
+ s=s+1
+ V=c(a1,a2,b1,b2,c1,c2,d1,d2,e1,e2,f1,f2,g1,g2,h1,h2)
+ cat(s,V,"\n")
+ M=rbind(M,V)
+ }}}}}}}}

With the print option, we end up with

5461 4 13 6 11 8 5 9 2 12 10 14 3 15 7 16 1
5462 4 13 6 11 8 5 9 2 12 10 14 7 15 1 16 3
5463 4 13 6 11 8 5 9 2 12 10 14 7 15 3 16 1

i.e.

> nrow(M)
[1] 5463

possible pairs (the list can be found here, where numbers are the same as the one in the csv file). Which was the probability mentioned in acomment in the article mentioned previously dailymail.co.uk/…. So the probability to have exactly the same output after the practise and the official draws was (in %)

> 100/nrow(M)
[1] 0.01830496

Which is not that small when we think about it….

And if someone has a mathematical expression for this probability, I am interested. The only reliable method I found was to list all possible pairs (the csv file is available if someone wants to check). But I am not satisfied….

# Actuariat 1, ACT2121, huitième cours

Pour le huitième cours d’actuariat 1 (ACT2121, préparation à l’examen P de la SOA), on continuera les exercices commencés la semaine passée. Je mets toutefois en ligne quelques exercices supplémentaires, pour ceux qui souhaitent s’entraîner davantage (le fichier est en ligne ici). Pour rappel (?) l’examen final aura lieu dans 2 semaines la semaine prochaine, et portera sur l’ensemble de la matière. Comme toujours, 30 questions, 3 heures, et on commence à 13 heures (dois-je le préciser ?). Cette fois, je fournis la table “officielle” de la SOA.