# Monty Hall problem, with Thompson sampling

We all know the Monty Hall problem. Recently, Jason Rosenhouse published a book on that topic (entitled The Monty Hall Problem, The Remarkable Story of Math’s Most Contentious Brain Teaser). The game is more or less described by the following question

Suppose you’re on a game show, and you’re given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what’s behind the doors, opens another door, say No. 3, which has a goat. He then says to you, “Do you want to pick door No. 2?” Is it to your advantage to switch your choice?

While I was preparing some slides for a lecture on Bayesian modeling and thinking, I wanted to find an illustration of what is sometimes called the Bayesian brain, that can be related to updates of beliefs, when we experience. And I was looking for examples of Thompson sampling. And actually, it is possible to learn that switching is the optimal strategy, in the Monty Hall problem, just by playing sequentially the game, and learning from previous strategies. The following code is used, to choose the door with the price (the car), and the one we first select

set.seed(1) n = 5000 listdoor = matrix(1:3,3,n) door = listdoor win = sample(1:3,size=n,replace=TRUE) pick1 = sample(1:3,size=n,replace=TRUE)

Then, the presenter picks one, that is neither the car, nor the one we chose initially. The following trick can be used, to get the list of available choices

door[win+(0:(n-1))*3] = NA door[,1:10] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] NA NA NA 1 NA NA 1 NA 1 NA [2,] 2 2 NA NA 2 2 2 NA NA 2 [3,] 3 NA 3 3 NA NA NA 3 NA NA door[pick1+(0:(n-1))*3] = NA door[,1:10] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] NA NA NA 1 NA NA 1 NA 1 NA [2,] 2 2 NA NA 2 2 2 NA NA 2 [3,] 3 NA 3 3 NA NA NA 3 NA NA

Then, the presenter picks one

presenter = apply(door,2, function(x) sample(x[!is.na(x)],size=1)) &gt; presenter[win != pick1] = apply(door,2,function(x) x[!is.na(x)])[win != pick1] presenter = unlist(presenter) presenter[1:10] [1] 3 2 3 1 2 2 2 3 1 2

Now, let us consider the  Monty Hall problem. We have two possible strategies. The first one is to keep the door we chose, initially

pick2a = pick1 gaina = (pick2a==win) mean(gaina) [1] 0.3392

As expected, on average, we win with (about) 1 chance out of 3. The second one is to (always) pick the other door (the one left). The code is close to the one we used before

door = listdoor door[pick1+(0:(n-1))*3] = NA door[presenter+(0:(n-1))*3] = NA pick2b = apply(door,2,function(x) x[!is.na(x)]) gainb = (pick2b==win) mean(gainb) [1] 0.6608

If you know Monty Hall problem the probability to win is now 2 chance out of 3 (which is what the maths tells us). That is what we have with simulations.

Now, what if we don’t know how to do the maths, and we don’t want to compute it? We can use Thompson sampling to explore, and exploit. In a general context, we have to choose among On a le choix entre $K$ alternatives (here $K=2$, since we can either keep our initial choice, or pick the other one), and the output is $\boldsymbol{X}=(X_1,\cdots, X_K)$, where $X_k\sim\mathcal{B}(\theta_k)$, but $\theta_k$ is unknow, and we will play the game, and learn. From previous computations, we know that $\theta_1=1/3$ while $\theta_2=2/3$.

We use some prior distribution, $\theta_k\sim\mathcal{B}eta(\alpha_k,\beta_k)$, since the Beta distribution is the conjugate of the Bernoulli. At time $t$, we draw $K$ (independent) Beta variables $B_k\sim\mathcal{B}eta(\alpha_k,\beta_k)$, and pick $k^\star = \displaystyle{\underset{k=1,\cdots,K}{\text{argmax}}\{B_k\}}$.  Here the code will be

set.seed(2) X = cbind(pick2a == win,pick2b == win)*1 AB1 = AB2 = tirage = matrix(NA,n,2) choix = rep(NA,n) k=1 AB1[k,] = AB2[k,] = c(1,1) for(k in 1:(n-1)){ tirage[k,] = c(rbeta(1,AB1[k,1],AB1[k,2]), rbeta(1,AB2[k,1],AB2[k,2])) choix[k] = which.max(tirage[k,]) if(choix[k] == 1){ AB1[k+1,] = AB1[k,] + c(X[k,1],1-X[k,1]) AB2[k+1,] = AB2[k,] } if(choix[k] == 2){ AB1[k+1,] = AB1[k,] AB2[k+1,] = AB2[k,] + c(X[k,2],1-X[k,2]) }}

Before showing some graphs, let us check that indeed, we select more the second strategy (which is here to select the other door)

AB1[n,] [1] 5 13 AB2[n,] [1] 3292 1693

Indeed, since the average of a Beta distribution, $\mathcal{B}eta(\alpha,\beta)$ is $\alpha/(\alpha+\beta)$

AB2[n,1]/(sum(AB2[n,])) [1] 0.6603811

i.e. the probability to win, with this second strategy is about $2/3$ (as obtained previously). We can visualize this on the animation below, with, in red the first strategy (keep your initial choice), in green the second one (select the other door), 0 and 1 respectively if we win, or not. Then we can visualize the evolution of $\alpha_2$ and $\beta_2$ on topc, and $\alpha_1$ and $\beta_1$ below (the index is time $t$). Finallly, we have the two variables $B_1$ and $B_2$ drawn,

Of course, another simulation would have given different $B_1$‘s and $B_2$‘s, but finally, we learn that the second strategy is better, and we learn it quite fast…

Here is another one (just to confirm)

So clearly, even if we don’t know which is the optimal strategy (keep our initial choice, or switch), a player who played that game about 30 times should be able to understand that switching should be a better strategy.

# Tests, Power and Significance

