Combien de faces a un dé? Inférence bayésienne dans JAGS

Problème

Je voudrais faire une inférence sur un système analogue à mourir avec un nombre inconnu de côtés. Le dé est lancé plusieurs fois, après quoi je voudrais déduire une distribution de probabilité sur un paramètre correspondant au nombre de côtés du dé, θ.

Intuition

Si après 40 rouleaux vous aviez observé 10 rouges, 10 bleus, 10 verts et 10 jaunes, il semble que θ devrait culminer à 4, et les biais de roulement de chaque côté sont des distributions centrées sur 1/4.

θ a une borne inférieure triviale, étant le nombre de côtés différents observés dans les données.

La limite supérieure est encore inconnue. Il pourrait y avoir un cinquième côté qui aurait probablement un faible biais. Plus vous observez de données dépourvues de cinquième catégorie, plus la probabilité postérieure de θ = 4 est élevée.

Approche

J'ai utilisé JAGS pour des problèmes similaires (via R et rjags), ce qui semble approprié ici.

En ce qui concerne les données, disons obs <- c(10, 10, 10, 10)correspondre aux observations de l'exemple ci-dessus.

Je pense que les observations devraient être modélisées avec une distribution multinomiale obs ~ dmulti(p, n), où p ~ ddirch(alpha)et n <- length(obs).

θ est lié au nombre de catégories impliquées par alpha, alors comment puis-je modéliser alphapour englober différents nombres possibles de catégories?

Des alternatives?

Je suis assez nouveau dans les analyses bayésiennes, donc peut-être aboyer complètement le mauvais arbre, existe-t-il des modèles alternatifs qui pourraient fournir des informations différentes sur ce problème?

Merci beaucoup! David

Il s'agit d'un problème intéressant appelé «échantillonnage des espèces», qui a fait l'objet de beaucoup d'attention au fil des ans et englobe de nombreux autres problèmes d'estimation (comme le marquage-recapture). Il suffit de dire que JAGS ne vous aidera pas dans ce cas - JAGS ne peut pas gérer les chaînes de Markov avec une dimension variable sur plusieurs itérations. Il faut recourir à un schéma MCMC conçu pour de tels problèmes, comme le MCMC à saut réversible.

Voici une approche qui convient au modèle spécifique que vous décrivez, que j'ai rencontrée pour la première fois dans le travail de Jeff Miller ( arxived ).

