Packages R nécessaires

L’execution des codes R ci-après fait appel aux packages tidyverse,rjags,ggextra,knitr,ggpubr,purrr, grid à installer préalablement

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
rm(list=ls())
library(tidyverse)
## ── Attaching packages ───────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ──────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(rjags)
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
library(ggExtra)

Lecture des données

On commence par introduire les données des trois groupes sous la forme d’une matrice \(X\) avec, en ligne le groupe d’étudiants et en colonne les résultats pour la première pêche (\(C_1\)), la quantité de poissons pris lors de la seconde pêche, mais pas la première \(C_{20}\) et enfin la quantité de poissons capturés les deux fois

X <- matrix(c(99,56,62,116,56,25,82,74,39),3,3, byrow = TRUE)
rownames(X)  <- c("G1","G2","G3")
colnames(X)<- c("C1","C20","C21")
print(X,digits=1)
##     C1 C20 C21
## G1  99  56  62
## G2 116  56  25
## G3  82  74  39

Comparaisons des deux estimateurs

On construit d’abord les fonctions petersen et schnabel qui évaluent les estimeurs de Lincoln-Petersen et de Schnabel-Chapman ainsi que leur écart-type conditionnel à \(C_1, C_2\) si l’argument but vaut ESTIM. Pour Lincoln-Petersen, on utilise que \(C_{21}|C_1,C_2\) est une hypergéométrique, déjà codée en R, ce qui permet sa simulation (1000 tirages). Pour Schnabel-Chapman, on calcule l’estimateur de la variance de Wittes.

petersen<-function(C1=99,C20 = 56, C21=62, but="ESTIM"){
  C2=C20+C21
  nu.hat <- C1*C2/C21
  pi.hat <- (C1+C2)/(2*nu.hat)
  if(but=="ESTIM"){
  sig.nu.hat= var(C1*C2/rhyper(nn = 1000,m = C1,n = nu.hat-C1,k = C2))^0.5
  return(list(pi.hat=pi.hat,nu025.hat=nu.hat-2*sig.nu.hat,
              nu.hat=nu.hat,nu975.hat=nu.hat+2*sig.nu.hat,
              sig.nu.hat=sig.nu.hat))
  } else   return(list(pi.hat=pi.hat, nu.hat=nu.hat))

}
##__________________________________________
schnabel<-function(C1=99,C20 = 56, C21=62, but="ESTIM"){
  C2=C20+C21
  nu.hat <- ((C1+1)*(C2+1)/(C21+1))-1
  pi.hat <- (C1+C2)/(2*nu.hat)
  if(but=="ESTIM"){
  sig.nu.hat=sqrt((C1+1)*(C20+C21+1)*(C1-C21)*(C2-C21)/(C21+1)/(C21+1)/(C21+2))
  return(list(pi.hat=pi.hat,nu025.hat=nu.hat-2*sig.nu.hat,
              nu.hat=nu.hat,nu975.hat=nu.hat+2*sig.nu.hat,
              sig.nu.hat=sig.nu.hat))
  } else   return(list(pi.hat=pi.hat, nu.hat=nu.hat)) 
}

Pour refaire la figure 6, on a juste besoin de simulations pour évaluer la variabilité d’échantillonnage des estimateurs.

nu=200
pi1=pi2=0.1 
repet=100000
C1=rbinom(repet,nu,pi1)
C21=rbinom(repet,C1,pi2)
C20=rbinom(repet,nu-C1,pi2)
NchapLincoln=petersen(C1 = C1, C20 = C20,C21 = C21,but = "Histo")$nu.hat
NchapSchnabel=schnabel(C1 = C1, C20 = C20,C21 = C21,but = "Histo")$nu.hat
print(paste("l'estimateur de Peterson-Lincoln est non défini avec une fréquence de ",round(sum(is.infinite(NchapLincoln))/repet,2)))
## [1] "l'estimateur de Peterson-Lincoln est non défini avec une fréquence de  0.13"
print(paste("Moyennes de nu, Lincoln et Schnabel"))
## [1] "Moyennes de nu, Lincoln et Schnabel"
print(c(nu, mean(NchapLincoln[!is.infinite(NchapLincoln)],na.rm=T),
        mean(NchapSchnabel,na.rm=T)),2)