In the mathematical statistics course today, we started talking about tests, and decision rules. To illustrate all the concepts introduced today, we considered the case where we have a sample $\boldsymbol{x}=\{x_1,\cdots,x_n\}$ with $U_i\sim\mathcal{U}([0,\theta])$. And we want to test

$H_0:\theta\leq \theta_0$  against $H_1:\theta> \theta_0$

In the course, we’ve seen that we could use a test based on the order statistics $x_{n:n}=\max\{x_i\}$.  The test would be

$\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(x_{n:n})$

i.e. if $\psi(\boldsymbol{x})=1$ we choose $H_1$, and if $\psi(\boldsymbol{x})=0$, we choose $H_0$.

From the definition of the first order risk,

$\alpha=\sup_{\theta\in\Theta_0}\left\lbrace{\mathbb E}_{\theta}\lbrack \psi(\boldsymbol{X})\rbrack\right\rbrace={\mathbb E}_{\theta_0}\lbrack \psi(\boldsymbol{X})\rbrack$

we can easily get that

$c=\theta_0\cdot(1-\alpha)^{\frac{1}{n}}$

Thus, the power is then

$p(\theta)=\left(1-\left(\frac{\theta_0}{\theta}\right)^n(1-\alpha)\right)\boldsymbol{1}_{(c\ ;\ +\infty)}(\theta)$

To visualize it, use the following parameters

n=5
alpha=.1
theta0=1

Then

C1=theta0*(1-alpha)^(1/n)
theta=seq(0,2,by=.01)
P1=(1-(theta0/theta)^n*(1-alpha))*(theta>C1)
plot(theta,P1,type="l",lwd=2,col="blue",xlab="",ylab="Power")

Note that, so far, we did never consider the maximum of our sample. Assume that the maximum is $x_{n:n}$, then we can compute the $p$-value,

$p=\mathbb{P}(X_{n:n}>x_{n:n})=1-\left(\frac{x_{n:n}}{\theta_0}\right)^n$

Here it is

PV=(1-theta^n)*(theta<=1)
plot(theta,PV,type="l",lwd=2,col="blue",xlab="",ylab="p-value")

Now, why not consider another test, based on the minimum (since we have the distribution of the minimum of a sample from a uniform distribution). The test is the same as before

$\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(x_{1:n})$

but here, the threshold is

$c=\theta_0\cdot (1-\alpha^{\frac{1}{n}})$

The power of the test is here

$\quad p(\theta)=\left(1-\frac{\theta_0}{\theta}(1-\alpha^{\frac{1}{n}})\right)^n \boldsymbol{1}_{( c\ ;\ +\infty)}(\theta)$

This test has the same significance level (by construction), but the power of the test is clearly lower than the one we got using the maximum of our sample, when $\theta\in\Theta_1$

C2=theta0*(1-alpha^(1/n))
P2=(1-(theta0/theta)*(1-alpha^(1/n)))^n*(theta>C2)
lines(theta,P2,type="l",lwd=2,col="red")

Why not consider a test based on $\overline{x}$? The problem is that we need the distribution (more specifically the survival distribution) of $\overline{X}$. We can compute it, numerically. But that might be painful. An alternative is to consider some approximation, based on the central limit theorem, i.e.

$2\overline{X}\sim\mathcal{N}\left(\theta,2^2 \frac{\theta^2}{12n}\right)$

Our test is based on $\psi(\boldsymbol{x})=\boldsymbol{1}_{(c\ ;\ +\infty)}(2\overline{x})$, and to get the same significance as before, use

$c=\Phi_{\mu,\sigma^2}^{-1}(1-\alpha)=\mu+\sigma \Phi^{-1}(1-\alpha)$

The power of the test is then

$p(\theta)=1-\Phi_{\theta,\sigma^2}(c)\cdot \boldsymbol{1}_{(c,+\infty)}(\theta)$

Here it is

mu=2*(theta0/2)
s2=2^2*(theta0^2/12)/n
C3=qnorm(1-alpha,mu,sqrt(s2))
(P=1-pnorm(C3,theta,sqrt(s2)))*(theta>C3)
lines(theta,P)

Observe here that the test based on the maximum is not more powerful than the one based on the average (I just wonder if it could be due to the Gaussian approximation…).

# Test, valeur critique et p-value

Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la $p$-value dans la lecture de la sortie d’un test.

• Les erreurs dans un test statistique

Mais avant, rappelons qu’un test est une prise de décision: accepter ou rejeter une hypothèse. Et qu’on peut commettre une erreur. Ou pour être plus précis, on peut commettre deux types d’erreur,
• accepter l’hypothèse alors que cette dernière est fausse
• rejeter l’hypothèse alors que cette dernière était vraie
Pour reprendre une terminologie plus médicale, un test de grossesse peut dire à une femme qu’elle n’est pas enceinte, alors qu’elle l’est; ou dire qu’elle l’est, alors qu’elle ne l’est pas (voir tous les exemples dans les exercices de probabilités de l’examen P de la SOA, ou le cours ACT2121).
Formellement, on a deux probabilités,
• la probabilité d’accepter à tort notre hypothèse (on parlera d’erreur de second espèce), $\beta$
• la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) $\alpha$
Dans un monde idéal on voudrait que les deux probabilités soient aussi petites que possibles… Mais c’est impossible, et le plus souvent, baisser une des probabilités se fait en augmentant l’autre. Les cas extrêmes étant
• avoir un test de grossesse qui déclare tout le monde enceinte: on ne rejette alors jamais à tort (on ne rejette jamais tout court en fait), mais on a un fort taux d’acceptation à tort,
• avoir un test de grossesse qui ne déclare personne enceinte: on n’accepte jamais à tort (car on n’accepte jamais) mais on a un fort taux de rejet à tort.
Bref, on a un arbitrage à faire entre deux types d’erreurs. Souvent, en pratique on va demander à contrôler l’erreur de première espèce (i.e. $\alpha$ de l’ordre de 5%), et on chercher a un test qui, à $\alpha$ donné, possède la plus faible erreur de première espèce. Voilà en gros pour la théorie: on se donne un seuil de significativité $\alpha$, qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse $H_0$ est vraie, l’alternative étant une hypothèse $H_1$.

 $H_0$ vraie $H_1$ vraie accepter $H_0$ OK erreur type 2 rejeter $H_0$ erreur type 1 OK
