Introduction au calcul Bayésien
- Théorème de Bayes
Le modèle Bayésien et théorie de la décision
- Notions de la théorie de décision.
  - Exemples : $\Theta=\mathcal{E}=\mathbb{R}$
- Le cas bayésien.
  - Risque Bayésien
  - Estimateur Bayésien
Le choix de la loi a priori
- Lois a priori conjuguées.
- Lois a priori de Jeffrey.
Modèle Bernoulli

Introduction au calcul Bayésien

(Séance 1, 3/9/2016)

Théorème de Bayes

Soit $(\Omega, \mathcal{A},\mathbb{P})$ un espace de probabilité. Soient $X$ et $Y$ deux évènements de $\mathcal{A}$. Alors \[ \mathbb{P}\left( X\mid Y\right)= \displaystyle\frac{ \mathbb{P}(X\cap Y)}{ \mathbb{P}(Y)}\]

Exercice

On considère un test pour une infection. On suppose que la sensitivité du test est égale à 99$\%$ (la probabilité qu’un infecté soit détécté positif) et que la spécificité du test est égale à 97$\%$ (la probabilité qu’une personne saine soit détéctée négative). Sachant que la proportion des personnes infectées dans une population est égale à 1$\%$.

Question 1 Quelle est la propbabilité qu’une personne ayant eu un résultat positif au test soit détectée infectée.
Question 2 Cette personne passe ce test une deuxième fois et il est encore positif. Quelle est la la propbabilité qu’elle soit infectée.

Solution avec `R`

Question 1

On simule une population avec un taux d’infection égal à 1$\%$.

> N=10000 # Taille de la population
> p=0.01 # Probabilité d'infection
> population=rbinom(N,1,prob = p)
> N_infected=sum((population==1));N_infected ## Nombre d'infectés
[1] 101
> infected=which(population==1)

On simule le test

> sensitivity=.99
> specificity=.97
> test=rbinom(N,1,(population==1)*sensitivity+(population==0)*(1-specificity))
> N_test=sum((test==1));N_test ## Nombre de positifs.
[1] 393

La matrice de confusion test x population

> confusion=xtabs(~population+test);confusion
          test
population    0    1
         0 9607  292
         1    0  101

Donc la probabilité qu’une personne ayant eu un résultat positif au test soit détectée infectée se calcule de la façon suivante

> confusion[2,2]/(confusion[1,2]+confusion[2,2])
[1] 0.2569975

Question 2

Simulant une deuxième fois le test

> test2=rbinom(N,1,(population==1)*sensitivity+(population==0)*(1-specificity))

Considérons le croisement entre les trois variables population, test et test2

> confusion2=xtabs(~population+test2+test);confusion2
, , test = 0

          test2
population    0    1
         0 9298  309
         1    0    0

, , test = 1

          test2
population    0    1
         0  280   12
         1    0  101

Donc la probabilité qu’une personne infectée sachant qu’elle a réagit positivement deux fois au test est

> confusion2[2,2,2]/(confusion2[2,2,2]+confusion2[2,1,2])
[1] 1

Solution théorique

On note par $M$ l’événement qu’une personne soit infecetée et par $P$ qu’elle soit détectée positive. Alors \[\mbox{sensitivité}=\mathbb{P}(P\mid M)=\displaystyle\frac{\mathbb{P}(P\cap M)}{\mathbb{P}(M)}=.99\] et \[\mbox{spécificité}=\mathbb{P}(\overline{P}\mid \overline{M})=\displaystyle\frac{\mathbb{P}(P\cap M)}{\mathbb{P}(\overline{M})}=.97\]

Donc la la probabilité qu’une personne ayant eu un résultat positif au test soit détectée infectée se calcule

\[ \mathbb{P}(M\mid P)=\displaystyle\frac{\mathbb{P}(P\cap M)}{\mathbb{P}(P)}=\displaystyle\frac{\mathbb{P}(P\mid M)\mathbb{P}(M)}{\mathbb{P}(P\mid \overline{M})\mathbb{P}(\overline{M})+\mathbb{P}(P\mid M)\mathbb{P}(M)}=\displaystyle\frac{.99 \times .01}{(1-.97)\times (1-.01)+ .99 \times .01}=0.25\]

