Category Archives: Bayesian

Confidence vs. Credibility Intervals

Tomorrow, for the final lecture of the Mathematical Statistics course, I will try to illustrate – using Monte Carlo simulations – the difference between classical statistics, and the Bayesien approach.

The (simple) way I see it is the following,

  • for frequentists, a probability is a measure of the the frequency of repeated events, so the interpretation is that parameters are fixed (but unknown), and data are random
  • for Bayesians, a probability is a measure of the degree of certainty about values, so the interpretation is that parameters are random and data are fixed

Or to quote Frequentism and Bayesianism: A Python-driven Primer,  a Bayesian statistician would say “given our observed data, there is a 95% probability that the true value of \theta falls within the credible region” while a Frequentist statistician would say “there is a 95% probability that when I compute a confidence interval from data of this sort, the true value of \theta will fall within it”.

To get more intuition about those quotes, consider a simple problem, with Bernoulli trials, with insurance claims. We want to derive some confidence interval for the probability to claim a loss. There were https://latex.codecogs.com/gif.latex?n = 1047 policies. And 159 claims.

Consider the standard (frequentist) confidence interval. What does that mean that

\overline{x}\pm%201.96%20\sqrt{\frac{\overline{x}(1-\overline{x})}{n}}

is the (asymptotic) 95% confidence interval? The way I see it is very simple. Let us generate some samples, of size n, with the same probability as the empirical one, i.e. \widehat{\theta} (which is the meaning of “from data of this sort”). For each sample, compute the confidence interval with the relationship above. It is a 95% confidence interval because in 95% of the scenarios, the empirical value lies in the confidence interval. From a computation point of view, it is the following idea,

> xbar <- 159
> n <- 1047
> ns <- 100
> M=matrix(rbinom(n*ns,size=1,prob=xbar/n),nrow=n)

I generate 100 samples of size https://latex.codecogs.com/gif.latex?n. For each sample, I compute the mean, and the confidence interval, from the previous relationship

> fIC=function(x) mean(x)+c(-1,1)*1.96*sqrt(mean(x)*(1-mean(x)))/sqrt(n)
> IC=t(apply(M,2,fIC))
> MN=apply(M,2,mean)

Then we plot all those confidence intervals. In red when they do not contain the empirical mean

> k=(xbar/n<IC[,1])|(xbar/n>IC[,2])
> plot(MN,1:ns,xlim=range(IC),axes=FALSE,
+ xlab="",ylab="",pch=19,cex=.7,
+ col=c("blue","red")[1+k])
> axis(1)
> segments(IC[,1],1:ns,IC[,2],1:
+ ns,col=c("blue","red")[1+k])
> abline(v=xbar/n)

Now, what about the Bayesian credible interval ? Assume that the prior distribution for the probability to claim a loss has a https://latex.codecogs.com/gif.latex?\mathcal{B}(\alpha,\beta) distribution. We’ve seen in the course that, since the Beta distribution is the conjugate of the Bernoulli one, the posterior distribution will also be Beta. More precisely

https://latex.codecogs.com/gif.latex?\mathcal{B}\left(\alpha+\sum%20x_i,\beta+n-\sum%20x_i\right)

Based on that property, the confidence interval is based on quantiles of that (posterior) distribution

> u=seq(.1,.2,length=501)
> v=dbeta(u,1+xbar,1+n-xbar)
> plot(u,v,axes=FALSE,type="l")
> I=u<qbeta(.025,1+xbar,1+n-xbar)
> polygon(c(u[I],rev(u[I])),c(v[I],
+ rep(0,sum(I))),col="red",density=30,border=NA)
> I=u>qbeta(.975,1+xbar,1+n-xbar)
> polygon(c(u[I],rev(u[I])),c(v[I],
+ rep(0,sum(I))),col="red",density=30,border=NA)
> axis(1)

What does that mean, here, that we have a 95% credible interval. Well, this time, we do not draw using the empirical mean, but some possible probability, based on that posterior distribution (given the observations)

> pk <- rbeta(ns,1+xbar,1+n-xbar)

In green, below, we can visualize the histogram of those values

> hist(pk,prob=TRUE,col="light green",
+ border="white",axes=FALSE,
+ main="",xlab="",ylab="",lwd=3,xlim=c(.12,.18))

And here again, let us generate samples, and compute the empirical probabilities,

> M=matrix(rbinom(n*ns,size=1,prob=rep(pk,
+ each=n)),nrow=n)
> MN=apply(M,2,mean)

Here, there is 95% chance that those empirical means lie in the credible interval, defined using quantiles of the posterior distribution. We can actually visualize all those means : in black the mean used to generate the sample, and then, in blue or red, the averages obtained on those simulated samples,