• La “valeur critique”

La notion de valeur critique a été introduite dans Neyman & Pearson (1928). Cette valeur dépend de la forme de l’hypothèse alternative, en particulier savoir si le test est bilatéral, unilatéral à gauche, ou unilatéral à droite. Pour un test donné, la valeur critique peut-être vue comme la valeur limite a partir de laquelle on pourra rejeter $H_0$ avec un seuil de significativité donné.

• La p-value

La $p$-value a été introduite dans Gibbons & Pratt (1975), meme si on peut retrouve l’idée beaucoup plus tôt, comme Pearson (1900), qui propose de calculer “the probability that the observed value of the chi-square statistic would be exceeded under the null hypothesis“. La $p$-value est la probabilité, sous $H_0$, d’obtenir une statistique aussi extrême (pour ne pas dire aussi grande) que la valeur observée sur l’échantillon. Aussi, pour un seuil de significativité $\alpha$ donné, on compare $p$ et $\alpha$, afin d’accepter, ou de rejeter $H_0$,
• si $p\leq\alpha$, on va rejeter l’hypothèse $H_0$ (en faveur de $H_1$)
• si $p>\alpha$, on va rejeter $H_1$ (en faveur de $H_0$).
On peut alors interpréter la $p$-value comme le plus petit seuil de significativité pour lequel l’hypothèse nulle est acceptée. Gibbons & Pratt (1975) reviennent longuement sur les interprétations, et surtout les mauvaises interprétations, de cette $p$-value.

• Valeur critique versus p-value

Si on formalise un peu, on peut vouloir tester $H_0:\theta=\theta_0$ contre $H_1:\theta>theta_0$ (par exemple). De manière très générale, on dispose d’une statistique de test $T$ qui a pour loi, sous $H_0$, $F_{\theta_0}(\cdot)$ (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme $H_1:\theta\neq\theta_0$, c’est juste plus pénible parce qu’il faut travailler sur $\vert T\vert$, et calculer des probabilités à gauche, ou à droite. Donc pour notre exemple, on va prendre un test unilatéral.
Dans l’approche classique (telle que présentée dans tous les cours de statistiques), on se donne un seul d’acceptation $\alpha$ petit (disons 5%), et on cherche une valeur critique $T_{1-alpha}$ telle que

$\mathbb{P}(T\geq T_{1-\alpha}\vert\theta=\theta_0)=\alpha$

Pour ceux qui se souviennent de leur cours de stats, cela peut faire penser à la puissance du test, définie par

$\pi(\theta\vert \alpha)=\mathbb{P}(T\geq T_{1-\alpha}\vert \theta)=1-F_{\theta}(T_{1-\alpha})$

Formellement, la $p$-value associée au test $T$ est la variable aléatoire $P$ définie par
$P=1-F_{\theta_0}(T)$.
Donc effectivement, la $p$-value et la puissance sont liées, puisque

$\mathbb{P}(P\leq \alpha\vert \theta)=\pi(\theta\vert \alpha)$

autrement dit, la puissance peut-être vue comme la fonction de répartition de la $p$-value.

• Intérêt computationnel de la p-value

D’un point de vue computationnel, la $p$-value est l’outil le plus important pour interpréter la sortie d’un test. Commençons par un test simple, comme une comparaison de moyennes. On cherche ici à tester $H_0:\mu_X=\mu_Y$ contre $H_1:\mu_X>\mu_Y$ pour des moyennes calculées sur deux groupes. Pour reprendre l’exemple abordé dans un précédant billet, on a les notes obtenues en ACT6420 par deux groupes différents. Et on veut savoir s’ils sont vraiment différents (ci-dessous le nombre de bonnes réponses, sur 40 questions, on travaillera ensuite sur la note sur 100)

image manquante

La statistique de test est ici

$T = \frac{\overline{X} - \overline{Y}}{\displaystyle{ \sqrt{ {s_X^2 \over n_X} + {s_Y^2 \over n_Y} }}}$

et sous $H_0$, $T$ va suivre une loi de Student à $\nu$ degrés de liberté, où $\nu$ est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)),

$\nu = {{\left( {s_X^2 \over n_X} + {s_Y^2 \over n_Y}\right)^2 } \over {{s_X^4 \over n_X^2 \cdot \left({n_X-1}\right)}+{s_Y^4 \over n_Y^2 \cdot \left({n_Y-1}\right)}}}$

Numériquement, on a ici

> Xbar=mean(X)
> Ybar=mean(Y)
> Sx2=var(X)
> Sy2=var(Y)
> nX=length(X)
> nY=length(Y)
> (T=(Xbar-Ybar)/sqrt(Sx2/nX+Sy2/nY))
[1] -2.155754

et pour les degrés de liberté

> (nu=(Sx2/nX+Sy2/nY)^2/(Sx2^2/nX^2/(nX-1)+
+ Sy2^2/nY^2/(nY-1)))
[1] 36.35279

La valeur critique est obtenue en lisant dans les tables,

(car ici on a des probabilité pour un test bilatéral dans la table) comme on apprenait dans les cours de statistique au siècle passé. D’un point de vue informatique, on cherche à savoir si on est à gauche, ou à droite de la valeur critique

> qt(.05,df=nu)
[1] -1.687865

image manquante

On peut aussi calculer la $p$-value,

> pt(T,df=nu)
[1] 0.01889768