## [1] 200 225 178
print(paste("Risque quadratique de Lincoln et Schnabel"))
## [1] "Risque quadratique de Lincoln et Schnabel"
print(c(mean((NchapLincoln[!is.infinite(NchapLincoln)]-nu)^2,na.rm=T),
        mean((NchapSchnabel-nu)^2,na.rm=T)), dig=2 )
## [1] 15127 10027

On peut alors réaliser une représentation de la répartition empirique des deux estimateurs

d<-tibble(`Lincoln-Petersen`=NchapLincoln,`Schnabel-Chapman`=NchapSchnabel)
d %>% gather(key=estimateurs,value = valeurs) %>% ggplot()+geom_density(aes(x=valeurs, color=estimateurs),adjust=2.5)+labs(y="densité empirique")+geom_vline(xintercept=nu)

Dans le même esprit, on peut refaire la figure 7 de l’article illustrant la décroissance du biais relatif avec \(nu\).

nu <- seq(100, 1000, by=50)
repet <- 20000
pi1 <- 0.3 ; pi2 <- 0.3
biaisP <- numeric(length(nu)) ; 
biaisS <- numeric(length(nu))
for (k in 1:length(nu)) {
  C1 <- rbinom(repet, nu[k], pi1)
  C21 <- rbinom(repet, C1, pi2)
  C20 <- rbinom(repet, nu[k]-C1, pi2)
  C2 <- C21 + C20
  NchapLincoln <- C1*C2/C21
  NchapSchnabel <- (C1+1)*(C2+1)/(C21+1)-1
  biaisP[k] <- (mean(NchapLincoln[C21!=0])-nu[k])/nu[k]
  biaisS[k] <- (mean(NchapSchnabel[C21!=0])-nu[k])/nu[k]
}
d<-tibble(`Lincoln-Petersen`=biaisP,`Schnabel-Chapman`=biaisS, nu=nu)
d %>% gather(key=Estimateur,value = valeurs,-nu) %>% ggplot()+geom_line(aes(x=nu, y= valeurs,color=Estimateur))+labs(y="biais relatif")+geom_hline(yintercept=0)

Pour le calcul du tableau 3, on va utiliser la fonction map du package purr pour répéter le calcul sur chaque groupe.

d_grouped <- (X) %>% data.frame() %>% bind_cols(Gpe=c("G1","G2","G3")) %>% 
  group_by(Gpe) %>% nest 
d_grouped<- d_grouped %>% 
  mutate( petersen=purrr::map(data, ~{petersen(.x$C1,.x$C20,.x$C21)} ),
          schnabel=purrr::map(data, ~{schnabel(.x$C1,.x$C20,.x$C21)} ),
         ) 
results<-d_grouped %>%
  mutate(nu.p=map_dbl(petersen,"nu.hat"),
         sig.nu.p=map_dbl(petersen,"sig.nu.hat"),
         nu.s=round(map_dbl(schnabel,"nu.hat")),
         sig.nu.s=map_dbl(schnabel,"sig.nu.hat")
         ) %>% 
  select(-petersen,-schnabel,-data)

knitr::kable(results, format="markdown", digits=c(0,0,1,0,1))
Gpe nu.p sig.nu.p nu.s sig.nu.s
G1 188 10.1 188 9.9
G2 376 59.1 368 51.8
G3 238 23.8 236 21.4 
#knitr::kable(results, format="latex", digits=c(0,0,1,0,1))

Mise en évidence de la composante hypergéométrique de la vraisemblance du modèle CMR

Un détour par la loi conditionnelle de \(C_{21}\) sachant \(C_{1}\) et \(C_{2}\) permet d’introduire la loi hypergéométrique comme autre description possible des résultats d’une expérience de CMR basée sur deux pêches successives.

Commençons par écrire la loi d’échantillonnage selon le modèle binomial séquentiel: \[\begin{align*} \lbrack C_{1},C_{20},C_{21}|\nu,\pi_1,\pi_2] & = [C_{1}|\nu, \pi_1][C_{20}|C_{1},\nu, \pi_2][C_{21}|C_{1},\pi_2]\\ & =\frac{\pi_{1}^{C_{1}}(1-\pi_{1})^{\nu-C_{1}}\nu!}{(\nu-C_{1}% )!C_{1}!}\\ & \times\frac{\pi_{2}^{C_{20}}(1-\pi_{2})^{\nu-C_{1}-C_{20}}(\nu -C_{1})!}{(\nu-C_{1}-C_{20})!C_{20}!}\\ & \times\frac{\pi_{2}^{C_{21}}(1-\pi_{2})^{C_{1}-C_{21}}C_{1}!}% {(C_{1}-C_{21})!C_{21}!}% \end{align*}\]