> abline(v=qbeta(c(.025,.975),1+xbar,1+
+ n-xbar),col="red",lty=2)
> points(pk,seq(1,40,length=ns),pch=19,cex=.7)
> k=(MN<qbeta(.025,1+xbar,1+n-xbar))|
+ (MN>qbeta(.975,1+xbar,1+n-xbar))
> points(MN,seq(1,40,length=ns),
+ pch=19,cex=.7,col=c("blue","red")[1+k])
> segments(MN,seq(1,40,length=ns),
+ pk,seq(1,40,length=ns),col="grey")

More details and exemple on Bayesian statistics, seen with the eyes of a (probably) not Bayesian statistician in my slides, from my talk in London, last Summer,

Bayesian Wizardry for Muggles

Monday, I will be giving the closing talk of the R in Insurance Conference, in London, on Bayesian Computations for Actuaries, or to be more specific, Getting into Bayesian Wizardry… (with the eyes of a muggle actuary). The animated version of the slides (since we will spend some time on MCMC algorithm, I thought that animated graphs could be more informative) can be downloaded from here.

Those slides are based on the chapter writen with Ben Escoto for the Computational Actuarial Science with R book, and some previous work.

Bayes, 250 years ago

Exactly 250 years ago, on December 23rd, while some people were still shopping to get presents for their kids, Mr Price took some time to write a letter to the Royal Society of London,

He sent an essay towards solving a problem in the doctrine of chances, writen by an old friend of his (who just died), Mr. Bayes,

A few days ago, Stephen Stigler came in Montréal to give a talk on Bayes, and I learnt that this title is actually not the one Price chose for the essay. The true title of the essay sent to the Royal Academy was

i.e.

The story of the true title of Bayes’s essay can be found on http://arxiv.org/, There are interesting connexions between Bayes theory and Christianism, actually. But I won’t spend time on December 23rd, talking about probabilities… I still have presents to find…

Bayes, credit scoring and terrorism

Once again, my neighbor Corey did publish a very interesting post on his blog http://bayesianbiologist.com/… on how likely the NSA program will catch a terrorist (a real one). I was working on something similar last weeks, with Stéphane Tufféry, for our chapter, entitled Statistical Learning in Actuarial Science. The idea was to show credit scoring techniques, from logistic regression, classification trees, random forests, etc. Of course, it is more boring, since we talk about loans and not terrorism. In credit scoring, we consider possible loans, and we have to predict if someone is more likely to be a bad guy or a good guy. The idea is the same: based on some covariates, we need to build a score function, that can be related to the probability of being bad. The higher the score, the more likely the person will be a bad guy. Then, of course, we have to discuss errors, namely false positive (good guys that we think are bad) and false negative (bad guys that we think are good). From the company, you do not want to have bad guys in your portfolio, and from everyone else point of view (since everyone believes he is with the good one, this is a classical optimistic bias), we do not want to be confused with those bad guys. Then we can spend hours on classification curves, and criteria to assess if our classifier is good or not, etc. While I was writing the introduction of the chapter, I remember that I found it hard to find proper words (to describe that 0/1 problem). But I did use (like everyone else) the terms good and bad. Like in terrorism. Except that to use this terminology (bad and good), we have to be more specific. In credit scoring, a bad guy is someone who did not pay back, at least once, for instance. But in terrorism, I think it is more difficult to say what a terrorist is.

I mean, in France, we did experiment terrorism too, a few years ago. In December 1996, I was in a RER train, going South, and we reached Cité Universitaire when a bomb did explode in Port Royal. The train following mine I guess. I remember that a couple of days after, I was traveling Paris, in bus, carrying with me a nice plant of… a plant that you’re not supposed to grow. Say I was carrying sandwiches, according to Ted Mosby. So in order to avoid troubles (since I was not suppose to have this kind of plant species), I put it in a large box. I remember that people were starring at me, and some actually asked me what was in the box. So for some reasons, people try to build there own terrorist classifier, based on what they think might be covariates. And dirty trousers, not well shaved, long hair (yes, I used to have long hair) and box in the bus were obviously some of them. Note that I don’t blame them, I do the same! After reading Corey’s post this morning, I took the bus. And I saw someone with a ninja sword.

At first, my terrorist classifier put her (yes, I try to have a gender-free terrorist model) in the bad guy class. Then I understood it was an umbrella. So I put her in the super cool geeky category (that only a few can reach).

When I started to teach non-life insurance in Paris, the last part of the course was dedicated to large risks, natural catastrophes, and a hot topic: terrorism. I was giving this course (probably my best experience, ever) in tandem with François Bucchini, who was working by that time for AXA France. The two of us were giving the course together, interacting: I was the boring guy doing the maths, and François was sharing his experience. And by that time, he was involved in the creation of GAREAT, a market structure, launched in France in 2002, to propose reinsurance against terrorism (for French companies). And one of the first claim was from the CAV (which is a pun for Comité d’Action Viticole) considered as a terrorist group. So, as he told us, be careful of prejudices when you think about terrorism. Cool wine drinkers can be dangerous terrorists…