Si on regarde, sous R, il existe des fonctions de tests, pour comparer des moyennes. Et dans ce cas, la sortie est

> t.test(X,Y,alternative = "less")

Welch Two Sample t-test

data:  X and Y
t = -2.1558, df = 36.353, p-value = 0.0189
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf -1.772507
sample estimates:
mean of x mean of y
48.75000  56.91667

Autrement dit, on a automatiquement la $p$-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une $p$-value (et que l’on vérifié au préalable les conditions d’application d’un test), on peut faire tous les tests que l’on veut !
Si on veut faire un peu plus compliqué, on peut regarder la distribution des notes, et se demander si une loi $\mathcal{N}(60,15^2)$ serait possible (par exemple, ça sera notre hypothèse $H_0$, l’hypothèse alternative étant que ce n’est pas cette loi). Pour faire ce test, il existe le test de Kolmogorov-Smirnov. La statistique de test est ici

$T=\sup\{\vert \widehat{F}_n(x)-F_0(x)\vert ,x\in\mathbb{R}\}$

$F_0(\cdot)$ est la fonction de répartition de la loi $\mathcal{N}(60,15^2)$, et $\widehat{F}_n(\cdot)$ est la fonction de répartition empirique

$\widehat{F}_n(x)=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(x_i\leq x)$

La loi de $T$ n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les $p$-values automatiquement,

> ks.test(Y, "pnorm", 60, 15)

One-sample Kolmogorov-Smirnov test

data:  Y
D = 0.1421, p-value = 0.5796
alternative hypothesis: two-sided

Aussi, on peut accepter ici l’hypothèse nulle. On peut d’ailleurs faire un petit dessin pour s’en convaincre,

> Femp=function(x) mean(Y<=x)
> plot(0:100,Vectorize(Femp)(0:100),type="s")
> lines(0:100,pnorm(0:100,60,15),col="red")

image manquante

Et ça va nous servir dans ce cours ? A priori oui… parce qu’on parlera du test de Student (pour tester si une variable dans une régression est significative), du test de Fisher (pour tester si plusieurs variables dans une régression sont significatives, ou plus généralement si une contrainte – linéaire – sur les coefficients peut être acceptée), du test de Chow (pour tester des ruptures dans un modèle linéaire, mais c’est un test de Fisher un peu déguisé), du test d’Anderson-Darling (pour tester si des résidus sont Gaussiens), du test de Breuch-Pagan voire le test de White (pour tester si les résidus peuvent être considérés de variance constante), du test de Durbin-Watson (pour tester s’il n’y a pas d’auto-corrélation dans la série des résidus), du test de Dickey-Fuller (pour tester si une série temporelle est – ou n’est pas – stationnaire), des tests de Franses (pour tester si une série peut être considérée comme saisonnière, ou pas), du test de Ljung-Box (pour tester si un bruit est un bruit blanc)… Et j’en oublie un paquet. Donc quand il est dit (dans le plan de cours) que le cours de statistique est un prérequis, il ne s’agit pas de l’avoir suivi, mais bel et bien de l’avoir compris, car on passera notre temps à utiliser des notions entrevues dans ce cours.

# (nonparametric) copula density estimation

Today, we will go further on the inference of copula functions. Some codes (and references) can be found on a previous post, on nonparametric estimators of copula densities (among other related things).  Consider (as before) the loss-ALAE dataset (since we’ve been working a lot on that dataset)

> library(MASS)
> library(evd)
> X=lossalae
> U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1))

The standard tool to plot nonparametric estimators of densities is to use multivariate kernels. We can look at the density using

> mat1=kde2d(U[,1],U[,2],n=35)
> persp(mat1$x,mat1$y,mat1$z,col="green", + shade=TRUE,theta=s*5, + xlab="",ylab="",zlab="",zlim=c(0,7)) or level curves (isodensity curves) with more detailed estimators (on grids with shorter steps) > mat1=kde2d(U[,1],U[,2],n=101) > image(mat1$x,mat1$y,mat1$z,col=
+ rev(heat.colors(100)),xlab="",ylab="")
> contour(mat1$x,mat1$y,mat1$z,add= + TRUE,levels = pretty(c(0,4), 11)) Kernels are nice, but we clearly observe some border bias, extremely strong in corners (the estimator is 1/4th of what it should be, see another post for more details). Instead of working on sample $(U_i,V_i)$ on the unit square, consider some transformed sample $(Q(U_i),Q(V_i))$, where $Q:(0,1)\rightarrow\mathbb{R}$ is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the $\mathcal{N}(0,1)$ distribution. Then, we can estimate the density of the transformed sample, and using the inversion technique, derive an estimator of the density of the initial sample. Since the inverse of a (general) function is not that simple to compute, the code might be a bit slow. But it does work, > gaussian.kernel.copula.surface <- function (u,v,n) { + s=seq(1/(n+1), length=n, by=1/(n+1)) + mat=matrix(NA,nrow = n, ncol = n) + sur=kde2d(qnorm(u),qnorm(v),n=1000, + lims = c(-4, 4, -4, 4)) + su<-sur$z
+ for (i in 1:n) {
+     for (j in 1:n) {
+ 	Xi<-round((qnorm(s[i])+4)*1000/8)+1;
+ 	Yj<-round((qnorm(s[j])+4)*1000/8)+1
+ 	mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))*
+ 	dnorm(qnorm(s[j])))
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

Here, we get

Note that it is possible to consider another transformation, e.g. the quantile function of a Student-t distribution.

