¿Cuántos lados tiene un dado? Inferencia bayesiana en JAGS

Problema

Me gustaría hacer alguna inferencia en un sistema análogo a morir con un número desconocido de lados. El dado se tira varias veces, después de lo cual me gustaría inferir una distribución de probabilidad sobre un parámetro correspondiente al número de lados que tiene el dado, θ.

Intuición

Si después de 40 rollos ha observado 10 rojos, 10 azules, 10 verdes y 10 amarillos, parece que peak debería alcanzar un máximo de 4, y los sesgos de rodar a cada lado son distribuciones centradas en 1/4.

θ tiene un límite inferior trivial, que es el número de lados diferentes observados en los datos.

El límite superior es aún desconocido. Podría haber un quinto lado que probablemente tenga un sesgo bajo. Cuantos más datos observe que carecen de una quinta categoría, mayor será la probabilidad posterior de θ = 4.

Acercarse

He usado JAGS para problemas similares (a través de R y rjags) que parece apropiado aquí.

Con respecto a los datos, digamos obs <- c(10, 10, 10, 10)que corresponden a las observaciones en el ejemplo anterior.

Creo que las observaciones deberían modelarse con una distribución multinomial obs ~ dmulti(p, n), donde p ~ ddirch(alpha)y n <- length(obs).

θ está vinculado al número de categorías implicadas alpha, entonces, ¿cómo puedo modelar alphapara abarcar diferentes números posibles de categorías?

¿Alternativas?

Soy bastante nuevo en los análisis bayesianos, por lo que podría estar ladrando por completo el árbol equivocado, ¿hay modelos alternativos que puedan proporcionar diferentes puntos de vista sobre este problema?

¡Muchas gracias! David

Este es un problema interesante denominado 'muestreo de especies', que ha recibido mucha atención a lo largo de los años, y abarca muchos otros problemas de estimación (como la recuperación de marcas). Baste decir que JAGS no lo ayudará en este caso: JAGS no puede manejar cadenas de Markov con una dimensión variable a través de iteraciones. Uno debe recurrir a un esquema MCMC diseñado para tales problemas, como el MCMC de salto reversible.

Aquí hay un enfoque que es adecuado para el modelo específico que está describiendo, que encontré por primera vez en el trabajo de Jeff Miller ( presentado ).

Parte I (Pregunta original)

Una suposición que haré es que una observación de una categoría dada implica la existencia de categorías de menor rango. Es decir, observar un lanzamiento de dado en el lado 9 implica la existencia de los lados 1-8. No tiene por qué ser así, las categorías pueden ser arbitrarias, pero supongo que sí en mi ejemplo. Esto significa que los valores 0 son observables, en contraste con otros problemas de estimación de especies.

Digamos que tenemos una muestra 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\})$$

donde es la categoría máxima observada, es el número (desconocido) de categorías, y todas iguales a 0. El parámetro es finito y necesitamos Un previo para ello. Cualquier previo discreto y adecuado con soporte en funcionará; tomemos, por ejemplo, un Poisson truncado a cero:$m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

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

Un previo conveniente para las probabilidades multinomiales es el Dirichlet,

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

Y para asumir de manera simple .$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Para hacer el problema más manejable, marginamos los pesos:

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

Lo que en este caso conduce a la bien estudiada distribución multinomial de Dirichlet . El objetivo es entonces estimar el posterior condicional,

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

Donde supongo explícitamente que y son hiperparámetros fijos. Es fácil ver eso:$\tilde{\alpha}$$\lambda$

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

Donde donde . Esta serie infinita debería converger bastante rápido (siempre que la cola de la anterior no sea demasiado pesada), por lo que es fácil de aproximar. Para el Poisson truncado, tiene la forma:$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!}$$

Llevando 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 tiene soporte en . No hay necesidad de MCMC en este caso, ya que las series infinitas en el denominador de la regla de Bayes se pueden aproximar sin demasiado esfuerzo.$[m, \infty)$

Aquí hay un ejemplo descuidado en 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")

Su intuición es correcta: el escaso muestreo entre categorías genera una mayor incertidumbre sobre el número total de categorías. Si desea tratar como un parámetro desconocido, necesitaría usar MCMC y actualizaciones alternativas de y .$\tilde{\alpha}$$n$$\tilde{\alpha}$

Por supuesto, este es un enfoque para la estimación. Encontrarás fácilmente otros (de sabores bayesianos y no bayesianos) con un poco de búsqueda.

Parte II (respuesta al comentario)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ es un vector multinomial parcialmente observado con probabilidades correspondientes : $\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}$$

Donde , e pero, de lo contrario, los índices son arbitrarios. Como antes, el problema es inferir el verdadero número de categorías , y comenzamos con un anterior en como un Poisson truncado a cero: $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}^+$$

También como antes, tratamos las probabilidades multinomiales como Dirichlet distribuido con un hiperparámetro simétrico , es decir, para un 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}$$

Integrar (marginar) sobre el vector de probabilidades da el 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})$$

Aquí es donde divergemos del modelo en la Parte I anterior. En la Parte I, había un orden implícito para las categorías: por ejemplo, en un dado de lados, las categorías (lados) tienen un orden implícito y la observación de cualquier categoría implica el existencia de categorías menores . En la Parte II, tenemos un vector aleatorio multinomial parcialmente observado que no tiene un orden implícito. En otras palabras, los datos representan una partición desordenada de los puntos de datos en categorías observadas. Denotaré la partición desordenada que resulta de aumentada por categorías no observadas, como .$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

La probabilidad de que la partición desordenada esté condicionada a un número verdadero de categorías , se puede encontrar considerando el número de permutaciones de categorías que resultan en la misma partición: $n$

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

Y esto se puede integrar sobre para dar: $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 la regla de Bayes para recuperar la parte 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)}$$

Simplemente conéctese desde las definiciones anteriores. Nuevamente, el denominador es una serie infinita que convergerá rápidamente: en este modelo simple, no es necesario que MCMC dé una aproximación adecuada.

Al modificar el código R de la 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))