On montre d’abord que la loi de \(C_{1},C_{20},C_{21}|\nu,\pi_1,\pi_2\) est la même en adoptant le modèle multinomial de bilan, décrit dans la section 4 équation 1 par: \[C_{1}|\nu, \pi_1 \sim dbin(\pi_{1},\nu)\\ C_{20}|C_{1}, \nu, \pi_2 \sim dbin(\pi_{2},\nu-C_{1})\\ C_{21}|C_{1}, \pi_2 \sim dbin(\pi_{2},C_{1})\] En effet, pour le modèle multinomial à 4 catégories, il vient: \[\begin{align*} \lbrack C_{1},C_{20},C_{21}|\nu,\pi_1,\pi_2] & =\frac{\pi_{1}^{C_{1}}(1-\pi_{1}% )^{\nu-C_{1}}\nu!}{(\nu-C_{1})!C_{1}!}\times\frac{\pi_{2}^{C_{20}}% (1-\pi_{2})^{\nu-C_{1}-C_{20}}(\nu-C_{1})!}{(\nu-C_{1}-C_{20})!C_{20}!}\\ & \times\frac{\pi_{2}^{C_{21}}(1-\pi_{2})^{C_{1}-C_{21}}C_{1}!}% {(C_{1}-C_{21})!C_{21}!}\\ & =\frac{\nu!}{C_{21}!(C_{1}-C_{21})!C_{20}!(\nu-C_{1}-C_{20})!}\\ & \times\pi_{1}^{C_{1}-C_{21}+C_{21}}(1-\pi_{1})^{\nu-C_{1}% -C_{20}+C_{20}}\pi_{2}^{C_{20}+C_{21}}(1-\pi_{2})^{C_{1}-C_{1}% +\nu-C_{21}-C_{20}}\\ & \propto\left( \pi_{1}\pi_{2}\right) ^{C_{21}}\left( \pi _{1}(1-\pi_{2})\right) ^{C_{1}-C_{21}}\\ & \times((1-\pi_{1})\pi_{2})^{C_{20}}((1-\pi_{1})(1-\pi _{2}))^{\nu-C_{1}-C_{20}}% \end{align*}\]

Pour plus de clarté, on omet momentanément le conditionnement sur \(\nu\), \(\pi_1\) et \(\pi_2\). La loi conditionnelle de \(C_{21}\) sachant \(C_{1}\) et \(C_{2}\) est facilement mise en évidence, en montrant comment la lire à partir de la conjointe de \(C_{1}\), \(C_{20}\) et \(C_{21}\) : \[\begin{align*} \lbrack C_{1},C_{20},C_{21}|\nu,\pi_1,\pi_2] & =[C_{1},C_{2}-C_{21},C_{21}|\nu,\pi_1,\pi_2]\\ \lbrack C_{21}|C_{1},C_{2},\nu,\pi_1,\pi_2] & \propto\frac{1}{C_{21}!(C_{1}-C_{21} )!C_{20}!(\nu-C_{1}-C_{20})!}\\ & \propto\frac{1}{C_{21}!(C_{1}-C_{21})!(C_{2}-C_{21})!(\nu-C_{1} -C_{2}+C_{21})!}\\ & \propto\frac{C_{1}!(\nu-C_{1})!}{C_{21}!(C_{1}-C_{21})!(C_{2}-C_{21} )!(\nu-C_{1}-C_{2}+C_{21})!}\\ & \propto C_{C_{1}}^{C_{21}}C_{\nu-C_{1}}^{C_{2}-C_{21}}=\frac{C_{C_{1} }^{C_{21}}C_{\nu-C_{1}}^{C_{2}-C_{21}}}{C_{\nu}^{C_{2}}}\\ \end{align*}\] car \({\displaystyle\sum\limits_{y}}C_{N}^{y}C_{B}^{K-y}=C_{N+B}^{K}\) d’après la formule du binome . On reconnaît une loi hypergéométrique de paramètres \((\nu,C_{1},C_{2})\). En utilisant la décomposition en produits suivante \([C_1,C_2,C_{21}|\nu,\pi_1,\pi_2]=[C_1|\nu,\pi_1][C_2|\nu,\pi_2][C_{21}|C_1,C_2,\nu]\), on tombe finalement sur la modélisation traditionnelle alternative d’une expérience de CMR, fondée sur un tirage sans remise des poissons recapturés lors de la 2ème pêche via une loi hypergéométrique:

\[\begin{align*} C_{1}|\nu,\pi_1 & \sim dbin(\pi_{1},\nu)\\ C_{2}|\nu,\pi_2 & \sim dbin(\pi_{2},\nu)\\ C_{21}|C_{1},C_{2},\nu & % \sim hypergeometrique(C_{2},C_{1},\nu-C_{1}) \sim hypergeometrique(\nu,C_{1},C_{2}) \end{align*}\]

Estimateur du maximum de vraisemblance

Sur le même modèle que les fonctions petersen et schnabel, on construit une fonction qui réalise l’estimation du max de vraisemblance et qui construit un intervalle de confiance à 95% par approximation du chi-deux sur la vraisemblance profilée.

##__________________________________________
loglikHaricot <-function(nu,pi,C1,C21,C20){
  dbinom(C1,nu,pi,log=T)+
    dbinom(C21,C1,pi,log=T)+
    dbinom(C20,nu-C1,pi,log=T)
}
##__________________________________________
maxloglikHaricot <-function(C1=99,C20 = 56, C21=62, numax=1000){
  nutest= seq(C1+C20 , numax, by=1)
  C2=C20+C21
  logliktest= loglikHaricot(nu=nutest,pi = (C1+C2)/(2*nutest) ,
                            C1=C1,C21=C21,C20=C20)
  imax=which.max(logliktest)
  nu.hat=nutest[imax]
  pi.hat=(C1+C2)/(2*nu.hat)
  
  #likelihood ratio test
  irange = 2*(logliktest[imax]-logliktest)-qchisq(p=0.95,df=1)
  i1=which.min(irange^2)
  if(i1 < imax) {i2 = which.min(irange[imax:length(irange)]^2)+(imax-1)
  nu025.hat = nutest[i1]
  nu975.hat = nutest[i2] } else {
    # i2 correspond à q05
    i2=which.min(irange[1:imax]^2) 
    nu025.hat = nutest[i2]
    nu975.hat = nutest[i1] }
    return(list(pi.hat=pi.hat,nu025.hat=nu025.hat,
              nu.hat=nu.hat,nu975.hat=nu975.hat))
  
}

Le tableau 4 est alors obtenu en appelant les fonctions petersen, schnabel et maxloglikHaricot , puis en rangeant les résultats dans un tableau imprimé joliment grâce à la fonction R kable.

d_grouped <- (X) %>% data.frame() %>% bind_cols(Gpe=c("G1","G2","G3")) %>% 
  group_by(Gpe) %>% nest 
d_grouped<- d_grouped %>% 
  mutate( petersen=purrr::map(data, ~{petersen(.x$C1,.x$C20,.x$C21)} ),
          schnabel=purrr::map(data, ~{schnabel(.x$C1,.x$C20,.x$C21)} ),
          maxloglik = purrr::map(data, ~{maxloglikHaricot(.x$C1,.x$C20,.x$C21, numax=1000)} )) 
results<-d_grouped %>%
  mutate(#pi.p= map_dbl(petersen,"pi.hat"),
         nu025.p=round(map_dbl(petersen,"nu025.hat")),
         nu.p=round(map_dbl(petersen,"nu.hat")),
         nu975.p=round(map_dbl(petersen,"nu975.hat")),
         #pi.s= map_dbl(schnabel,"pi.hat"),
         nu025.s=round(map_dbl(schnabel,"nu025.hat")),
         nu.s=round(map_dbl(schnabel,"nu.hat")),
         nu975.s=round(map_dbl(schnabel,"nu975.hat")),
         #pi.ml= map_dbl(maxloglik,"pi.hat"),
         nu025.ml=round(map_dbl(maxloglik,"nu025.hat")),
         nu.ml=round(map_dbl(maxloglik,"nu.hat")),
         nu975.ml=round(map_dbl(maxloglik,"nu975.hat")),
           ) %>% 
  select(-petersen,-schnabel,-maxloglik,-data)