> student.kernel.copula.surface =
+  function (u,v,n,d=4) {
+  s <- seq(1/(n+1), length=n, by=1/(n+1))
+  mat <- matrix(NA,nrow = n, ncol = n)
+ sur<-kde2d(qt(u,df=d),qt(v,df=d),n=5000,
+ lims = c(-8, 8, -8, 8))
+ su<-sur$z + for (i in 1:n) { + for (j in 1:n) { + Xi<-round((qt(s[i],df=d)+8)*5000/16)+1; + Yj<-round((qt(s[j],df=d)+8)*5000/16)+1 + mat[i,j]<-su[Xi,Yj]/(dt(qt(s[i],df=d),df=d)* + dt(qt(s[j],df=d),df=d)) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + } Another strategy is to consider kernel that have precisely the unit interval as support. The idea is here to consider the product of Beta kernels, where parameters depend on the location > beta.kernel.copula.surface= + function (u,v,bx=.025,by=.025,n) { + s <- seq(1/(n+1), length=n, by=1/(n+1)) + mat <- matrix(0,nrow = n, ncol = n) + for (i in 1:n) { + a <- s[i] + for (j in 1:n) { + b <- s[j] + mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) * + dbeta(b,v/by,(1-v)/by)) / length(u) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + } On those two graphs, we can clearly observe strong tail dependence in the upper (right) corner, that cannot be intuited using a standard kernel estimator… # Sondages et probabilités en trois graphiques Petite précision liminaire: je suis nul en sondage. Je n’ai jamais rien compris aux sondages. C’est la matière où j’ai eu ma pire note quand j’étais étudiant. Mais si on voit un sondage comme l’estimation d’une probabilité de succès, dans des tirages de lois binomiales (ou multinomiales), alors je me sens plus à l’aise pour commenter des choses lues ici ou là. Car en ces temps d’élections (qui approchent), on lit beaucoup de choses sur les sondages ou sur les intervalles de confiance des sondages, voire sur des comparaison de sondages… mais rares sont les justifications. Donc on va creuser un peu (histoire de voir ce qui peut être dit, dans un cadre largement simplifié). • “Il y a 1 chance sur 4 qu’un sondage donne une variation d’un point en + ou en – sur une proportion inchangée” comme le disait @adelaigue (sur Twitter). Faire un sondage politique (si on résume simplement) c’est estimer la probabilité de voter pour A, sachant que l’on a alors une probabilité pour qu’une personne prise au hasard ne vote par pour A. Sur personnes interrogées, le nombre de personnes qui déclarera voter pour A est , qui suit une loi binomiale de paramètres et . Or par approximation (ou convergence) de la loi binomiale par une loi Gaussienne, on peut dire que Or l’estimateur naturel de est , i.e. Aussi, on en déduit l’intervalle de confiance classique, tel que appartient à cet intervalle avec une probabilité , désigne la fonction quantile de la loi normale centrée réduite; ou encore, en remplaçant les quantités inconnues par des estimateurs, Les bornes de l’intervalle de confiance de niveau sont alors C’est ce qu’on peut visualiser sur la figure ci-dessous, Comme le notait @adelaigue, si on souhaite avoir (car on parle d’un point de base 100), il faudrait que Classiquement, est aux alentours de 1000. Et ce qu’on cherche c’est une probabilité d’appartenir à un intervalle, i.e. , Graphiquement, en fonction de on a alors p=seq(0,1,by=.0025) plot(p,2*(1-pnorm(.01/sqrt(p*(1-p)/1000)))) (ca sera notre premier graphique). Effectivement, entre 35% et 65%, il y a environ 50% de chances d’être autour de la fréquence observée, à un point près. Ou 50% de chances d’avoir plus d’un point de différence. Autrement dit, si on a observé avec 1000 personnes une fréquence de 40%, il y a une chance sur 4 que la probabilité soit entre 39% et 40%, une chance sur 4 qu’elle soit entre 40% et 41%, et surtout une chance sur 4 qu’elle soit inférieure à 39%, et 1 chance sur 4 qu’elle soit supérieure à 41%. Comme le notait @adelaigue, sur les bords, il semble que la formule ne soit pas valide. Mais si on prend justement les cas limites (probabilité très proche de 0, voire franchement nulle), est-ce qu’on peut dire quelque-chose ? • peut-on dire quelque chose (malgré tout) quand on ne peut rien dire ? ou plus précisément, si en interrogeant personnes et que personne n’a déclaré vouloir voter pour A, quelle serait l’intervalle de confiance de la probabilité de voter pour A ? J’avais évoqué ce point dans un très vieux billet. La première piste est la réponse apporté par Pearson-Clopper (Pearson & Clopper (1934)), basée sur l’utilisation de la loi binomiale, et non plus d’une approximation. library(Hmisc) bornesup=function(N){ binom.test(x=0, n=N, p = 0, alternative="less")$conf.int[2]*100 }
BS=Vectorize(bornesup)
vn=seq(10,500,by=10)
plot(vn,BS(vn))

Une alternative est de faire du bayésien. On utilise le fait que la loi conjuguée de la loi binomiale est la loi beta (comme évoqué dans des précédents billets). Si on suppose que la loi a priori de est une loi alors la loi a posteriori, sachant que sur lancés, on n’a observé que des 0, devient

bornesup=function(N){ qbeta(.95,1,1+N)*100 }
BS=Vectorize(bornesup)
lines(vn,BS(vn),lwd=2,col="blue")

La courbe ci-dessus correspond au calcul dans le cas du modèle de Pearson-Clopper, en rouge, et avec le quantile de la loi Beta en bleu (qui sont confondues sur ce graphique). Autrement dit, avec 300 personnes interrogées, si aucune ne dit vouloir voter pour A, l’intervalle de confiance pour est [0%;1%] (environ). Ca sera notre second graphique.

• si on a trois alternatives, quelles sont les vraies probabilités de gagner ?

c’était déjà sur cette question que je concluais mon précédant billet, mais il est possible d’exploiter davantage cette idée de modèle bayésien avec un a priori suivant une loi beta (comme au dessus) ou de Dirichlet (si on a plus de deux alternatives). La première chose que l’on peut noter, c’est que s’il y a trois alternatives (A, B et C), et que la probabilité de voter pour A est 40%, la probabilité que A gagne… dépend des autres probabilités. On ne peut pas faire comme si A était seul, contre tous les autres:
– si les probabilités de voter pour B et C sont de l’ordre de 30%, on se dit que A a de fortes chances de gagner (oui, pour l’instant ce n’est pas très quantitatif).
– si les probabilités de voter pour B et C sont respectivement 55% et 5%, alors on se dit que A a de faibles chances de gagner.
Plus précisément, si on souhaite savoir si A va l’emporter ou pas, et que 200 personnes ont été interrogées, on peut essayer de faire une abaque, avec des courbes de niveau… dans le simplexe. Pour rappel, les scores obtenus par A, B et C se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%))