Le modèle Bayésien et théorie de la décision

(Séance 2, 10/9/2015)

Notions de la théorie de décision.

Soit $\mathcal{M}$ un modèle paramétrique : Soit $\underline{X}=(X_1,\ldots,X_n)$ un $n-$échantillon de loi $\mathbb{P}_\theta$ dépendant d’un paramètre $\theta\in \Theta \subseteq \mathbb{R}^p$.

Rappelons un estimateur de $\theta$ est une fonction du $n-$échantillon $\underline{X}$ : \[\widehat{\theta}(\underline{X})\in \mathcal{E}\supseteq \Theta\]
Comment évaluer la qualité d’un estimateur $\widehat{\theta}$ ?
Fonction de Perte : \[L\,:\,\Theta \times \mathcal{E} \longrightarrow \mathbb{R}^+\] telle que $L(\theta,\theta)=0$ pour tout $\theta\in \Theta$. Donc l’erreur de l’estimation de $\widehat{\theta}$ est mesurée par $L(\theta,\widehat{\theta})$.
Fonction de Risque \[\theta\longmapsto R(\theta,\widehat{\theta})=\mathbb{E}\left[L(\theta,\widehat{\theta})\right]\]

Exemples : $\Theta=\mathcal{E}=\mathbb{R}$

Perte quadratique

$L(\theta,\widehat{\theta})=(\theta-\widehat{\theta})^2$ et le risque associée est le risque quadratique \[ R(\theta,\widehat{\theta})=\mathbb{E}\left[(\theta-\widehat{\theta})^2\right]\]

Perte absolue

$L(\theta,\widehat{\theta})=|\theta-\widehat{\theta}|$ et le risque associée est le risque quadratique \[ R(\theta,\widehat{\theta})=\mathbb{E}\left[|\theta-\widehat{\theta}|\right]\]

Intuitivement plus le risque de l’estimateur est faible plus l’estimateur $\widehat{\theta}$ est considéré comme performant.

Exercice

Soit $\underline{X}=(X_1,\ldots,X_n)$ telle que pour tout $i=1,\ldots,n$, $X_i\sim \mathcal{N}(\theta,1)$. On considère $\overline{X}$ comme estimateur de $\theta$. Calculer le risque quadratique de $\overline{X}$.

Le cas bayésien.

Risque Bayésien

Dans le cas bayésien on suppose qu’on dispose d’une certaine information a priori sur le paramètre $\theta$ qu’on la traduit avec une distribution de probabilité $\theta$ notée $\pi(\theta)$.

Notre objectif est de mettre à jour cette information une fois qu’on a observé l’échantillon des données.

Si on suppose que $\underline{X}=(X_1,\ldots,X_n)$ est issu d’une loi de probabilité de densité $f(x\mid \theta)$. Alors la fonction vraisemblance \[L(\theta,\underline{x})=\prod_{i=1}^n f(x_i\mid \theta)\] On calcule alors ce qu’on appelle la loi a posteriori de $\theta$ : \[\pi(\theta\mid \underline{x}) \propto L(\theta,\underline{x}) \pi(\theta)\]

Exemples

Calculer la loi a posteriori dans les cas suivants

$f$ est la loi Bernoulli$(\theta)$ et $\pi\sim \mbox{Beta}(a,b)$.
$f$ est la loi $\mathcal{N}(\theta,1)$ et $\pi\sim \mathcal{N}(0,1)$.

Risque Bayésien On considère une fonction $L(\dot,\dot)$ et modèle statistique paramètrique $\mathcal{M}=\left\{\mathcal{X},\mathcal{A},\mathbb{P}_\theta, \theta\in \Theta \right\}$. On considère une loi a priori $\pi$ sur le paramètre $\theta$. Le risque bayésien d’une estimateur $\widehat{\theta}(\underline{x})$ \[ \mathcal{R}_\pi\left(\widehat{\theta}(\underline{x})\right)=\displaystyle\int_\Theta R\left( \theta, \widehat{\theta}(\underline{x})\right) \pi(d\theta)\]