knitr::kable(results, format="markdown")
Gpe nu025.p nu.p nu975.p nu025.s nu.s nu975.s nu025.ml nu.ml nu975.ml
G1 168 188 209 168 188 208 172 189 213
G2 251 376 501 264 368 472 297 385 534
G3 191 238 284 193 236 278 205 242 299 
#knitr::kable(results, format="latex")

La piste bayésienne

2 pêches

On réalise l’inférence bayésienne en profitant de la commodité du langage BUGS (on déclare le modèle dans la variable model). On utilise une loi a priori bêta(2,1) pour \(\pi\), un choix acceptable pour encoder le jugement d’un étudiant annonçant que son évaluation moyenne personnelle de l’efficacité \(\pi\) est de l’ordre de \(2/3\) mais qu’il n’en est pas très sûr : environ une chance sur deux pour l’intervalle \([0.4,0.8]\) ou son complémentaire. Et une loi a priori uniforme pour \(\nu\) , dans un souci de traduire une certaine équiprobabilité, mais sur un intervalle limité supérieurement à \(1000\) haricots .

numax=1000
model<-"model
{
#priors
   nuReal ~ dunif(0,numax)
   nu_prior ~ dunif(0,numax)
   nu <- round(nuReal)
   pi ~ dbeta(2,1)
   pi_prior ~ dbeta(2,1)
#Vraisemblance
   C1 ~ dbin(pi,nu)
   PasPris <- nu - C1
   C21 ~ dbin(pi,C1)
   C20 ~ dbin(pi,PasPris)
}
"
C1=99; C20=56;  C21=62;
data_jags<-list(C1=C1,C20=C20, C21=C21, numax=numax)
params <- c("pi", "pi_prior", "nu", "nu_prior")
min_nu <- C1 + C20 + 1
inits <- function(){
  list(pi = rbeta(1, 1, 1), pi_prior = rbeta(1, 1, 1),
       nuReal = runif(1, min_nu, numax))
}
n_iter <- 25000 # TOTAL Number of iterations including Burn in
n_burn <- 5000 #Number of iterations for Burn in
n_thin <- 10 #thin
n_chains <- 3 #nb chaines
modele <- jags.model(file = textConnection(object = model),
                     data = data_jags, 
                     inits = list(inits(),inits(),inits()), 
                     n.adapt =  n_burn,
                     n.chains = n_chains, quiet = FALSE)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 3
##    Unobserved stochastic nodes: 4
##    Total graph size: 13
## 
## Initializing model
update(modele,n_burn)

Un point de vigilance important : Bien que de nombreux outils d’inférence bayésienne, relativement simples d’utilisation, existent aujourd’hui, ces outils n’intègrent aucun garde-fou permettant de s’assurer de leur bonne utilisation pratique : un minimum de connaissances en statistique bayésienne est donc requise pour une utilisation sereine, e.g., comment s’assurer, dans la mesure du possible, du bon comportement de l’algorithme d’inférence bayésienne mis en oeuvre? C’est l’objet du package coda qui calcule la taille effective d’échantillons MCMC, et met en oeuvre le test de Gelman-Rubin pour s’assurer que les chaînes ont convergé.

echantillons_mcmc <- coda.samples(model = modele, 
                                       variable.names = params,
                                        n.iter = n_iter, 
                                        thin = n_thin)
# Have a look on posterior statistics
summary(echantillons_mcmc) 
## 
## Iterations = 10010:35000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 2500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##              Mean       SD  Naive SE Time-series SE
## nu       190.8473  10.5239 0.1215196      0.1272015
## nu_prior 494.1731 290.7658 3.3574743      3.4358545
## pi         0.5703   0.0398 0.0004596      0.0004648
## pi_prior   0.6664   0.2347 0.0027103      0.0027083
## 
## 2. Quantiles for each variable:
## 
##              2.5%      25%      50%      75%    97.5%
## nu       173.0000 183.0000 190.0000 197.0000 214.0000
## nu_prior  22.3689 243.4486 490.5462 744.0597 975.3858
## pi         0.4899   0.5437   0.5714   0.5977   0.6454
## pi_prior   0.1503   0.5048   0.7074   0.8626   0.9874
# Have a look on MCMC chains
gelman.plot(echantillons_mcmc[, params])

gelman.diag(echantillons_mcmc[, params])
## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## pi                1       1.00
## pi_prior          1       1.00
## nu                1       1.00
## nu_prior          1       1.01
## 
## Multivariate psrf
## 
## 1
effectiveSize(echantillons_mcmc[, params])
##       pi pi_prior       nu nu_prior 
## 7346.520 7510.693 6858.606 7165.587