Partie I (question d'origine)

Une hypothèse que je ferai est qu'une observation d'une catégorie donnée implique l'existence de catégories de rang inférieur. Autrement dit, l'observation d'un jet de dé sur le côté 9 implique l'existence des côtés 1-8. Il n'a pas a être de cette façon - les catégories pourrait être arbitraire - mais je suppose que oui dans mon exemple. Cela signifie que 0 valeurs sont observables, contrairement à d'autres problèmes d'estimation des espèces.

Disons que nous avons un échantillon multinomial,

$$Y = \{y_1, y_2, \dots, y_m, y_{m+1}, \dots, y_{n} \} \sim \mathcal{M}(\{p_1, p_2, \dots, p_m, p_{m+1}, \dots, p_n\})$$

où est la catégorie maximale observée, est le nombre (inconnu) de catégories, et tous égaux à 0. Le paramètre est fini, et nous avons besoin un avant pour cela. Tout préalable discret et approprié avec support sur fonctionnera; prenons par exemple un Poisson tronqué à zéro:n { y m + 1 , … , y n } n [ 1 , ∞ $m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

$$n \sim \mathcal{P}(\lambda), n > 0$$

Un a priori pratique pour les probabilités multinomiales est le Dirichlet,

$$P = \{ p_1, \dots, p_n \} \sim \mathcal{D}(\{ \alpha_1, \dots, \alpha_n \})$$

Et pour supposer .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Pour rendre le problème plus traitable, nous marginalisons les poids:

$$p(Y|\tilde{\alpha}, n) = \int_P p(Y|P, n)p(P|\tilde{\alpha}, n) dP$$

Ce qui, dans ce cas, conduit à la distribution multinomiale de Dirichlet bien étudiée . Le but est alors d'estimer le postérieur conditionnel,

$$p(n|Y, \tilde{\alpha}, \lambda) = \frac{ p(Y|n, \tilde{\alpha}) p(n|\lambda) }{ p(Y|\tilde{\alpha}, \lambda) }$$

Où je suppose explicitement que et sont des hyperparamètres fixes. Il est facile de voir que:$\tilde{\alpha}$$\lambda$

$$p(Y|\tilde{\alpha}, \lambda) = \sum_{n=1}^\infty p(Y|n, \tilde{\alpha}) p(n|\lambda)$$

Où où . Cette série infinie devrait converger assez rapidement (tant que la queue du prieur n'est pas trop lourde), et est donc facile à estimer. Pour le Poisson tronqué, il a la forme:n < $p(Y|n, \tilde{\alpha}) = 0$$n < m$

$$p(Y|\tilde{\alpha}, \lambda) = \frac{1}{(e^\lambda - 1)} \sum_{n=m}^\infty \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!}$$

Menant à:

$$p(n|Y,\tilde{\alpha}, \lambda) = \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!} \cdot \left(\sum_{j=m}^\infty \frac{\Gamma(j\tilde{\alpha})\prod_{i=1}^j \Gamma(y_i + \tilde{\alpha})}{\Gamma(j\tilde{\alpha} + \sum_{i=1}^j y_i) \Gamma(\tilde{\alpha})^j} \cdot \frac{\lambda^j}{j!}\right)^{-1}$$

Qui a un support sur . Il n'y a pas besoin de MCMC dans ce cas car la série infinie du dénominateur de la règle de Bayes peut être approximée sans trop d'effort.$[m, \infty)$

Voici un exemple bâclé dans R:

logPosteriorN <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        posterior <- lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { posterior <- posterior + (j-m)*lgamma(alpha) } 
        else if( j < m ) { posterior = -Inf }
        prior + posterior
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

## with even representation of sides
Y <- c(10, 10, 10, 10)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

## with uneven representation of sides
Y <- c(1, 2, 1, 0, 0, 2, 1, 0, 1)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

Votre intuition est correcte: un échantillonnage épars entre les catégories conduit à une plus grande incertitude sur le nombre total de catégories. Si vous souhaitez traiter comme un paramètre inconnu, vous devez utiliser MCMC et des mises à jour alternatives de et . n ˜ $\tilde{\alpha}$$n$$\tilde{\alpha}$

Bien sûr, il s'agit d'une approche d'estimation. Vous en trouverez facilement d'autres (de saveurs bayésiennes et non bayésiennes) avec un peu de recherche.

Partie II (réponse au commentaire)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ est un vecteur multinomial partiellement observé avec les probabilités correspondantes : $\Omega = \{\omega_1, \dots, \omega_m, \omega_{m+1}, \dots, \omega_n\}$

$$\mathrm{Pr}(Y|\Omega, n) = \frac{\Gamma(\sum_{i=1}^n y_i + 1)}{\prod_{i=1}^n \Gamma(y_i + 1) } \prod_{i=1}^n \omega_i^{y_i}$$

Où , et mais sinon les indices sont aléatoires. Comme précédemment, le problème est d'inférer le vrai nombre de catégories , et nous commençons par un a priori sur tel qu'un Poisson tronqué à zéro: $y \in \mathbb{N}$$y_1 \dots y_m > 0$$y_{m+1} \dots y_n = 0$$n$$n$

$$\mathrm{Pr}(n|\lambda) = \frac{\lambda^{n}}{(\exp\{\lambda\} - 1)n!},~n \in \mathbb{Z}^+$$

De même que précédemment, nous traitons les probabilités multinomiales comme Dirichlet distribué avec un hyperparamètre symétrique , c'est-à-dire pour un donné , $\Omega$$\tilde{\alpha}$$n$

$$\mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n\tilde{\alpha})}{\Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \omega_i^{\tilde{\alpha}-1}$$

L'intégration (marginalisation) sur le vecteur de probabilités donne le Dirichlet multinomial:

$$\mathrm{Pr}(Y|\tilde{\alpha}, n) = \int \mathrm{Pr}(Y|\Omega, n) \mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n \tilde{\alpha})} {\Gamma(\sum_{i=1}^n y_i + n \tilde{\alpha}) \Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})$$

Voici où nous nous écartons du modèle de la partie I ci-dessus. Dans la partie I, il y avait un ordre implicite aux catégories: par exemple, dans un dé à côtés, les catégories (côtés) ont un ordre implicite et l'observation de toute catégorie implique la existence de catégories inférieures . Dans la partie II, nous avons un vecteur aléatoire multinomial partiellement observé qui n'a pas d'ordre implicite. En d'autres termes, les données représentent une partition non ordonnée des points de données en catégories observées. Je dénoterai la partition non ordonnée qui résulte de augmentée de catégories non observées, comme .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

La probabilité de la partition non ordonnée conditionnée à un vrai nombre de catégories , peut être trouvée en considérant le nombre de permutations de catégories qui aboutissent à la même partition: $n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) = \frac{n!}{(n-m)!} \mathrm{Pr}(Y|\tilde{\alpha}, n)$$

Et cela peut être intégré sur pour donner: $n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda) = \sum_{j=m}^{\infty} \mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) \mathrm{Pr}(n|\lambda)$$

Utilisation de la règle de Bayes pour récupérer le postérieur:

$$\mathrm{Pr}(n|\mathcal{P}[Y], \tilde{\alpha}, \lambda) = \frac{\mathrm{Pr}(\mathcal{P}[Y]|n, \tilde{\alpha}) \mathrm{Pr}(n|\lambda)}{\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda)}$$

Branchez-vous simplement à partir des définitions ci-dessus. Encore une fois, le dénominateur est une série infinie qui convergera rapidement: dans ce modèle simple, il n'est pas nécessaire que MCMC donne une approximation adéquate.

En modifiant le code R de la partie I:

logPosteriorN_2 <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        likelihood <- lchoose(j, m) + lgamma(m + 1) + lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { likelihood <- likelihood + (j-m)*lgamma(alpha) } 
        else if( j < m ) { likelihood = -Inf }
        prior + likelihood
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

Y_1 <- rep(10, 15)
pos_1 <- logPosteriorN_2(50, Y_1, 6, 1)
plot(1:50, exp(pos_1))