Estimateur Bayésien

L’estimateur $\widehat{\theta}^B(\underline{X})$ est appelé estimateur de Bayes associé à la fonction de perte $L$ et à l’a priori $\pi$ c’est l’estimateur qui minimise le risque bayésien \[ \mathcal{R}_\pi\left(\widehat{\theta}^B(\underline{X})\right)= \min_{\widehat{\theta}(\underline{x})\in \mathcal{E}} \mathcal{R}_\pi\left(\widehat{\theta}(\underline{x})\right)\]
Remarque : L’expression de l’estimateur bayésien dépend du choix de la loi a priori $\pi$ et de la fonction perte $L$.
Théorème : Supposons que
1. Il existe un estimateur $\widehat{\theta}_0(\underline{X})$ tel que $\mathcal{R}_\pi\left(\widehat{\theta}_0(\underline{X})\right) <\infty$.
2. Pour presque tout $x \in \mathcal{X}$ il existe \[ \widehat{\theta}(\underline{x}) =\mbox{argmin}_{y\in \mathcal{E}} \mathbb{E}\left[ L(\theta,y)\mid X=x \right]\]
3. La fonction $x\longrightarrow \widehat{\theta}(\underline{x})$ est mesurable.

Alors \[\widehat{\theta}(\underline{X})=\widehat{\theta}^B(\underline{X})\]

