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))
