Multinomial with small counts and many cells I have data distributed as multinomial, with $N= 1,000$ observations divided among $n=10,000$ cells, with unequal probabilities $p_i$, for $i=1,...,n$.
I would like to assess how likely it is that an observed $x$ arose from this distribution. I cannot use the Pearson's chi-squared test, as my expected cell counts are way too small. I don't think Fisher's exact test can be used in this setting either. How can I go about computing this p-value?

Some background, per the comments:
I have a dataset consisting of several thousand observations $x_i \in \mathbb{Z}_{\geq 0}^n$, where each $x_i$ is a sample from a multinomial distribution, possibly with different probabilities. My overall goal is to determine if all the $x_i$ derived from the same distribution. I thought to approach this by estimating $p$ as a weighted average of all the observations in the dataset and then determining the probability that $x_i$ could have been generated from Multinomial$(p, n_i)$. It is not clear to me how to accomplish this, and that was what inspired this question.

As @whuber noted, the rule of thumb that expected cells counts should >= 5 is overly conservative for a multinomial distribution. It is, however, appropriate for a hypergeometric distribution and my understanding is that this is how it was originally introduced (at least by Cochran). For more on this, see here.
In any case, I agree with the suggestion above that you can simulate the p-value by sampling from the contingency table under the null. This way you don't have to worry about whether or not the approximation to the Chi squared distribution is adequate.
The R chisq.test function includes an option to simulate p-values. It only offers hypergeometric sampling. For small tables, this can produce very different values than multinomial sampling, but if I'm not mistaken for a table of this size, there should be little difference.
However, if the the probabilities vary significantly, you might should be concerned about the Chi squared statistic. Because it weights by expected values, it can return anomalous results in some cases. An alternative that might work better under conditions of significant variance would be to simply use an unweighted Euclidean distance. Ward et al. suggest the following:
$$ RMS = \sqrt{N^{-1} \sum_{i} \sum_j (q_{ij} - p_{ij})^2} \tag{11} $$
See here.