Corollaire Considérons le cas $\Theta=\mathcal{E}=\mathbb{R}$
1. Si $L(\theta,\theta')=(\theta-\theta')^2$ est la fonction perte quadratique, alors l’estimateur bayésien est donné par la moyenne a posteriori \[ \widehat{\theta}^B(\underline{X})=\mathbb{E}(\theta\mid X)\] De plus $\mathcal{R}_\pi\left(\widehat{\theta}^B(\underline{X})\right)=\mathbb{E}\left[\mathbb{V}(\theta\mid X)\right]$
2. Si $L(\theta,\theta')=|\theta-\theta'|$, alors l’estimateur bayésien est donné par la médiane a posteriori \[ \widehat{\theta}^B(\underline{X})=\mbox{Médiane}_{\theta|X}\]
Théorème Si un estmateur bayésien est sans biais alors il a un risque bayésien nul.

Exercice

Soit $X=(X_1,\ldots,X_n)$ un échantillon i.i.d de loi normale de moyenne $\theta$ inconnue et de variance $\sigma^2$. Considérons la loi a priori $\pi=\mathcal{N}(0,\sigma_0^2)$ 1. Calculer la loi a posteriori sur $\theta$. 2. En déduire l’expression de l’estimateur bayésien associé à la fonction perte quadratique. 3. Calculer le risque bayésien de l’estimateur bayésien. 4. Répondre aux questions 1,2,3 avec la fonction perte $L(\theta,\theta')=(\theta-\theta')^2 \exp(\theta^2)$. On considérera le cas où $\sigma^2=1$.

Le choix de la loi a priori

(Séance 3, 17/9/2016)

Pour construire une loi a priori on doit distinguer 3 cas possibles :

On dispose d’une information a priori partielle, donnée par l’expérience.
On dispose d’une quantité limité : lois conjuguées.
On dispose d’une quantité limité ou même aucune information sur $\theta$ : la loi a priori de Jeffrey.

Lois a priori conjuguées.

Une loi a priori conjuguée sur $\theta$ pour un modèle statistique donné si la loi a posteriori sur $\theta$ appartient à la même famille de la loi a priori.

Exercice.

Calculer dans les cas suivants la loi a priori conjuguée sur $\theta$ et calculer la loi a posteriori :

$\mathcal{N}(\theta,\sigma^2)$ où $\sigma^2$ est connue.
$\mbox{Poisson}(\theta)$
$\Gamma(\nu,\theta)$ où $\nu$ est connue.
$\mbox{Binomiale}(N,\theta)$
$\mathcal{N}(\mu,\theta^{-1})$ où $\mu$ est connue.

Lois a priori de Jeffrey.

On suppose que le $n-$échantillon $X_1,\ldots,X_n$ de densité $f(\,\bullet\mid \theta)$. La loi de Jefferey se construit à partir seulement la connaissance de la densité $f(\,\bullet\mid \theta)$. Cette famille de lois est basée sur le calcul de l’information de Fisher : \[ I(\theta)=\mathbb{E}_\theta\left[ \left( \displaystyle\frac{\partial \mbox{ log }f(X\mid \theta)}{\partial \theta} \right)^2\right]\] Dans le cas de modèles réguliers. \[ I(\theta)=\mathbb{E}_\theta\left[ \displaystyle\frac{\partial^2 \mbox{ log }f(X\mid \theta)}{\partial \theta^2} \right]=\mathbb{E}_\theta\left[ \displaystyle\frac{\partial^2 \mbox{ log }l_n( \theta)}{\partial \theta^2} \right]\] où $l_n$ est le log-vraisemblance. Ces lois sont souvent des lois a priori non-informatives.

Définition la loi a priori de Jeffrey qu’on note $^ c’est la loi qui admet la densité suivante

Si $\Theta\subseteq \mathbb{R}$ \[ q(\theta)\propto I^{-1/2}(\theta)\]
Si $\Theta\subseteq \mathbb{R}^d$ \[ q(\theta)\propto |I(\theta)|^{-1/2}\]

Exercice

Reprendre les modèles statistiques de l’exercice précedent et calculer les lois a priori de Jeffrey correspondantes.

Exercice

On considère $\underline{x}=\left(x_{1},\ldots,x_{n}\right)$ un échantillon i.i.d d’une variable $X$ de loi de densité \[ f\left(x\mid\theta,\lambda\right)\propto\exp\left(-\theta x^{2}+\lambda x\right) \]

définie sur tout $\mathbb{R}$.

Montrer que la loi de $X$ est une loi $\mathcal{N}\left({\displaystyle \frac{\lambda}{2\theta}},\,{\displaystyle \frac{1}{2\theta}}\right)$
En déduire la constante de normalisation de la densité $f$.
Ecrire la fonction vraisemblance $L\left(\theta,\lambda\mid\underline{x}\right).$
Supposons que $\theta$ est connue.
1. En déduire une loi a priori conjuguée pour $\lambda$
2. Calculer la loi a posteriori de $\lambda$
3. Calculer l’estimateur de Bayes de $\lambda$ en considèrant la fonction de perte quadratique.
4. Reprendre les questions a,b et c. en considérant la loi a priori de Jeffrey.
On suppose maintenant que $\lambda$ est connue
1. En déduire une loi a priori conjuguée pour $\theta$
2. Calculer la loi a posteriori de $\theta$

Exercice

On considère le code suivant

>  p = seq(0.05, 0.95, by = 0.1) 
>  prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0) 
>  prior = prior/sum(prior) 
>  Likelyhood=function(p,s,f) p^s*(1-p)^f 
>  Lp=sapply(p,function(x) Likelyhood(x,11,16)) 
>  post=Lp*prior 
>  post=post/sum(post) 
>  round(post,4)  
 [1] 0.0000 0.0023 0.1291 0.4768 0.3338 0.0559 0.0021 0.0000 0.0000 0.0000

Quel est le modèle Bayésien considéré dans ce code?
Calculer l’estimateur bayèsien du param`etre considèré dans ce code.
Donner une loi a priori conjuguée du mode`ele Bayésien.

Modèle Bernoulli

$y_1,\ldots,y_n$ un $n-$échantillon de loi Bernoulli $\mathcal{B}(p)$.
La loi a priori conjuguée

\[g(p) = \displaystyle\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha-1}\, (1-p)^{\beta-1}\]

De moyenne

\[E(p)=\displaystyle\frac{\alpha}{\alpha+\beta}\]
De Variance

