Premiers pas avec OpenBUGS

Exemple : lancement de satellites

• Entre 1980 et 2002 il y a eu 11 lancements de satellites, 3 succés et 8 checs.

• Objectif : Estimer la probabilité de réussite d’un nouveau lanceur, fournir un intervalle de confiance/crédibilité.

Approche fréquentiste (Rappel)

• Distribution des données (le modèle statistique)

  • \(y\) est le nombre de succès parmi \(n\) lancements de satelittes,
  • \(\theta\) est la probabilité de succès d’un lancement \[y\mid \theta\;\sim\; \mathcal{B}(n,\theta);\]

• Vraisemblance et estimateur du maximum de vraisemblance

  • \(l(y\mid \theta)=C_n^y \theta^y (1-\theta)^{n-y}\), \(\log\left( l(y\mid \theta)\right)=y \log \theta +(n-y) \log (1-\theta)\).

  • \(\widehat{\theta}=\mbox{argmax}_\theta \log\left( l(y\mid \theta)\right)\)

  \[\widehat{\theta} =\displaystyle\frac{y}{n}\]

  • Calcul d’intervalle de confiance (asymptotique \(n\) grand) : \((\widehat{\theta}\pm z_{\alpha/2} \sqrt{\displaystyle\frac{1}{n}\widehat{\theta}(1-\widehat{\theta})}\).

Approche Bayésienne

• On considère une loi a priori Beta sur \(\theta\) \[\pi(\theta)\propto \theta^{a-1}(1-\theta)^{b-1}.\]

• Calcul de la loi a posteriori \[\begin{array}{rcl}\pi(\theta\mid y)&\propto &l(y\mid \theta) \pi(\theta) & \propto & \theta^{y+a-1}(1-\theta)^{n-y+b-1}.\end{array}\] Donc \(y\mid \theta\sim \mbox{Beta}(y+a,n-y+b)\).

• La distribution, moyenne, variance, intervalle de crédibilité a posteriori ?

• Solution OpenBuGS

OpenBugs

• Présentation OpenBUGS fait partie du projet BUGS (Bayesian inference Using Gibbs Sampler) qui vise à rendre simple la pratique des méthodes MCMC aux statisticiens. Il a été développé par l’unité MRC Biostatistics de l’université de Cambridge.

• OpenBUGS est un logiciel libre et gratuit.

• Utilisation possible

  • Via une interface clique-bouton qui permet de contrôler l’analyse,
  • En utilisant des modèles définis par des interfaces graphiques, appelés DoddleBUGS,
  • Via d’autres logiciels tels que R (en particulier via le package R2OpenBUGS)

Algorithme de Gibbs

• Supposons que le vecteur paramètre \(\theta\) est constitué de \(k\) composantes, i.e. \(\theta^T=(\theta_1,\ldots,\theta_k)\). Les étapes de l’algorithme de Gibbs sont :

  1. Choisir des valeurs initiales pour \(\theta_1^{(0)},\ldots,\theta_k^{(0)}\).
  2. Générer \(\theta_1^{(1)}\) selon \(\pi\left(\theta_1\mid \theta_2^{(0)},\ldots,\theta_k^{(0)},y\right)\).
     
  3. Générer \(\theta_2^{(1)}\) selon \(\pi\left(\theta_2\mid \theta_1^{(1)},\theta_3^{(0)}\ldots,\theta_k^{(0)},y\right)\).
     …
  4. Générer \(\theta_k^{(1)}\) selon \(\pi\left(\theta_k\mid \theta_1^{(1)},\theta_2^{(1)}\ldots,\theta_{k-1}^{(1)},y\right)\).

Répéter l’étape 2 plusieurs milliers de fois.

Exemple de code OpenBUGS

Le modèle :

model{
for(i in 1:n){
y[i] ~ dbern(p)
}
p~dbeta(a,b)
}

Les données

> list(n=11,y=c(0,0,0,0,0,0,1,0,0,1,0),a=2,b=3)

Les valeurs initiales.

> list(p=.2)
> list(p=.9)
> list(p=.4)

Les étapes de l’exécution sous OpenBUGS

1. Vérifier le modèle,
2. Charger les données,
3. Spécifier le nombre de chaînes MCMC,
4. Compiler le modèle,
5. Charger les conditions initiales,
6. Générer les réalisations ``burnin’’,
7. Spécifier les paramètres/quantités àconserver,
8. Générer les réalisations àconserver,
9. Vérifier la convergence et mise en place des résultats.

D’autres Exemples:

Jet d’une pièce de monnaie.

On jette une pièce de monnaie 8 fois dont la probablité d’avoir Pile est de 0.5. Quelle est la probabilitée de voir Pile au plus 2 fois.

Le modèle mathématique

\[Y\sim\mbox{Binomial}(0.5,8)\] Il s’agit d’estimer \(\mathbb{P}(Y\leq 2)\).

Le modèle sous OpenBugs

model{
  Y ~ dbin(0.5,8)  
  P <- step(2.5-Y)
}

Comment faire fonctionner ce code dans OpenBUGS

1. Ouvrir une nouvelle fen^etre dans OpenBUGS, File>New
2. Copier/Coller le code ci-dessus.
3. Model>Specification
4. Select model et clique sur check model
5. clique sur compile ensuite gen inits.
6. Model>Update et clique sur update
7. Inference> Samples et introduit le paramètre P dans set
8. Tous les boutons de la denrière interface vont s’allumer, donc vous pouvez avoir votre estimation.

En effet, la fonction step va prendre la valeur si2.5-Yest $\geq 0$ et 0 sinon. DoncPest éegal à siYest plus petit o\u éegal à 1.

Utilisation de OpenBUGS pour simuler une distribution

Une loi de Student

Supposons qu’on voudrait simuler une loi de student de moyenne 10, de préecision 2 et de ddl 4.

model{
y ~ dt(10,2,4)
}

Une transformation de la loi Normale

On considère \(Z\sim\mathcal{N}(0,1)\), on veut simuler \(Y=(2Z+1)^3\) et estimer \(\mathbb{P}(Y>10)\)

model{
  z ~ dnorm(0,1)
  y <- pow(2*z+1,3)
  p <- step(y-10)
}

Combien de réparations?

Supposons que le coût de réeparation a une distribution Gamma de moyenne 100 TND et d’écart-type 50 TND. On veut répondre à la question, combien de fois peut-on réparer avec 1000 TND?

1. Montrer qu’une Gamma(4,0.04) a une moyenne 100 et un éecart-type 50 TND
2. Ecrire le code OpenBUGS permettant de calculer ce nombre.

On considère \(Y_i\), \(i=1,\ldots,I\) variables aléeatoires de loi Gamma(4,0.04) et il s’agit de chercher \(M\) tel que \[ \displaystyle\sum_{i=1}^M Y_i <1000\].

Le code OpenBUGS

On prendra \(I=20\), on va calculer la probabilité de faire plus de 12 réparations

model{
  for(i in 1:20) {
    Y[i] ~ dgamma(4,0.04)
  }
  cum[1] <- Y[1]
  for(i in 2:20) {
    cum[i] <- cum[i-1]+Y[i]
  }
  for(i in 1:20) {
    cum.step[i] <- i*step(1000- cum[i])
  }
  number <- ranked(cum.step[],20) ## le maximum dans cum.step
  P <- step(number-12)
}

Rappelons que

• rank(v,s) s’applique au vecteur v et et au nombre s et calcule le nombre composantes dans v qui sont plus petites que s

Utilisation de OpenBUG à partir de R

• Pour utiliser OpenBUGS il faut installer dans R le package R2OpenBUGs.

• Considérons le problème suivant: On jette deux pièces de monnaie dont ne connaît pas la probbilité d’avoir Pile. Les deux pièces sont jetées 86 fois. Pour la première on a obtenu 83 fois Pile et pour la deuxième on obtenu 72 fois Pile. Comparer les deux probabilités d’avoir Pile.

• D’abord écrire les données dans R.

> donnees=list(r=c(83,72),n=c(86,86))

• Le modèle sous OpenBUGS

model{
  for(i in 1:2){
    r[i] ~ dbin(p[i],n[i])
  }
  for (i in 1:2){ p[i] ~ dunif(0,1)}
  delta <- p[1] - p[2]
  delta.up <- step(delta)
  lambda <- log( (p[1]/(1-p[1])) / (p[2]/(1-p[2])) );
  lambda.up <- step(lambda)
}

• Le code pour générer le modèle à partir de R sans avoir recours à OpenBUGS:

On crée un fichier nommé exempleBinom.txt qui contient le code BUGS

> sink('exempleBinom.txt')
> cat("
+ model{
+   for(i in 1:2){
+     r[i] ~ dbin(p[i],n[i])
+   }
+   for (i in 1:2){ p[i] ~ dunif(0,1)}
+   delta <- p[1] - p[2]
+   delta.up <- step(delta)
+   lambda <- log( (p[1]/(1-p[1])) / (p[2]/(1-p[2])) );
+   lambda.up <- step(lambda)
+ }
+ ",fill=T)
> sink()

On crée une variable avec le nom du fichier.

> filename<-"exempleBinom.txt"

Les valeurs initiales

> inits<-function(){
+   inits=list(p=c(.5,.5))
+ }

Les paramètres à estimer

> params <- c('p','delta.up','lambda.up')

Enfin l’exécution sous R

> library(R2OpenBUGS)
> outBinom <-bugs(donnees,inits,params,filename,codaPkg=F,n.thin =1,
+                n.iter=10000,debug=F,n.chains = 3,
+                working.directory=getwd())

Le résultat

> outBinom
Inference for Bugs model at "exempleBinom.txt", 
Current: 3 chains, each with 10000 iterations (first 5000 discarded)
Cumulative: n.sims = 15000 iterations saved
          mean sd 2.5% 25% 50%  75% 97.5% Rhat n.eff
p[1]       1.0  0  0.9 0.9 1.0  1.0   1.0    1 15000
p[2]       0.8  0  0.7 0.8 0.8  0.9   0.9    1 15000
delta.up   1.0  0  1.0 1.0 1.0  1.0   1.0    1 15000
lambda.up  1.0  0  1.0 1.0 1.0  1.0   1.0    1 15000
deviance   9.3  2  7.3 7.9 8.7 10.1  14.7    1 15000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 1.8 and DIC = 11.0
DIC is an estimate of expected predictive error (lower deviance is better).

• Notons que le paramètre lambda est le logarithme du odds ratio, Si lambda est :

  • proche de 0, on a la même chance d’avoir Pile avec les deux pièces. ;
  • supérieur à 0, il est plus fréquente d’avoir Pile avec la pièce 1 qu’avec la pièce 2 ;
  • bien supérieur à 0, il est beaucoup plus fréquente d’avoir Pile avec la pièce 1 ququ’avec la pièce 2 ;
  • inférieur à 0, il est moins fréquente d’avoir Pile avec la pièce 1 qu’avec la pièce 2 ;
  • proche de -\(\infty\), il est beaucoup moins fréquente d’avoir Pile avec la pièce 1 qu’avec la pièce 2 ;

• Les paramètres estimés dans l’objet outBinom

  • deviance : Si \(y=(y_1,\ldots,y_n)\) est le vecteur des données dont on suppose généré par une loi de probabilité de densité \(f(y\mid \theta)\). Alors la fonction vraisemblance \[l(y\mid \theta)=\displaystyle\prod_{i=1}^n f(y_i\mid \theta)\] Donc la deviance s’exprime par \[D(\theta)=-\log l(y\mid \theta)\] Dans notre cas : \(y_1\sim\)Binom\((N,p_1)\) et \(y_2\sim\)Binom\((N,p_2)\), donc \[l(y\mid p_1,p_2)\propto \left(\displaystyle\frac{p_1}{1-p_1}\right)^{y_1}\,\left(\displaystyle\frac{p_2}{1-p_2}\right)^{y_2}\]

  • DIC : Deviance Information Criteria
    • La moyenne \(\bar{D}=\mathbf{E}^\theta[D(\theta)]\) est une mesure de l’ajustement du modèle aux données; plus \(\bar{D}\) est grand, plus l’ajustement est mauvais.
    • pD est une estimation du nombre effective des paramètres du modèle : \(pD=\bar{D}-D(\bar{\theta})\), plus les valeurs du \(p_D\) est élévé, plus le modèle ajuste au mieux les données. C’est pour cela qu’on définit une pénalisation à la vraisemblance.
      La Deviance Information Criteria se calcule \[DIC=pD+\bar{D}=D(\bar{\theta})+2pD\] Des modèles avec des valeurs faibles du DIC sont les plus préférés. Comme \(\bar{D}\) décroit quand le nombre de paramètres dans le modèle croit, le \(p_D\) va compenser cette effet en favorisant les modèles avec un petit nombre de paramètres.

• Estimation des densités a posteriori

Les probabilités

> library(reshape2)
> dt=melt(outBinom$sims.matrix[,1:2])
> library(ggplot2)
> gr<-ggplot(data=dt,aes(x=value,fill=Var2,colour=Var2))+geom_density(alpha=.5)
> gr

delta.up

> dt=data.frame(outBinom$sims.matrix[,3:4])
> gr<-ggplot(data=dt,aes(x=delta.up))+geom_density(alpha=.5,fill="pink")
> gr

lambda.up

> gr<-ggplot(data=dt,aes(x=lambda.up))+geom_density(alpha=.5,fill="pink")
> gr

Exercice : Comparaison de deux modèles bayésiens

On considère les modèles bayésiens suivants :

• Modèle M1 : \(y\sim\mbox{Binomial}(n,\tau)\) où \(\tau\sim \mbox{Beta}(a,b)\)

• Modèle M2 :\(y\sim\mbox{Binomial}(n,\tau)\) où \(\tau=\displaystyle\frac{1}{p+q+1}\) où \(p,q\sim \mbox{Uniform}[0,1]\).

Ecrire les deux codes BUGS et R pour générer la loi a posteriori des modèles bayésiens M1 et M2. Comparer les DIC, pD et la déviance des deux modèles, Conclure?