Actually, I would love to see covariates used by the NSA to predict if you’re a bad guy, or a potentially dangerous terrorist. Let us have a guess… You have asked for a visa for Pakistan? or Afghanistan? or Libya (not Libya, not yet bad guys still have good friends there)? You have a NRA membership? You bought some heavy metal on iTunes? You still have a stop acta sticker on your blog? you have a blog? you wrote a post including the word terrorist in it?

Note: I am supposed to be in Chicago next week. Si if I cannot enter in the U.S., we’ll probably know more about potential covariates.

Bayes is playing Russian roulette

There was (once again) a nice puzzle inhttp://www.futilitycloset.com/. Bayes and a good friend are playing Russian roulette. The revolver has six chambers. He puts two bullets in two adjacent chambers, spin the cylinder, hold the gun to his friend’s head, and pull the trigger. It clicks. So it is now Bayes’s turn: he can choose either to spin the cylinder again or leave it as it is. Which is better? Hopefully, Bayes knows his theorem: if he does spin it, the probability of getting killed is 2 out of 6 (four empty chambers out of six), but if he does not, since his friend is still alive, then the hammer should be next to one of the four cylinders in red, below


So here, there is 3 chance out of 4 to survive, i.e. the probability of getting killed is 1 out of 4 (while it was 1 out of 3 when spinning). So Bayes should not spin. And as always, it is possible to see it is a more general result: more generally, in a revolver with http://freakonometrics.blog.free.fr/public/perso5/bullet01.gif chambers, it there are http://freakonometrics.blog.free.fr/public/perso5/bullet02.gif bullets in http://freakonometrics.blog.free.fr/public/perso5/bullet02.gif adjacent chambers,  if the first player survives, the probability of getting killed is k over http://freakonometrics.blog.free.fr/public/perso5/bullet01.gif, when spinning, while it would be 1 over http://freakonometrics.blog.free.fr/public/perso5/bullet03.gif if we don’t. Not spinning is better if and only if

http://freakonometrics.blog.free.fr/public/perso5/bullet04.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso5/bullet05.gif

So you’d better not spin, unless there was one bullet in the revolver, i.e. http://freakonometrics.blog.free.fr/public/perso5/bullet06.gif… or http://freakonometrics.blog.free.fr/public/perso5/bullet07.gif (in that case, it might not be a good idea actually to play the game).

La formule de Bayes et les avis d’experts

Pour répondre à une question de Laurence (qui voulait que je réexplique la formule de Bayes), je voulais faire une digression par un exemple amusant, emprunté à Daniel Kahnmann et Amos Tversky. Mais tout d’abord revenons à la formule de Bayes, qui nous dit que la fonction

http://freakonometrics.blog.free.fr/public/perso4/bayyes-1.gif

est une mesure de probabilité,

http://freakonometrics.blog.free.fr/public/perso4/bayyes-2.gif

que l’on appellera probabilité conditionnelle. Avec cet outils, on peut essayer de résoudre le problème suivant

La question est alors

Daniel Kahnmann et Amos Tversky ont présenté ce problème à beaucoup d’étudiants, et comme ils le notent

Qu’en est-il ? Comme toujours, un peu de formalisme. La couleur C du taxi prend deux valeurs, Bleue (B) ou Vert (G). Quant au témoin (W), il peut être correct (C) ou parfois se tromper (F). Les informations dont on dispose sont

http://freakonometrics.blog.free.fr/public/perso4/bayeees04.gif et http://freakonometrics.blog.free.fr/public/perso4/bayeees05.gif

(pour la couleur des taxis) et pour le fait de se tromper

http://freakonometrics.blog.free.fr/public/perso4/bayeees02.gif et http://freakonometrics.blog.free.fr/public/perso4/bayeees03.gif

D’après la formule de Bayes

http://freakonometrics.blog.free.fr/public/perso4/bayeees07.gif

http://freakonometrics.blog.free.fr/public/perso4/bayeees08.gif

On peut réécrire la première relation

http://freakonometrics.blog.free.fr/public/perso4/bayeees09.gif

et en faisant pareil pour la seconde, on en déduit

http://freakonometrics.blog.free.fr/public/perso4/bayeees06.gif

donc finalement le rapport des probabilités d’avoir un taxi bleu, ou vert, s’écrit

http://freakonometrics.blog.free.fr/public/perso4/bayeees10.gif

Comme on sait que http://freakonometrics.blog.free.fr/public/perso4/bayeees11.gif on en déduit

http://freakonometrics.blog.free.fr/public/perso4/bayees12.gif

On notera que malgré notre avis d’expert, la probabilité que la voiture soit bleue est toujours plus grande que la probabilité que la voiture soit verte. Même avec une faible probabilité de se tromper (15%) on ne lui accorde guère de crédit… Quel est le crédit que l’on doit accorder à un expert ? Les experts n’ont-il pas tendance à sur-estimer leurs connaissances… C’est un peu ce que je lisais dans l’article d’Eric Van den Steen, Overconfidence by Bayesian-Rational Agents [à suivre]