\[\mbox{Var}(p)=\displaystyle\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}\]

De densité

> curve(dbeta(x,4,3),from = 0,to = 1,col="red",lwd=2,xlab = "",ylab="")

Calcul d’une loi a priori : On croit que la médiane et 90ième percentile de $p$ sont respectivement égales à 0,3 et 0,5.

>  library(LearnBayes)
> quantile2=list(p=.9,x=.5)
> quantile1=list(p=.5,x=.3)
> beta.select(quantile1,quantile2)
[1] 3.26 7.19

La loi a posteriori \[g(p\mid \mbox{ données })\propto p ^{s+\alpha-1} (1-p)^{N-s+\beta-1}\]
Représentation des lois a priori, a posteriori et vraisemblance.

>  a  =  beta.select(quantile1,quantile2)[1]
>  b  =  beta.select(quantile1,quantile2)[2]
> s  =  11
> f  =  16
>  curve(dbeta(x,a+s,b+f),  from=0,  to=1,
+  xlab="p",ylab="Density",lty=1,lwd=4)
>   curve(dbeta(x,s+1,f+1),add=TRUE,lty=2,lwd=4)
>   curve(dbeta(x,a,b),add=TRUE,lty=3,lwd=4)
>   legend(.7,4,c("Prior","Likelihood","Posterior"),
+          lty=c(3,2,1),lwd=c(3,3,3))

Calcul de la probabilité $P(p\geq .5\mid \mbox{ données })$.

> 1 - pbeta(0.5,a+s,b+f)
[1] 0.0690226

Estimer une valeur future $\widetilde{y}$ (le nombre de fois de la réalisations d’événement de probabilité $p$ quand l’expérience s’est répétée $n$ fois). On calcule la densité prédictive de $\widetilde{y}$ :

\[f(\widetilde{y})=\int f(\widetilde{y}\mid p)g(p) dp.\]

Deux type densités prédictives : a priori si $g$ est la loi a priori ou a posteriori si $g$ est la loi a posteriori

Exemple

On lance une pièce de monnaie $n$ fois, on observe $y$ fois piles, on considère une loi a priori sur $p$. Supposons qu’on veut jeter $m$ fois la pièce et on veut prédire le nombre $\widetilde{y}$ d’apparitions de la face Pile.

Solution fréquentiste : estimer $p$ par $\widehat{p}=y/n$, et prédire $\widetilde{y}$ par $\widetilde{y}=m\widehat{p}=my/n$ ?
Solution bayésienne : calculer la densité prédictive a posteriori de $\widetilde{y}$.

\[\begin{array}{rcl} g(\widetilde{y}\mid y)&=&\int g(\widetilde{y},\, p\mid y) dp \\ &=& \int g(\widetilde{y}\mid \, p) \, g(p\mid y) dp \\ \end{array}\]

La loi Beta-Binomiale

On a déjà vu que si $p\sim \mbox{Beta}(a,b)$ alors $p\mid y\sim Beta(y+a,n-y+b)$. La loi densité prédictive a posteriori s’écrit

\[\begin{array}{rcl} g(\widetilde{y}\mid y)& = & \displaystyle\frac{1}{B(y+a,n-y+b)} \displaystyle\int_0^1 C_m^{\widetilde{y}} p^{\widetilde{y}}\\ & & (1-p)^{m-\widetilde{y}}\, p^{y+a}(1-p)^{n-y+b}\, dp\\ & & \\ &=& C_m^{\widetilde{y}} \displaystyle\frac{B(\widetilde{y}+y+a,m-\widetilde{y}+n-\widetilde{y}+b)}{B(y+a,n-y+b)} \end{array}\] On dit que $\widetilde{y}$ suit une loi Beta-Binomiale de paramètres $m,y+a,n-y+b$.

Application

On observe pour $n=40,$ $y=13$ et on veut prédire les probabilités d’observer le nombre de Piles après $m=13$