Les probabilités sont les suivantes, avec n=200 (personnes qui se seraient exprimés)

où la région rouge correspond à une (quasi) certitude de gagner, pour A, et la région blanche, au contraire, une probabilité infinitésimale de gagner.

Ces résultats ont été obtenus par simulation,
vx=v[i]
vz=v[j]
vy=1-vz-vx
alpha=c(vx,vy,vz)
alpha=c(alpha,1-sum(alpha))*200
RD=rdirichlet(ns, alpha)
PROBA[i,j]=mean(RD[,1]>=apply(
cbind(RD[,3],RD[,2]),1,max))

Pour ceux que ne sont pas très à l’aise avec les excursions dans le simplexe, on peut se ramener dans un repère orthonormé. Ici, les scores obtenus par A et B se lisent sur le triangle ci-dessous (le point représenté correspond au triplet (40%,35%,25%), le même qu’auparavant

On a alors les probabilités de gagner suivantes, avec à droite la représentation sur le simplexe, en zoomant autour de notre triplet,

(ça sera le troisième et dernier graph). Si on regarde plus en détails, l’intervalle de confiance sur la probabilité de gagner dépend de la région où on se trouve dans le simplexe (un peu comme l’intervalle de confiance dépendait de la probabilité de l’emporter, dans la première partie de ce billet).

Maintenant, comme le rappelait @Christian dans un commentaire, tout cela est la version très simpliste de la pratique des sondages. Dans les vrais sondages (et pas des exercices de probabilités comme on vient d’en faire ici), on utilise des quotas, on fait du redessement, et c’est beaucoup plus complexe (enfin… c’est ce que j’essaye de me dire, pour me convaincre que la mauvaise note que j’avais obtenue lorsque j’étais étudiant était motivée).

# the Dirichlet distribution

In the course, since we are still introducing some concepts of dependent distributions, we will talk about the Dirichlet distribution, which is a distribution over the simplex of . Let denote the Gamma distribution with density (on )

Let denote independent random variables, with . Then where

has a Dirichlet distribution with parameter

Note that has a distribution in the simplex of ,

and has density

We will write .

The density for different values of can be visualized below, e.g. , with some kind of symmetry,

or and , below

and finally, below,

Note that marginal distributions are also Dirichlet, in the sense that if

then

if , and if , then ‘s have Beta distributions,

See Devroye (1986) section XI.4, or Frigyik, Kapila & Gupta (2010) .This distribution might also be called multivariate Beta distribution. In R, this function can be used as follows

> library(MCMCpack)
> alpha=c(2,2,5)
> x=seq(0,1,by=.05)
> vx=rep(x,length(x))
> vy=rep(x,each=length(x))
> vz=1-x-vy
> V=cbind(vx,vy,vz)
> D=ddirichlet(V, alpha)
> persp(x,x,matrix(D,length(x),length(x))

(to plot the density, as figures above). Note that we will come back on that distribution later on so-called Liouville copulas (see also Gupta & Richards (1986)).

# Exchangeability, credit risk and risk measures

Exchangeability is an extremely concept, since (most of the time) analytical expressions can be derived. But it can also be used to observe some unexpected behaviors, that we will discuss later on with a more general setting. For instance, in a old post, I discussed connexions between correlation and risk measures (using simulations to illustrate, but in the context of exchangeable risk, calculations can be performed more accurately). Consider again the standard credit risk problem, where the quantity of interest is the number of defaults in a portfolio. Consider an homogeneous portfolio of exchangeable risk. The quantity of interest is here

or perhaps the quantile function of the sum (since the Value-at-Risk is the standard risk measure). We have seen yesterday that – given the latent factor – (either the company defaults, or not), so that

i.e. we can derive the (unconditional) distribution of the sum

so that the probability function of the sum is, assuming that

Thus, the following code can be used to calculate the quantile function

> proba=function(s,a,m,n){
+ b=a/m-a
+ choose(n,s)*integrate(function(t){t^s*(1-t)^(n-s)*
+ dbeta(t,a,b)},lower=0,upper=1,subdivisions=1000,
+ stop.on.error =  FALSE)$value + } > QUANTILE=function(p=.99,a=2,m=.1,n=500){ + V=rep(NA,n+1) + for(i in 0:n){ + V[i+1]=proba(i,a,m,n)} + V=V/sum(V) + return(min(which(cumsum(V)>p))) } Now observe that since variates are exchangeable, it is possible to calculate explicitly correlations of defaults. Here i.e. Thus, the correlation between two default indicators is then Under the assumption that the latent factor is beta distributed we get Thus, as a function of the parameter of the beta distribution (we consider beta distributions with the same mean, i.e. the same margin distributions, so we have only one parameter left, with is simply the correlation of default indicators), it is possible to plot the quantile function, > PICTURE=function(P){ + A=seq(.01,2,by=.01) + VQ=matrix(NA,length(A),5) + for(i in 1:length(A)){ + VQ[i,1]=QUANTILE(a=A[i],p=.9,m=P) + VQ[i,2]=QUANTILE(a=A[i],p=.95,m=P) + VQ[i,3]=QUANTILE(a=A[i],p=.975,m=P) + VQ[i,4]=QUANTILE(a=A[i],p=.99,m=P) + VQ[i,5]=QUANTILE(a=A[i],p=.995,m=P) + } + plot(A,VQ[,5],type="s",col="red",ylim= + c(0,max(VQ)),xlab="",ylab="") + lines(A,VQ[,4],type="s",col="blue") + lines(A,VQ[,3],type="s",col="black") + lines(A,VQ[,2],type="s",col="blue",lty=2) + lines(A,VQ[,1],type="s",col="red",lty=2) + lines(A,rep(500*P,length(A)),col="grey") + legend(3,max(VQ),c("quantile 99.5%","quantile 99%", + "quantile 97.5%","quantile 95%","quantile 90%","mean"), + col=c("red","blue","black", +"blue","red","grey"), + lty=c(1,1,1,2,2,1),border=n) +} e.g. with a (marginal) default probability of 15%, > PICTURE(.15) On this graph, we observe that the stronger the correlation (the more on the left), the higher the quantile… Note that the same graph can be plotted with on the X-axis the correlation, Which is quite intuitive, somehow. But if the marginal probability of default decreases, increasing the correlation might decrease the risk (i.e. the quantile function), > PICTURE(.05) (with the modified code to visualize the quantile as a function of the underlying default correlation) or even worse, > PICTURE(.0075) And it because all the more counterintuitive that the default probability decreases ! So in the case of a portfolio of non-very risky bond issuers (with high ratings), assuming a very strong correlation will lower risk based capital ! # Beta kernel and transformed kernel This Thursday I will give a talk at Laval University, on “Beta kernel and transformed kernel : applications to copula density estimation and quantile estimation“. This time, I will talk at the department of Mathematics and Statistics (13:30 at the pavillon Adrien-Pouliot). “Because copulas have bounded support (the unit square in dimension 2), standard kernel based estimators of densities are (multiplicatively) biased on borders and in corners of the support. Two techniques can be used to avoid that underestimation: Beta kernels and Transformed kernel. We will describe and discuss those two techniques in the first part of the talk. Then, we will see that it is possible to combine those two techniques to get nice estimator of several quantities (e.g. quantiles): transform the data to get on the unit interval – using a transformed kernel – then estimate the (transformed) quantile on [0,1] using a beta kernel, then get back on the initial support. As we will see on simulations, that technique can be better than standard quantile estimators, especially when data are heavy tailed.” Slides can be downloaded here. • kernel based density estimation Kernel based estimation are a popular (and natural) technique to estimate densities. It is simply and extension of the moving histogram: so we count how many observations are a the neighborhood of the point where we want to estimate the density of the distribution. Then it is natural so consider a smoothing function, i.e. instead of a step function (either observations are close enough, or not), it is possible to give weights to observations, which will be a decreasing function of the distance, With a smooth kernel, we have a smooth estimation of the density Then it is possible to play on the bandwidth, either to get a more accurate estimation of the density, but not that smooth (small bias but large variance), or a smoother one (large bias, but small variance), In R, it is simply > X=rnorm(100) > (D=density(X)) Call: density.default(x = X) Data: X (100 obs.); Bandwidth 'bw' = 0.3548 x y Min. :-3.910799 Min. :0.0001265 1st Qu.:-1.959098 1st Qu.:0.0108900 Median :-0.007397 Median :0.0513358 Mean :-0.007397 Mean :0.1279645 3rd Qu.: 1.944303 3rd Qu.:0.2641952 Max. : 3.896004 Max. :0.3828215 > plot(D$x,D$y) • Beta kernel The idea of Beta kernel is to consider kernels having support [0,1]. In the univariate case, where is the density of a Beta distribution, i.e. For additional material, I have uploaded some R code to fit copula densities using beta kernels, library(copula) beta.kernel.copula.surface = function (u,v,bx,by,p) { s = seq(1/p, len=(p-1), by=1/p) mat = matrix(0,nrow = p-1, ncol = p-1) for (i in 1:(p-1)) { a = s[i] for (j in 1:(p-1)) { b = s[j] mat[i,j] = sum(dbeta(a,u/bx,(1-u)/bx) * dbeta(b,v/by,(1-v)/by)) / length(u) } } return(data.matrix(mat)) } Then we can used it to see what we get on a simulated sample library(copula) COPULA = frankCopula(param=5, dim = 2) X = rcopula(n=1000,COPULA) p0 = 26 Z= beta.kernel.copula.surface(X[,1],X[,2],bx=.01,by=.01,p=p0) u = seq(1/p0, len=(p0-1), by=1/p0) persp(u,u,Z,theta=30,col="green",shade=TRUE, box=FALSE,zlim=c(0,6)) (yes, the surface is changing… to illustrate the impact of the bandwidth on the estimation). • transformed kernel estimation I the talk, I will also mention the transformed Kernel estimate, as introduced in the book on L1 density estimation by Luc Devroye and Laszlo Györfi (the book can be downloaded here). I probably spend a few minutes on the original chapter, in order to provide another application of that techniques (not only to estimate copula densities, but here to estimate quantiles of heavy tailed distribution). In the univariate case, the R code is the following (here I consider two transformation, the quantile function of the Gaussian distribution, and the quantile function of the Student distribution with 3 degrees of freedom), set.seed(1) sample=rbeta(100,4,3) transfN = function(x){ Y=qnorm(sample) f=density(Y,from=-4,to=4,n=2001) ny=sum(f$x<=qnorm(x));
g=f$y[ny]/dnorm(qnorm(x)) return(g) } df0=3 transfT = function(x){ Y=qt(sample,df=df0) f=density(Y,from=-4,to=4,n=2001) ny=sum(f$x<=qt(x,3));
g=f\$y[ny]/dt(qt(x,df=df0),df=df0)
return(g)
}

tN=Vectorize(transfN)
tT=Vectorize(transfT)

u=seq(.01,.99,by=.01)
vN=tN(u)
vT=tT(u)
plot(u,vN,type="l",lwd=3,col="blue")
lines(u,vT,lwd=3,col="green")
lines(u,dbeta(u,4,3),col="red",lty=2)

The density estimation is the following,

(the red dotted line is the true density, since we work on a simulated sample). Now, let us get back on the initial chapter,

In the book, this is introduced as follows,

The original idea we add it to use this kernel based estimator for copulas, i.e. since we can estimate densities in high dimension with unbounded support, using

the idea is to transform marginal observations,

and to use the fact that the associated copula density can be written

to derive an intuitive estimator for the copula density

An important issue is how do we choose the transformation

And Luc Devroye and Laszlo Györfi mention that this can be used to deal with extremes.

well, extremes are introduced through bumps (which is not the way I would have been dealing with extremes)

and note that several results can be derived on those bumps,

e.g.

Then, there is an interesting discussion about estimating the optimal transformation

and I will prove that this can be an extremely interesting idea, for instance to estimate quantiles of heavy tailed distribution, if we use also the beta kernel estimator on the unit interval. This idea was developed in a paper with Abder Oulidi, online here.

Remark: actually, in the book, an additional reference is mentioned,

but I have never been able to find a copy… if anyone has one, I’d be glad to read it…

# Will I ever be a bayesian statistician ? (part 1)

Last week, during the workshop on Statistical Methods for Meteorology and Climate Change (here), I discovered how powerful bayesian techniques could be, and that there were more and more bayesian statisticians. So, if I was to fully understand applied statisticians in conferences and workshops, I really have to understand basics of bayesian statistics. I have published some time ago some posts on bayesian statistics applied to actuarial problems (here or there), but so far, I always thought that bayesian was a synonym for magician. To be honest, I am a Muggle, and I have not been trained as a bayesian. But I can be an opportunist…

So I decided to publish some posts on bayesian techniques, in order to prove that it is actually not that difficult to implement.

As far as I understand it, in bayesian statistics, the parameter is considered as a random variable (which is also the case, in classical mathematical statistics). But here, here assume that this parameter does have a parametric distribution….
Consider a classical statistical problem: assume we have a sample i.i.d. with distribution . Here we note

since parameter  is a random variable. The idea is to assume that has a (so called a priori) distribution, e.g.

So far it was simple. The idea is then to consider the posterior distribution of , given the observations . Thus, we need to compute the distribution of which is here extremely simple (due to properties of the Gaussian distribution), i.e.

where

And them, it becomes extremely natural to consider as an estimator of given our sample data (and thus, we also have a confidence interval since we know the distribution of given the observations ).
In order to be sure that we understood, consider now a heads and tails problem, i.e. . Note, first, that \theta has support . So we need a distribution on that support. Why not a beta distribution ? E.g.

Thus,

and

From Bayes formula,

and we get easily

which is the density of a Beta distribution, i.e.

prior=dbeta(u,a,b)
posterior=dbeta(u,a+y,n-y+b)

The estimator proposed is then the expected value of that conditional distribution,

Note that

Further, it is possible to derive confidence intervals using quantiles of the posterior distribution.
On the graphs below, we consider the following heads/tails sample

A first idea is to consider a uniform prior distribution.

A second idea is to consider an asymmetric beta distribution. First, with an asymmetry on the left,

or on the right

Finally a third idea is simply to get back to the standard Gaussian approximation,

If we compare the four models, we obtain (the plain black line is the Gaussian approximated distribution for the empirical mean), and red lines are obtained from prior beta distributions

The code to generate those graphs is the following
a1=1; b1=1
D1[1,]=dbeta(u,a,b)
a2=4; b2=2
D2[1,]=dbeta(u,a,b)
a3=2; b3=4
D3[1,]=dbeta(u,a,b)
setseed(1)
S=sample(0:1,size=100,replace=TRUE)
COULEUR=rev(rainbow(120))
D1=D2=D3=D4=matrix(NA,101,length(u))
for(s in 1:100){
y=sum(S[1:s])
D1[s+1,]=dbeta(u,a1+y,s-y+b1)
D2[s+1,]=dbeta(u,a2+y,s-y+b2)
D3[s+1,]=dbeta(u,a3+y,s-y+b3)
D4[s+1,]=dnorm(u,y/s,sqrt(y/s*(1-y/s)/s))
plot(u,D1[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D1[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D2[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D2[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D3[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D3[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D4[1,],col="white",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D4[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D4[s+1,],col="black",lwd=2,type="l",
ylim=c(0,8),xlab="",ylab="")
lines(u,D1[1+i,],col="blue")
lines(u,D2[1+i,],col="red")
lines(u,D3[1+i,],col="purple")
points(y/s,0,pch=3,cex=2)
}

Here, we can see that computations are simple if the prior distribution has a distribution which is the conjugate of the observations’ distribution (see here for the list of prior and posterior standard distributions).
So far, I have two questions that naturally show up

• is it possible to start with a neutral prior distribution, non informative ?
• what if we are no longer working with conjugate distributions ?

Well, I guess I have to work a bit more to answer those questions…. to be continued

# Estimation de quantile par noyau beta

Le papier sur l’estimation de quantile par noyau beta, coécrit avec Abder Oulidi, est accepté pour publication dans Statistics and Computing, http://link.springer.com/…

In this paper we propose several nonparametric estimators of quantiles based on Beta kernel and applied to transformed data by the generalized Champernowne distribution initially fitted to the data. A Monte-Carlo based study will show that those estimators improve the efficiency of a traditional ones, not only for light tailed distributions, but also heavy tails, when the probability level is close to 1. We also compare these estimators with the Extreme Value Theory Quantile applying to Danish data on large fire insurance losses.