On stocke les résultats dans la variable echantillons_posterior sous forme de tableau et on produit la figure 8.

echantillons_posterior <- echantillons_mcmc %>% 
  map_dfr(# On applique aux 3 éléments de la liste la même function
    function(tab_mcmc){ # tab mcmc
      resultat <- tab_mcmc %>% # On prend le tableau mcmc
        as.data.frame() %>% # On transforme en data.frame
        mutate(num_echantillon = 1:nrow(tab_mcmc)) # On créé un indice de numéro
    },
  .id = "chaine") %>% # On créée une colonne chaine qui stocke le numero de chaine
  as_tibble()


p <- echantillons_posterior %>% rename(x=pi,y=nu) %>% 
  ggplot(aes(x = x,y=y)) +
  geom_point(col = "red") + 
  geom_density_2d(col = "black") +
  theme(legend.position="none")+
  labs(x=expression(pi),
       y=expression(nu),title="CMR: Population vs capturabilité")
ggExtra::ggMarginal(p, type = "densigram")

Et le tableau 5

sommaire<-echantillons_posterior %>% mutate(chaine=1:length(chaine)) %>% map_df(~{tibble(
                                     mean=mean(.x),
                                     sdt=sd(.x),
                                     q05=quantile(.x,probs=.05),
                                     q25=quantile(.x,probs=.25),
                                     q50=quantile(.x,probs=.5),
                                     q75=quantile(.x,probs=.75),
                                     q95=quantile(.x,probs=.95)
                                                )},
                                  .id="param") 
knitr::kable(sommaire[-1,], format="markdown",digits=2)
param mean sdt q05 q25 q50 q75 q95
nu 190.85 10.52 176.00 183.00 190.00 197.00 210.00
nu_prior 494.17 290.77 45.25 243.45 490.55 744.06 950.33
pi 0.57 0.04 0.50 0.54 0.57 0.60 0.63
pi_prior 0.67 0.23 0.22 0.50 0.71 0.86 0.97
num_echantillon 1250.50 721.74 125.95 625.75 1250.50 1875.25 2375.05

3 pêches

On reprend quasiment les mêmes instructions en ajoutant simplement quelques lignes pour modéliser la troisième pêche.

model2<-"model
{
#priors
   nuReal ~ dunif(0,500)
   nu_prior ~ dunif(0,500)
   nu <- round(nuReal)
   pi ~ dbeta(2,1)
   pi_prior ~ dbeta(2,1)
#Vraisemblance
   C1 ~ dbin(pi,nu)
   PasPris <- nu - C1
   C21 ~ dbin(pi,C1)
   C20 ~ dbin(pi,PasPris)
# troisieme peche
  C321 ~ dbin(pi,C21)
  UnMaisPasDeux<-C1-C21
  C301 ~ dbin(pi,UnMaisPasDeux)
  C320 ~ dbin(pi,C20)
  JamaisPris <- nu - C1-C20
  C300 ~ dbin(pi,JamaisPris)
}
"
C1=99; C20=56;  C21=62;C300=41;C321=43;C301=15;C320=38;
data_jags_bis<-list(C1=C1,C20=C20, C21=C21,C300=C300,C321=C321,C301=C301,C320=C320)
params <- c("pi", "pi_prior", "nu", "nu_prior")
min_nu <- C1 + C20 + 1+C300
inits <- function(){
  list(pi = rbeta(1, 1, 1), pi_prior = rbeta(1, 1, 1),
       nuReal = runif(1, min_nu, 500))
}
modelebis <- jags.model(file = textConnection(object = model2),
                     data = data_jags_bis, 
                     inits = list(inits(),inits(),inits()), 
                     n.adapt =  n_burn,
                     n.chains = n_chains, quiet = FALSE)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 4
##    Total graph size: 19
## 
## Initializing model
update(modelebis,n_burn)
echantillons_mcmc_bis <- coda.samples(model = modelebis, 
                                  variable.names = params,
                                  n.iter = 3*n_iter, 
                                  thin = n_thin)

