Quantos lados tem um dado? Inferência bayesiana no JAGS

Problema

Eu gostaria de fazer alguma inferência em um sistema análogo para morrer com um número desconhecido de lados. O dado é rolado várias vezes, após o qual eu gostaria de inferir uma distribuição de probabilidade sobre um parâmetro correspondente ao número de lados do dado, θ.

Intuição

Se, após 40 jogadas, você tivesse observado 10 vermelhos, 10 azuis, 10 verdes e 10 amarelos, parece que θ deveria ter um pico em 4, e os vieses de rolar cada lado são distribuições centradas em 1/4.

θ possui um limite inferior trivial, sendo o número de lados diferentes observado nos dados.

O limite superior ainda é desconhecido. Poderia haver um quinto lado que provavelmente teria um viés baixo. Quanto mais dados você observar sem uma quinta categoria, maior a probabilidade posterior de θ = 4.

Aproximação

Eu usei o JAGS para problemas semelhantes (via R e rjags), o que parece apropriado aqui.

Com relação aos dados, digamos obs <- c(10, 10, 10, 10)que correspondam às observações no exemplo acima.

Eu acho que as observações devem ser modeladas com uma distribuição multinomial obs ~ dmulti(p, n), onde p ~ ddirch(alpha)e n <- length(obs).

θ está vinculado ao número de categorias implícitas por alpha, então como posso modelar alphapara abranger diferentes números possíveis de categorias?

Alternativas?

Eu sou bastante novo nas análises bayesianas; portanto, pode estar latindo completamente a árvore errada; existem modelos alternativos que podem fornecer insights diferentes sobre esse problema?

Muito Obrigado! David

Esse é um problema interessante denominado 'amostragem de espécies', que recebeu muita atenção ao longo dos anos e abrange muitos outros problemas de estimativa (como a recuperação de marcas). Basta dizer que o JAGS não ajudará nesse caso - o JAGS não pode manipular cadeias de Markov com uma dimensão variável nas iterações. É preciso recorrer a um esquema MCMC projetado para problemas como o MCMC de salto reversível.

Aqui está uma abordagem adequada ao modelo específico que você está descrevendo, que eu encontrei pela primeira vez no trabalho de Jeff Miller ( arxived ).

Parte I (questão original)

Uma suposição que farei é que uma observação de uma determinada categoria implica a existência de categorias de menor classificação. Ou seja, observar um rolo de dado no lado 9 implica na existência dos lados 1-8. Não precisa ser assim - as categorias podem ser arbitrárias -, mas eu assumo isso no meu exemplo. Isso significa que 0 valores são observáveis, em contraste com outros problemas de estimativa de espécies.

Digamos que temos uma amostra 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\})$$

onde é a categoria máxima observada, é o número (desconhecido) de categorias e todos iguais a 0. O parâmetro é finito e precisamos um prior para isso. Qualquer prévia discreta e adequada com suporte em funcionará; Tomemos, por exemplo, um Poisson truncado com zero: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$$

Um prior conveniente para as probabilidades multinomiais é o Dirichlet,

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

E, por simplicidade, assuma .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Para tornar o problema mais tratável, marginalizamos os pesos:

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

Que neste caso lidera a bem estudada distribuição multinomial de Dirichlet . O objetivo é estimar o posterior condicional,

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

Onde eu estou assumindo explicitamente que e são hiperparâmetros fixos. É fácil ver que:$\tilde{\alpha}$$\lambda$

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

Onde onde . Essa série infinita deve convergir rapidamente (desde que a cauda do prior não seja muito pesada) e, portanto, fácil de aproximar. Para o Poisson truncado, ele tem a forma: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!}$$

Levando a:

$$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}$$

Que tem suporte em . Não há necessidade de MCMC nesse caso, pois as séries infinitas no denominador da regra de Bayes podem ser aproximadas sem muito esforço.$[m, \infty)$

Aqui está um exemplo desleixado em 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")

Sua intuição está correta: a amostragem esparsa entre categorias leva a uma maior incerteza sobre o número total de categorias. Se você deseja tratar como um parâmetro desconhecido, precisará usar o MCMC e atualizações alternativas de e .$\tilde{\alpha}$$n$$\tilde{\alpha}$

Obviamente, essa é uma abordagem para a estimativa. Você encontrará facilmente outros (de sabores bayesianos e não bayesianos) com um pouco de pesquisa.

Parte II (Resposta ao comentário)

Ω = { ω 1 , … , ω m$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ é um vetor multinomial parcialmente observado com probabilidades correspondentes : $\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}$$

Onde , e mas, caso contrário, os índices são abitrários. Como antes, o problema é inferir o número real de categorias , e começamos com um prior em como um Poisson truncado com zero: y 1 … y m$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}^+$$

Também como antes, tratamos as probabilidades multinomiais como Dirichlet distribuído com um hiperparâmetro simétrico , ou seja, para um dado , $\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}$$

A integração (marginalização) sobre o vetor de probabilidades fornece ao 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})$$

Aqui é onde divergimos do modelo na Parte I acima. Na Parte I, havia uma ordem implícita para categorias: por exemplo, em um dado de lados, as categorias (lados) têm uma ordem implícita e a observação de qualquer categoria implica a existência de categorias menores . Na Parte II, temos um vetor aleatório multinomial parcialmente observado, que não possui ordem implícita. Em outras palavras, os dados representam uma partição não ordenada dos pontos de dados em categorias observadas. Denotarei a partição não ordenada que resulta de aumentada por categorias não observadas, como .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

A probabilidade da partição não ordenada condicional a um número real de categorias pode ser encontrada considerando o número de permutações de categorias que resultam na mesma partição: $n$

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

E isso pode ser integrado sobre para fornecer: P r ( P [ Y ] $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)$$

Usando a regra de Bayes para recuperar o posterior:

$$\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)}$$

Basta conectar a partir das definições acima. Novamente, o denominador é uma série infinita que convergirá rapidamente: nesse modelo simples, não há necessidade de o MCMC fornecer uma aproximação adequada.

Modificando o código R da parte 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))