Il y a des jours avec, et des jours sans

Je ne sais pas si ça vous est déjà arrivé, mais quand on fait des crêpes, la probabilité d’en réussir une est proportionnelle à la probabilité que l’on vient de réussir les dernières. Autrement dit si j’ai réussi les premières, je suis dans un bon jour, et je réussi davantage celles qui suivent. Alors que si je commence par rater les premières, je suis dans un mauvais jour, et je les rate toutes (pour les sportifs, il est aussi possible de noter que, par exemple, au tennis, la probabilité de réussir son service est d’autant plus forte que j’ai réussi les derniers: c’est pareil).

  • Une première formalisation du problème

Supposons que sur les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais (qui servent d’étalonnage), http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif fois j’ai réussi ma crêpe, et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif fois, je l’ai ratée (trop de trous, trop irrégulière, trop cuite, que sais-je encore). Mon a priori sur la probabilité de ne pas avoir de chance (et de rater ma crêpe) est http://freakonometrics.blog.free.fr/public/perso3/crepes-04.gif. Supposons que la probabilité de réussir à l’instant http://freakonometrics.blog.free.fr/public/perso3/crepes-06.gif (on commence à compter après les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais) soit la fréquence empirique observée sur le passé. Autrement dit

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est le nombre de succès. En particulier http://freakonometrics.blog.free.fr/public/perso3/crepes-09.gif. On retrouve ici le fait que plus j’en ai raté, plus j’ai de chances de rater les prochaines… c’est la mauvaise main en quelque sorte.
Afin de mieux comprendre, supposons http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif. Autrement dit, sur mes deux premiers essais (qui ne servent qu’à étalonner), j’ai raté et réussi une fois. Juste avant l’instant 1, j’ai deux possibilité: réussir ou rater, avec la probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-11.gif. I.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-12.gif

Si on se place une crêpe plus loin, et que l’on regarde les probabilités conditionnelles

http://freakonometrics.blog.free.fr/public/perso3/crepes-13b.gif

avec quelque-chose de similaire si on conditionne par l’autre événement,

http://freakonometrics.blog.free.fr/public/perso3/crepes-14.gif

Et donc la probabilité non conditionnelle s’écrit

http://freakonometrics.blog.free.fr/public/perso3/crepes-15.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-16.gif
http://freakonometrics.blog.free.fr/public/perso3/crepes-18.gif

et

http://freakonometrics.blog.free.fr/public/perso3/crepes-19.gif

Moralité

http://freakonometrics.blog.free.fr/public/perso3/crepes-20.gif

Autrement dit, on poursuivant la récurrence, on a que http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est uniformément distribué sur http://freakonometrics.blog.free.fr/public/perso3/crepes-21.gif
Étonnant, non ? Autrement dit, si j’ai commencé par un échec et une réussite, la probabilité d’avoir 50 réussites sur 100 essais sera très exactement http://freakonometrics.blog.free.fr/public/perso3/crepes-23.gif. Qui est la même probabilité que d’avoir 1 réussite, ou 99 !
On a aussi que

http://freakonometrics.blog.free.fr/public/perso3/crepes-22.gif

Mais que se passe-t-il si http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif ne sont plus égaux à 1 ?

  • Lien avec l’approche bayésienne et la théorie de la crédibilité

Ce qui est amusant, c’est que la probabilité

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

est bien connue par les personnes qui connaissent l’assurance et la théorie de la crédibilité. Supposons que je note http://freakonometrics.blog.free.fr/public/perso3/crepes-24.gif la variable indicatrice correspondant à une réussite ou un succès au http://freakonometrics.blog.free.fr/public/perso3/crepes-26.gifème essai (avec probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-25.gif). D’un point de vu bayésien, la loi conjuguée de la loi binomiale (on va compter les réussites sur n essais) est une loi Beta. Rappelons que la loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif a pour densité

http://freakonometrics.blog.free.fr/public/perso3/crepes-28.gif

sur http://freakonometrics.blog.free.fr/public/perso3/crepes-29.gif, et que l’espérance de cette variable est http://freakonometrics.blog.free.fr/public/perso3/crepes-30.gif, i.e. http://freakonometrics.blog.free.fr/public/perso3/crepes-31.gif. Donc cette modélisation devrait parfaitement convenir à notre problème. Ce qui est intéressant, c’est que comme on travaille avec la loi conjuguée, la loi ex post est encore une loi Beta. Autrement dit, ici, si je suppose que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit a priori une loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif, et que la loi de http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif sachant http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit une loi binomiale http://freakonometrics.blog.free.fr/public/perso3/crepes-32.gif, alors la loi a posteriori de http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif sachant http://freakonometrics.blog.free.fr/public/perso3/crepes-33.gif est une loi http://freakonometrics.blog.free.fr/public/perso3/crepes-34.gif (c’est ce que l’on avait vu ici). Et en particulier