summary(echantillons_mcmc_bis)
## 
## Iterations = 10010:85000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 7500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##              Mean        SD  Naive SE Time-series SE
## nu       216.6025   6.07090 0.0404727       0.041120
## nu_prior 250.3625 144.17838 0.9611892       0.949298
## pi         0.5460   0.02453 0.0001635       0.000166
## pi_prior   0.6691   0.23562 0.0015708       0.001571
## 
## 2. Quantiles for each variable:
## 
##              2.5%      25%      50%      75%    97.5%
## nu       206.0000 212.0000 216.0000 220.0000 230.0000
## nu_prior  12.4547 125.0699 251.5266 375.2688 486.6905
## pi         0.4977   0.5295   0.5463   0.5627   0.5935
## pi_prior   0.1567   0.5041   0.7120   0.8677   0.9874
echantillons_posterior_bis <- echantillons_mcmc_bis %>% 
  map_dfr(# On applique aux 3 éléments de la liste la même function
    function(tab_mcmc){ # tab mcmc
      resultat <- tab_mcmc %>% # On prend le tableau mcmc
        as.data.frame() %>% # On transforme en data.frame
        mutate(num_echantillon = 1:nrow(tab_mcmc)) # On créé un indice de numéro
    },
    .id = "chaine") %>% # On créée une colonne chaine qui stocke le numero de chaine
  as_tibble()

echantillons_post=rbind(echantillons_posterior %>% mutate(nbpech="Deux"),
                        echantillons_posterior_bis %>% mutate(nbpech="Trois")
                        )

####
library(cowplot) 
# Main plot
df<-echantillons_post %>% select(pi,nu,nbpech) %>% mutate(Peches= nbpech)
pmain<-ggpubr::ggscatter(data = df, x = "pi", y = "nu",
                         add = "reg.line",                         # Add regression line
                         color = "Peches", palette = "jco",           # Color by groups 
                         shape = "Peches",                            # Change point shape by groups 
                         fullrange = TRUE,                         # Extending the regression line
                         rug = FALSE# Add marginal rug
)
# Marginal densities along x axis
xdens <- axis_canvas(pmain, axis = "x")+
  geom_density(data = df, aes(x = pi, fill = Peches),
               alpha = 0.7, size = 0.2)+
  ggpubr::fill_palette("jco")
# Marginal densities along y axis
# Need to set coord_flip = TRUE, if you plan to use coord_flip()
ydens <- axis_canvas(pmain, axis = "y", coord_flip = TRUE)+
  geom_density(data = df, aes(x = nu, fill = Peches),
               alpha = 0.7, size = 0.1)+
  coord_flip()+
  ggpubr::fill_palette("jco")
p1 <- insert_xaxis_grob(pmain, xdens, grid::unit(.2, "null"), 
                        position = "top")
p2<- insert_yaxis_grob(p1, ydens, grid::unit(.2, "null"), 
                       position = "right")
ggdraw(p2)

Lois conditionnelles complètes de \(\pi\) et \(\nu\)

Dans cette annexe, nous supposons que \(\pi_{1} =\pi_{2}=\pi\). Par ailleurs, \(\pi \sim dbeta(a,b)\) et \(\nu\) suit une loi discrète de support fini.

La loi jointe des inconnues et des observables est donnée par: \[\begin{align*} \lbrack\nu,\pi,C_{1},C_{20},C_{21}] & =[C_{21}|C_{1} ,\pi][C_{20}|C_{1},\nu ,\pi][C_{1}|\nu,\pi][\nu][\pi]\\ & =\frac{\nu!\pi^{C_{1}+C_{2}}(1-\pi)^{2\nu-C_{1}-C_{2}}[\nu][\pi ]}{C_{21}!(C_{1}-C_{21})!C_{20}!(\nu-C_{1}-C_{20})!} \end{align*}\]

On en déduit immédiatement (en focalisant sur le rôle de ses arguments \(\pi\) puis \(\nu\)) les lois conditionnelles complètes de \(\pi\) et de \(\nu\): \[\begin{align*} \pi|\nu,C_{1},C_{20},C_{21} & \sim dbeta(a+C_{1}+C_{2}% ,b+2\nu-C_{1}-C_{2})\\ \lbrack \nu|\pi, C_{1},C_{20},C_{21}] & = \frac{(1-\pi)^{2\nu}\nu![\nu]}{(\nu -C_{1}-C_{20})!}\times Cste \end{align*}\]

L’écriture d’un échantillonneur de Gibbs et la comparaison de ses sorties par rapport à celles de Jags sont laissées au lecteur.