# How to test the statistical hypothesis that data was generated from a multinomial distribution?

First, a quick definition of the multinomial distribution: Suppose that we were to roll a die with $k$ sides, and the probability that side $i$ comes up on top is $p_i$. If we run $n$ independent trials of rolling the die, and count the number of times $n_i$ that side $i$ comes up on top, then the random vector $(n_1, n_2, \ldots, n_k)$ is generated according to the multinomial distribution.

Now suppose that we have observed counts $(n_1, n_2, \ldots, n_k)$. Is there a way to statistically test an hypothesis of the assumption that the data was generated from the multinomial distribution with parameters $(p_1, p_2, \ldots, p_k)$?

This looks like a goodness-of-fit application to me. You should use the chi-square test. So if you have 3 sides and your probability vector is (0.1, 0.2, 0.7), and you have 100 trials. You would expect the outcomes be (10, 20, 70). Use the theoretical counts to compare with your observed counts in the chi-square test. Yes, it's called a Pearson's $\chi^2$ test. The expected frequencies are $E_i = N \times p_i$ (where $N$ is the total sample size) and the observed frequencies are what you called $n_1, ..., n_k$.

The existing answers point out Pearson's chi-squared test. Though this is usually a fine solution, it does rely on an approximation. Pearson derived that the chi-squared statistic approximately follows a chi-squared distribution, which he used to calculate p values. He did this derivation by assuming that the multinomial random vector is actually multivariate normal. So this works well when all the counts are pretty high, but not necessarily when they are low.

I want to point out four alternatives that do not rely on this approximation.

### 1. Enumeration

One is an exact multinomial test, with the p value calculated by enumeration. I mean, you calculate a badness-of-fit statistic from your data, either the chi-squared statistic or the likelihood ratio statistic. Then, you actually go through every possible multinomial random vector, check whether the test statistic for that vector exceeds the observed one, and add up the probabilities. The total probability of equaling or exceeding the observed value is your p-value. This method is implemented in the R package XNomial, as xmulti().

### 2. Monte Carlo simulation

Choose a measure of badness-of-fit to use as a test statistic. Again, probably either chi-squared or likelihood ratio statistic. Simulate samples from your theoretical distribution, and count what fraction of the time the statistic equals or exceeds the observed one. That fraction is your p-value. This method is also implemented in XNomial, as xmonte(). It's also really easy to code yourself, using rmultinom in R or the equivalent in another language.

### 3. Plug the highest observed value into the multinomial "CDF"

This only really works for a special case, where the null distribution is kind of "even", and you're testing whether the data is more "uneven" than it should be. For example, the astronomer Rudolph Wolf did an experiment where he rolled two dice 20,000 times. One die came up on side 6 most often--3932 times. The other came up on side 2 most often--3631 times. Just this information is sufficient to reject the null hypothesis that the dice are fair. Your test statistic is the number of observations of the most frequently observed side. Focusing on the second die, the probability that this number is less than observed is this, which looks kind of like a CDF:

$$P(X_1 \le 3630, X_2 \le 3630, X_3 \le 3630, X_4 \le 3630, X_5 \le 3630, X_6 \le 3630)$$

The probability that all sides were observed 3630 times or fewer, see? $X_i$ here means the number of times $i$ was observed. So this is the probability that the maximum of the $X$'s is 3630 or less. So, our p value--the probability that the maximum is 3631 or more--is one minus this probability.

This probability can be calculated in R using my package pmultinom. In fact it's the example code in the documentation for the "pmultinom" function. There's also another package called pmultinom. Theoretically, the other package is using an approximation whereas mine is an exact calculation, but I haven't gotten a chance to compare and see if the difference is large enough to matter.

So, this is a special case, and it can also be less powerful since you're only using the information from one possible outcome, but sometimes it's enough.

## 4. Approximate likelihood ratio test

And I guess there's a fourth possibility, which is to use a likelihood ratio test assuming 2 * log(likelihood ratio) has a chi-squared distribution. (where likelihood ratio is probability of the data given that the true probabilities are the observed frequencies, divided by the probability of the data given your hypothesized ones. This is upside-down relative to the definition on Wikipedia, but it makes more sense to me in this context, since it's high when the fit is bad, just like the chi-squared statistic.) This is also approximate, but maybe it works in different situations. Surprisingly I can't find an R function which does this. It should be very easy to code yourself. I'm not totally sure how many degrees of freedom to use for the chi-squared distribution though.

The application of the goodness of fit tests is the wrong solutions because they reject the hypothesis, which you like to confirm. The right technique the equivalence testing. Please, look at my papers https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2907258 and https://www.researchgate.net/publication/312481284_Testing_equivalence_of_multinomial_distributions. The implementation is also available on github. Do not hesitate to ask me if any questions occur.

• It is recommended that more context be provided for external links to make the answer more self contained, in order that the answer is still useful if a link breaks. See stats.stackexchange.com/help/how-to-answer – ReneBt Apr 5 at 10:37