http://freakonometrics.blog.free.fr/public/perso3/crepes-35.gif

Autrement dit on obtient une combinaison convexe entre la moyenne de la loi a priori (à gauche) et la fréquence empirique (à droite).
On retrouve ici ce que nous avions noté auparavant, en notant que http://freakonometrics.blog.free.fr/public/perso3/crepes-36.gif, et http://freakonometrics.blog.free.fr/public/perso3/crepes-37.gif.
On peut toutefois se demander si le résultat que l’on avait obtenu sur la loi de S_n reste valide ici…
Notons que

http://freakonometrics.blog.free.fr/public/perso3/crepes-38.gif

i.e. (moyennant quelques calculs)

http://freakonometrics.blog.free.fr/public/perso3/crepes-40.gif

Dans le cas où http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif, en réécrivant aussi le premier terme, on a

http://freakonometrics.blog.free.fr/public/perso3/crepes-41.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-42.gif

On retrouve bien la formule obtenue auparavant… Ce qui est intéressant aussi, c’est que si on note que

http://freakonometrics.blog.free.fr/public/perso3/crepes-43.gif

or

http://freakonometrics.blog.free.fr/public/perso3/crepes-44.gif

et donc

http://freakonometrics.blog.free.fr/public/perso3/crepes-45.gif

et par récurrence, on obtient que

http://freakonometrics.blog.free.fr/public/perso3/crepes-46.gif

et que de plus

http://freakonometrics.blog.free.fr/public/perso3/crepes-47.gif

Autrement dit, effectivement, si j’en ai raté beaucoup au début, ma probabilité d’en rater beaucoup est grande (et elle est d’ailleurs constante). Certains parleraient de no-free lunch.

Will I ever be a bayesian statistician ? (part 2)

