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.