> library(emdbook)
>  y=13
>  n=40
>  ### a priori beta 2,4
>  a=2
>  b=4
>  m=13
>  yy=0:m
>  ## Calcul de la loi predictive a posteriori.
>  predy=sapply(yy,function(x) dbetabinom(x,shape1 = y+a,shape2 = n-y+b,size = m))
>  names(predy)=yy
>  plot(yy,predy,xlab="y",ylab="loi predictive a posteriori",
+  type="h",lwd=3)

Exemple : combien de fois on a jeté un dé

Supposons qu’on a jeté plusieurs fois une pi`{e}ce de monnaie équilibrée et on vous dit qu’on a obtenu ``face’’ 13 fois.
Question : Pourrions-nous estimer le nombre de jets de la pi`{e}ce ?
Solution : On considère un a priori non-informative sur $n$ : $g(n)\propto 1$ et la vraisemblance est alors \[L(n)\propto C_n^{13} \left( \frac{1}{2} \right)^n.\] Ainsi la loi a posteriori s’exprime \[g(n\mid x)\propto L(n) g(n)= C_n^{13} \left( \frac{1}{2} \right)^n.\]

>  n=13:100
>  posteriori=function(n,p) choose(n,13)*p^13*(1-p)^(n-13)
>  ### le cas p=.5
>  postn=sapply( n,function(n) posteriori(n,.5))
>  postn=postn/sum(postn)
>  #### Graphique et comparer avec  2*13=26.
>  plot(n,postn,type="h",lwd=.8,col="red",xlab="n",
+ ylab="a posteriori")
>  abline(v=26,col="blue",lwd=1.5)

Comparons $n$ avec $x/p$ pour diff'erentes valeurs de $p$.

>  p=seq(0.1,0.8,by=0.1)
>  par(mfrow=c(2,4))
>  for(i in 1:8){
+  postn=c()
+  postn=sapply( n,function(n) posteriori(n,p[i]))
+  postn=postn/sum(postn)
+  plot(n,postn,type="h",lwd=.8,col="red",xlab="n",
+  ylab="a posteriori", main=paste("p=",p[i],sep=""))
+  abline(v=13/p[i],col="blue",lwd=1.5)
+  }

Cours Statistique Bayésienne

Dhafer Malouche

3ième année ESSAI, 2015-2016

Introduction au calcul Bayésien

(Séance 1, 3/9/2016)

Théorème de Bayes

Exercice

Solution avec `R`

Question 1

Question 2

Solution théorique

Le modèle Bayésien et théorie de la décision

(Séance 2, 10/9/2015)

Notions de la théorie de décision.

Exemples : \(\Theta=\mathcal{E}=\mathbb{R}\)

Perte quadratique

Perte absolue

Exercice

Le cas bayésien.

Risque Bayésien

Exemples

Estimateur Bayésien

Exercice

Le choix de la loi a priori

(Séance 3, 17/9/2016)

Lois a priori conjuguées.

Exercice.

Lois a priori de Jeffrey.

Exercice

Exercice

Exercice

Modèle Bernoulli

Exemple

La loi Beta-Binomiale

Application

Exemple : combien de fois on a jeté un dé

Cours Statistique Bayésienne

Dhafer Malouche

3ième année ESSAI, 2015-2016

Introduction au calcul Bayésien

(Séance 1, 3/9/2016)

Théorème de Bayes

Exercice

Solution avec R

Question 1

Question 2

Solution théorique

Le modèle Bayésien et théorie de la décision

(Séance 2, 10/9/2015)

Notions de la théorie de décision.

Exemples : \(\Theta=\mathcal{E}=\mathbb{R}\)

Perte quadratique

Perte absolue

Exercice

Le cas bayésien.

Risque Bayésien

Exemples

Estimateur Bayésien

Exercice

Le choix de la loi a priori

(Séance 3, 17/9/2016)

Lois a priori conjuguées.

Exercice.

Lois a priori de Jeffrey.

Exercice

Exercice

Exercice

Modèle Bernoulli

Exemple

La loi Beta-Binomiale

Application

Exemple : combien de fois on a jeté un dé

Solution avec `R`