A few weeks ago, I started a series of posts on the magic of bayesian statistics from the eyes of a muggle (see http://freakonometrics.hypotheses.org/2191). It might be time to go a bit further…. And today, I wanted to discuss the choice of the a priori distribution of the parameter (which was mentioned in the commentary of the previous post). As far as I understood, there are several houses with different ideas on how to choose it.

  • The conjugate house

The first idea (used in the previous post, here) is to consider an exponential distribution. To be formal, those distributions can be written as

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-10.gif

(in a form as general as possible), i.e.

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-11.gif

Here http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif is somehow the new parameter of the distribution. Then, a conjugate priorhttp://freakonometrics.blog.free.fr/public/perso2/bayes-prior-12.gif for the parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif of the exponential family is given by

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-15.gif

where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-16.gif (where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-17.gif is the dimension of http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif).

The conjugate prior is interesting since, when combined with the likelihood (and normalized), produces a posterior distribution which is of the same type as the prior. And a lot of standard distributions have a conjugate prior. E.g.

  • For a Bernoulli distribution, i.e. http://freakonometrics.blog.free.fr/public/perso2/conj-00c.gif are i.i. with distribution http://freakonometrics.blog.free.fr/public/perso2/conju-01.gif, assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, , then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conju-03.gif
  • For a binomial distribution,  http://freakonometrics.blog.free.fr/public/perso2/conju-05.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conj-06.gif
  • For a Negative Binomial distribution, http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-20.gif,  assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posterioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conj-08.gif
  • For a Poisson distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-09.gif, assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conj-10.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gifis still gamma, with parameter
http://freakonometrics.blog.free.fr/public/perso2/conj-13.gif
  • For a Geometric distribution,http://freakonometrics.blog.free.fr/public/perso2/conj-15.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-22.gif
  • For an Exponential distributionhttp://freakonometrics.blog.free.fr/public/perso2/conj-23.gif assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif ishttp://freakonometrics.blog.free.fr/public/perso2/bayes-prior-21.gif , then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-22.gif
  • For a Gaussian distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-28.gif assume that http://freakonometrics.blog.free.fr/public/perso2/conj-31.gif, then
http://freakonometrics.blog.free.fr/public/perso2/conj-35.gif

  • For a gamma distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-43.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is  http://freakonometrics.blog.free.fr/public/perso2/bayes-000000.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-23.gif
  • For a Pareto distribution, , assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-21.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters

http://freakonometrics.blog.free.fr/public/perso2/bayes----ooooo.gif

  • The non-informative or vague house

So far, the choice of the prior was not neutral, in the sense that the a priori of the statistician will have an influence on a posteriori distributions (we’ll discuss that point further later on). We could be interested by the case where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-12.gif is somehow t neutral. A famous example is the case of  http://freakonometrics.blog.free.fr/public/perso2/conju-01.gifdistribution. Between 1745 and 1784, Pierre Simon Laplace observed 393,386 birth of boys versus 377,555 birth of girls (or 251,527 boys versus 241,945 girls if we consider the initial article, for birth before 1770). He wanted to quantify the probability that p, the provability to have a boy, exceed 1/2. He assume that a priorip was uniform on the unit interval claiming that it was being as neutral as possible. But it is not that correct.
The idea of noninformative prior is that we should get an equivalent result when considering a transformed parameter. So assume that the parameter is no longer http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-03.gif, but http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif, where http://freakonometrics.blog.free.fr/public/perso2/bayes-----0000.gif (for some bijective transformation). The distribution (density) of http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif is then

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-01.gif

Let http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-02.gif denote Fisher information of parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-03.gif,

http://freakonometrics.blog.free.fr/public/perso2/bayes-priori-04.gif

Then, Fisher information of parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif is

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-06.gif

which can be written

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-07.gif

So if we want a distribution invariant by transformation of the parameter (http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-08.gif), it seems natural to consider

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-09.gif

or at least something proportional to that square root, since we want to get a density.
Thus, from Jeffrey’s principle, the prior distribution for a single parameter is noninformative if it is taken proportional to the square root of Fisher’s information measure. For those who want to go further, see Noninformative Priors Do Not Exist or  A Catalog of Noninformative Priors.

  • For the Poisson distribution, the Jeffreys prior for the rate parameter  is
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-poisson.gif
  • For the Bernoulli distribution, the Jeffreys prior for the probability parameter  is
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-binomial.gif

This is the arcsine distribution and is a beta distribution with parameters 1/2.

  • The expert house

The idea is quite simple. We need a prior distribution so that

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-32.gif

But assume that we have already seen a similar problem before. For instance, I remember Eric Parent mentioning the case of river flow models. If we have two similar rivers, then it might be interesting to use information on one river as an a prioriinformation for the second one, something like

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-33.gif

I guess it is also possible to use meta-regression to get an aggregation of experts opinion.

To go further on bayesian statistics, I suggest to go on Albus Dumbledore’s og,here, or the the blog of some PhD (and postdoc) students in Hogwarts, there. Or if you can wait, a dozen other posts will come soon (well, let’s hope so). The next one will probably be on a posteriori calculations (which is the natural step since we’ve seen a priori choice here).

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 http://freakonometrics.free.fr/blog/bayy1.png i.i.d. with distribution http://freakonometrics.free.fr/blog/bayy2.png. Here we note

http://freakonometrics.free.fr/blog/bayy3.png

since parameter http://freakonometrics.free.fr/blog/bayyyyy001.png is a random variable. The idea is to assume that http://freakonometrics.free.fr/blog/bayyyyy001.png has a (so called a priori) distribution, e.g.

http://freakonometrics.free.fr/blog/bayy4.png

So far it was simple. The idea is then to consider the posterior distribution of http://freakonometrics.free.fr/blog/bayyyyy001.png, given the observations http://freakonometrics.free.fr/blog/bayyyyyy02.png. Thus, we need to compute the distribution of http://freakonometrics.free.fr/blog/bayyyyyy03.png which is here extremely simple (due to properties of the Gaussian distribution), i.e.

http://freakonometrics.free.fr/blog/bayyyyyy04.png

where

http://freakonometrics.free.fr/blog/bayyyyyy05.png

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

http://freakonometrics.free.fr/blog/bayy7.png

Thus,

http://freakonometrics.free.fr/blog/bayy8.png

and

http://freakonometrics.free.fr/blog/bayy9.png

From Bayes formula,

http://freakonometrics.free.fr/blog/bayy10.png

and we get easily

http://freakonometrics.free.fr/blog/bayy11.png

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

http://freakonometrics.free.fr/blog/bayy12.png
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,

http://freakonometrics.blog.free.fr/public/perso/bayyyyyyyyyyy.png

Note that

http://freakonometrics.free.fr/blog/bayy13.png

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.

http://freakonometrics.free.fr/blog/bayes-cv-1.gif

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

http://freakonometrics.free.fr/blog/bayes-cv-3.gif

or on the right
http://freakonometrics.free.fr/blog/bayes-cv-2.gif

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

http://freakonometrics.free.fr/blog/bayes-cv-gauss.gif

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

http://freakonometrics.free.fr/blog/bayes-cv-all.gif

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

Crédibilité et calcul bayésien

On a commencé en cours la tarification a posteriori, en commençant par la tarification bayésienne, et en présentant ensuite l’approche d’Hans Bühlmann comme un cas particulier en projetant dans un sous espace de combinaisons linéaires (un peu comme l’économétrie linéaire était présenté en cours comme un cas particulier de calcul d’espérance conditionnelle).La classification a priori sera poursuivie un TD, sur des bases de données, afin d’orienter davantage le cours sur la pratique des GLM.
L’utilisation des techniques bayésiennes en actuariat ne s’est pas faite sans peine. Pour ceux qui souhaitent s’en convaincre, il suffit de relire les comptes rendus de la conférence à Trieste en 1962. Dans le discours de Bruno de Finetti (publié dans l’ASTIN Bulletin), on retrouve un point historique, présentant les nouveaux travaux d’Hans Bühlmann,

Et Bruno de Finetti de conclure en rappelant les apports d’Arthur Bailey en 1950,

  • La formule de Bayes, ou apprendre par l’expérience

La formule de Bayes ne dit pas grand chose quand on le regarde sous la forme suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-01.png

A la rigueur on peut aussi la réécrire en utilisant deux fois cette formule

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-02.png

Bon, là encore, on ne voit peut être toujours pas ce qu’on peut en faire…. L’idée géniale est de faire la lecture suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-03.gif

Autrement dit la loi a posteriori du paramètre, conditionnellement aux observations (en bleu) est proportionnelle au produit de  la loi a priori du paramètre (en vert) et de la vraisemblance des observations (en rouge, qui sera un simple produit de lois si on suppose l’indépendance conditionnelle au paramètre). Le paramètre en noir est juste un facteur de normalisation, histoire que la loi soit effectivement une mesure de probabilité. D’où parfois l’écriture

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-bayes-4.gif

Bref, on a réussit à inverser les deux termes dans la probability (on retrouve le terme “inverse probability” dans la littérature anglo-saxonne).
En cours, on avait vu que si on supposait le nombre annuel de sinistres Poisson, conditionnellement au paramètre, avec une indépendance (toujours conditionnellement à ce paramètre), et que l’on suppose que la loi a priori du paramètre est une loi Gamma, alors la loi a posteriori se calcule facilement et c’est une loi Gamma.
Considérons la petite simulation suivante: on suppose que le nombre de sinistres suit une loi de Poisson, et on cherche à estimer le paramètre de cette loi. On commence par une loi a priori Gamma, centrée sur 5, et on regarde évoluer la loi a posteriori au fur et à mesure que l’on observe les sinistres arriver (ou les années passer en l’occurence).

On peut également considerer l’évolution de l’espérance a posteriori du paramètre, ainsi qu’un intervalle de confiance (puisqu’on connaît la loi)

La “magie” des méthodes bayésiennes est que peu importe la loi a priori, on va converger vers la vraie valeur. Par exemple sur l’animation ci-dessous, on part d’une loi a priori de moyenne 2,

avec là aussi l’évolution de la moyenne a posteriori avec un intervalle de confiance,

La morale de l’histoire est que comme la loi a priori n’influence pas (j’exagère un peu car elle va influencer la vitesse de la convergence, ce qui n’est pas rien en pratique), autant prendre quelque chose de simple pour les calculs: la loi conjuguée de la loi exponentielle. La loi conjuguée d’une loi de la famille exponentielle de la forme

s’écrit

où et sont des constantes, et où  est une constante de normalisation (afin que cette fonction soit effectivement une mesure de probabilité).

  • L’hérésie de la statistique bayésienne

La principale hérésie des statisticiens bayésien a été de proposer une lecture des probabilités d’un point de vue subjectif. En particulier Jimmie Savage et Bruno de Finetti partait de l’hypothèse que chacun avait à sa disposition son propre espace de probabilité.C’est ce qui lui permet de se construire un a priori. Et au fur et à mesure que le temps passe, on apprend, et on peut dériver une loi a posteriori. Par exemple  je pars d’un a priori comme quoi un enfant sur deux qui nait est un garçon, en France. Je peux ensuite récupérer des données, et affiner ma connaissance. C’est ce qu’avait fait Laplace il y a fort fort longtemps….
L’hérésie est ici un
problème philosophique lié au concept même de probabilité. Mais on peut noter qu’à l’époque dans les années 60, les débats étaient assez virulants. Et pour Jimmie Savage et Bruno de Finetti, il n’est pas nécessaire d’avoir le formalisme mathématique d’Andreï Kolmogorov pour appréhender la notion de probabilité (c’est à dire la mesure d’ensembles dans un espace abstrait d’évènements, voir ici ou ). On notera que pour John Keynes (qui avait écrit un “traité de probabilités“) il existait une probabilité subjective, mais davantage liée à une information donnée par l’éducation. C’est un problème culturel plus que personnel. Mais finalement, il a surtout fallu attendre les travaux de Daniel Kahneman et Amos Tversky dans les années 70 pour que des psychologues se penchent sur la construction individuelle des probabilités.

  • La crédibilité au sens de Bühlmann

L’idée d’Hans Bühlmann est de proposer une formule fermée qui approche l’espérance conditionnelle, tout simplement en utilisant celle qui avait été obtenue dans le cas bayésien. En tarification, on s’intéresse à des primes pures, et il est finalement assez naturel de chercher à calculer

Hans Bühlmann a suggéré d’approcher cette grandeur par

avec classiquement

est appelé facteur de crédibilité, avec , et où respectivement,


Bref, on a un facteur de crédibilité qui dépend de la variance intra et de la variance inter, pour reprendre des termes connus en statistique. L’énorme avantage est que cette formule peut très facilement s’utiliser sur des données. En particulier, considérons un tableau de la forme suivante

où on regarde si un contrat a – ou n’a pas – eu de sinistre chacune des années. Formellement, on dispose d’un tableau d’observations,

(on notera que le nombre d’années d’observations pour chacun des contrats peut être différent). On peut alors obtenir les estimateurs empiriques des moments nécéssaires, i.e. on pose

correspondant à la moyenne empirique pour un contrat,

et

Le facteur de crédibilité pour le contrat i est alors

où tout naturellement

Intervalle de confiance pour une proportion

Je voulais poster un petit billet puisque c’est la deuxième fois que l’on me pose une question similaire. Pour déterminer un intervalle de confiance pour une proportion, on connaît tous la formule

La question qui m’avait été posée est simple: que se passe-t-il si n est petit, et si – empiriquement – on a trouvé une proportion de 100% ? Dit autrement: “on fait 5 lancers de “pile” ou “face” et que on obtient 5 fois “pile”: que peut-on dire sur p, la probabilité de tomber sur “pile” ?”

On utilisant la formule précédente, on est un peu bloqué…

  • La première réponse pourrait être fréquentiste, et consiste à proposer une lecture duale de l’intervalle de confiance, et des abaques.

Le graphique de droite montre, en trait violet, l’intervalle de confiance approché, c’est à dire avec l’intervalle de confiance indiqué ci-dessus. La courbe bleue, encadrant la zone verte correspond à l’intervalle de confiance calculé à partir des quantiles d’une loi binomiale de paramètre la fréquence empirique observée. Ici, on aurait envie de dire que p a 95% de chance d’être supérieure à 47.8%.

En fait, cette “lecture duale” des abaques revient à construire – me semble-t-il – l’intervalle de confiance de Clopper-Pearson. L’intervalle de confiance a ici la forme suivante

Notons que c’est cette technique qui est utilisée sous R, via la fonction binom.test du package de base. Ici, on obtient

> binom.test(x, n, p = 0.5,conf.level = 0.95)
Exact binomial test
data: x and n
number of successes = 5, number of trials = 5, p-value = 0.0625
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.4781762 1.0000000
sample estimates:
probability of success
1

…. ce qui correspond exactement à l’intervalle de confiance calculé précédemment. On retrouve d’ailleurs cette fonction dans binconf i.e.

> binconf(x = 5, n = 5, alpha = .05, method = "exact")
PointEst Lower Upper
1 0.4781762 1

qui donne (c’est rassurant) toujours la même chose. On peut aussi programme l’intervalle de confiance proposé par Heldge Blacker, qui donne ici un intervalle de confiance de la forme [50%; 100%]. Pour aller plus loin, je pourrais renvoyer au papier de Lawrence Brown, Tony Cai et Anirban DasGupta, dans Statistical Science sur ce sujet.

  • Une seconde réponse pourrait être une approche bayésienne, qui d’ailleurs correspond exactement à l’expérience de Bayes (1763) ou Laplace (1786).

On suppose que p est une variable aléatoire. On se donne une information a priori, et on peut se demander ce que devient la loi a posteriori de p, sachant que 5 “piles” ont été observés. On suppose que les tirages sont indépendants, et donc la probabilité d’avoir k fois “piles” sur n tirages suit une loi binomiale, conditionnellement à p, i.e.

Le plus simple est de considérer comme loi a priori la loi conjuguée de la loi binomiale, à savoir une loi beta. Dans ce cas, si p suit a priori une loi B(a,b), alors la loi a posteriori de p, sachant que x “pile” ont été observés est une loi B(a+x,b+n-x).

expérience décrite

La courbe en rouge est la densité de la loi a priori (uniforme ou beta de paramètres 2 et 2), et la courbe bleue, la loi a posteriori, avec un intervalle de confiance à 95%. Sur l’exemple de droite, par exemple, on a 95% de chance que p soit compris entre 47,3% et 96,8%.

Article (de vulgarisation) sur la crédibilité

Parution d’un court article sur la Théorie de la Crédibilité dans la revue Risques.

La théorie de la crédibilité est une jolie application de l’approche bayésienne de la statistique, puisque l’on cherche à déterminer la probabilité d’appartenir à une classe de risque conditionnellement à une expérience de sinistralité observée. La théorie développée par Hans Bühlmann propose d’ailleurs de jolies applications, en particulier sur les jeux de fléchettes, où l’on compare la moyenne sur un sous-groupe homogène (ou supposé l’être) avec la moyenne de la population totale. En faisant intervenir des variance intra et des variances inter, on arrive à un résultat relativement simple et élégant.