# How many sides does a die have? Bayesian inference in JAGS

### Problem

I would like to do some inference on a system analogous to die with an unknown number of sides. The die is rolled several times, after which I would like to infer a probability distribution over a parameter corresponding to the number of sides the die has, θ.

### Intuition

If after 40 rolls you had observed 10 red's, 10 blues's, 10 greens's and 10 yellows's, it seems that θ should peak at 4, and the biases of rolling each side are distributions centred on 1/4.

θ has a trivial lower bound, being the number of different sides observed in the data.

The upper bound is still unknown. There could be a fifth side which would probably have a low bias. The more data you observe lacking a fifth category, the higher the posterior probability of θ = 4.

### Approach

I've used JAGS for similar problems (via R and rjags) which seems appropriate here.

With respect to the data, lets say obs <- c(10, 10, 10, 10) to correspond to the observations in the example above.

I think the observations should be modelled with a multinomial distribution obs ~ dmulti(p, n), where p ~ ddirch(alpha) and n <- length(obs).

θ is linked to the number of categories implied by alpha, so how can I model alpha to encompass different possible numbers of categories?

### Alternatives?

I'm pretty new to Bayesian analyses so might be barking up the wrong tree entirely, are there alternative models which might provide different insights on this problem?

Many thanks! David

## migrated from stackoverflow.comSep 1 '15 at 17:15

This question came from our site for professional and enthusiast programmers.

## 1 Answer

This an interesting problem termed 'species-sampling', that has received a lot of attention over the years, and encompasses many other estimation problems (such as mark-recapture). Suffice it to say, JAGS will not help you in this case--JAGS cannot handle Markov chains with a variable dimension across iterations. One must recourse to an MCMC scheme designed for such problems, such as reversible jump MCMC.

Here's one approach that's suitable for the specific model you are describing, which I first encountered in the work of Jeff Miller (arxived).

Part I (Original question)

One assumption I'll make is that an observation a given category implies the existence of categories of a lesser rank. That is, observing a die roll on side 9 implies the existence of sides 1-8. It doesn't have to be this way--the categories could be arbitrary--but I'll assume so in my example. This means that 0 values are observable, in contrast to other species-estimation problems.

Let's say we have a multinomial sample,

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

where $m$ is the maximum category observed, $n$ is the (unknown) number of categories, and all $\{y_{m+1},\dots,y_{n}\}$ equal 0. The parameter $n$ is finite, and we need a prior for it. Any discrete, proper prior with support on $[1, \infty)$ will work; take for example a zero-truncated Poisson:

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

A convenient prior for the multinomial probabilities is the Dirichlet,

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

And for simplicitly assume $\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$.

To make the problem more tractable, we marginalize out the weights:

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

Which in this case leads the well-studied Dirichlet-multinomial distribution. The goal is then to estimate the conditional posterior,

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

Where I'm explicitly assuming $\tilde{\alpha}$ and $\lambda$ are fixed hyperparameters. It is easy to see that:

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

Where $p(Y|n, \tilde{\alpha}) = 0$ where $n < m$. This infinite series should converge rather quickly (as long as the tail on the prior isn't too heavy), and so is easy to approximate. For the truncated Poisson, it has the form:

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

Leading to:

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

Which has support on $[m, \infty)$. There isn't a need for MCMC in this case as the infinite series in the denominator of Bayes' rule can be approximated without too much effort.

Here's a sloppy example in 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")


Your intuition is correct: sparse sampling across categories leads to greater uncertainty about the total number of categories. If you wished to treat $\tilde{\alpha}$ as an unknown parameter, you would need to use MCMC and alternate updates of $n$ and $\tilde{\alpha}$.

Of course, this is one approach to estimation. You will readily find others (of Bayesian and non-Bayesian flavors) with a little searching.

Part II (Answer to comment)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$ is a partially-observed multinomial vector with corresponding probabilities $\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}$$

Where $y \in \mathbb{N}$, $y_1 \dots y_m > 0$ and $y_{m+1} \dots y_n = 0$ but otherwise the indices are abitrary. As before, the problem is to infer the true number of categories $n$, and we start with a prior on $n$ such as a zero-truncated Poisson: $$\mathrm{Pr}(n|\lambda) = \frac{\lambda^{n}}{(\exp\{\lambda\} - 1)n!},~n \in \mathbb{Z}^+$$

Also as before, we treat the multinomial probabilities $\Omega$ as Dirichlet distributed with a symmetric hyperparameter $\tilde{\alpha}$, i.e. for a given $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}$$

Integrating (marginalizing) over the vector of probabilities gives the multinomial Dirichlet: $$\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})$$

Here is where we diverge from the model in Part I above. In Part I, there was an implicit ordering to categories: for example, in an $n$-sided die, the categories (sides) have an implicit ordering and the observation of any category $i \in \{1 \dots n\}$ implies the existence of lesser categories $j < i$. In Part II, we have a partially-observed multinomial random vector which has no implicit ordering. In other words, the data represent an unordered partition of the data points into $m \leq n$ observed categories. I'll denote the unordered partition that results from $Y$ augmented by $n-m$ unobserved categories, as $\mathcal{P}[Y]$.

The probability of the unordered partition conditional on a true number of categories $n$, can be found by considering the number of permutations of categories which result in the same partition: $$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) = \frac{n!}{(n-m)!} \mathrm{Pr}(Y|\tilde{\alpha}, n)$$

And this can be integrated over $n$ to give: $$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda) = \sum_{j=m}^{\infty} \mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) \mathrm{Pr}(n|\lambda)$$

Using Bayes' rule to retrieve the 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)}$$

Just plug in from the definitions above. Again, the denominator is an infinite series that will converge quickly: in this simple model, there is no need for MCMC to give an adequate approximation.

By modifying the R code from Part 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))

• Many thanks for your very complete answer. (Sorry for my very slow response). I've come back to this type of question and am still working my way through the maths. In my system the categories are not ordinal, so the assumption that an observation a given category implies the existence of categories of a lesser rank is invalid. – davipatti Nov 24 '15 at 16:28
• @davipatti Answered in second part. – Nate Pope Jan 10 '16 at 18:00