Calculating the "Spreads" for Different Outcomes in Dice Rolls?

Suppose I roll a 6-sided die 100 times and observe the following data - let's say that I don't know the probability of getting any specific number (but I am assured that each "trial" is independent from the previous "trial").

Below, here is some R code to simulate this experiment:

# Set the probabilities for each number (pretend this is unknown in real life)
probs <- c(0.1, 0.2, 0.3, 0.2, 0.1, 0.1)

# Generate 100 random observations
observations <- sample(1:6, size = 100, replace = TRUE, prob = probs)

# Print the observations
print(observations)

[1] 2 4 2 2 4 6 2 2 6 6 3 4 6 4 2 1 3 6 3 1 2 5 3 6 4 6 1 3 4 2 6 2 4 1 3 3 3 5 2 5 2 3 5 1 4 6 1 6 4 2
[51] 2 3 2 3 3 5 6 5 4 3 2 3 2 1 2 3 2 2 5 3 2 1 1 1 3 3 2 4 4 3 1 4 4 6 3 3 5 5 2 2 1 3 2 1 6 3 4 3 3 3


As we know, the above experiment corresponds to the Multinomial Probability Distribution Function (https://en.wikipedia.org/wiki/Multinomial_distribution):

$$P(X_1 = x_1, X_2 = x_2, \dots, X_k = x_k) = \frac{n!}{x_1!x_2!\dots x_k!}p_1^{x_1}p_2^{x_2} \dots p_k^{x_k}$$

Using Maximum Likelihood Estimation (https://en.wikipedia.org/wiki/Maximum_likelihood_estimation MLE), the estimate for the probability for getting any number on this die is given by (e.g. what is the probability that this above die gives you a "3"?):

$$\hat{p}_{i,\text{MLE}} = \frac{x_i}{n}$$

Next, the Variance for each of these parameters can be written as follows :

$$\text{Var}(\hat{p}_{i,\text{MLE}}) = \frac{p_i(1 - p_i)}{n}$$

From here, I am interested in estimating the "spreads" of these probabilities - for example, there might be a 0.2 probability of getting a "6" - but we can then "bound" this estimate and say there is a 0.2 ± 0.05 probability of rolling a 6. Effectively, this "bounding" corresponds to a Confidence Interval (https://en.wikipedia.org/wiki/Confidence_interval).

Recently, I learned that when writing Confidence Intervals for "proportions and probabilities", we might not be able to use the "classic" notion of the Confidence Interval (i.e. parameter ± z-alpha/2*sqrt(var(parameter))), because this could result in these bounds going over "1" and below "0", thus violating the fundamental definitions of probability.

Doing some reading online, I found different methods that might be applicable for writing the Confidence Intervals for the parameters of a Multinomial Distribution.

• Bootstrapping (https://en.wikipedia.org/wiki/Bootstrapping_(statistics)): By virtue of the Large Law of Large Numbers (https://en.wikipedia.org/wiki/Law_of_large_numbers), Bootstrapping works by repeatedly resampling your observed data and using this MLE formulas to calculate the parameters of interest on each of these re-samples. Then, you would sort the parameter in estimates in ascending order and take the estimates corresponding to the 5th and 95th percentiles. These estimates from the 5th and 95th percentiles would now correspond to the desired Confidence Interval. As I understand, this is an approximate method, but I have heard that the Law of Large Numbers argues that for an infinite sized population and an infinite number of resamples, the bootstrap estimates will converge to the actual values. It is important to note that in this case, the "Sequential Bootstrap" approach needs to be used such that the chronological order of the observed data is not interrupted.

• Delta Method (https://en.wikipedia.org/wiki/Delta_method): The Delta Method uses a Taylor Approximation (https://en.wikipedia.org/wiki/Taylor%27s_theorem) for the function of interest (i.e. MLE variance estimate). Even though this is also said to be an approximate method (i.e. the Delta Method relies on the Taylor APPROXIMATION), there supposedly exists mathematical theory (e.g. https://en.wikipedia.org/wiki/Continuous_mapping_theorem) which can demonstrate that estimates from the Delta Method "converge in probability" to the actual values. This being said, I am not sure how the Delta method can directly be used to calculate Confidence Intervals.

• Finally, very recently I learned about the Wilson Interval (https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval), which is said to be more suitable for writing Confidence Intervals in the case of proportions and probabilities. In the case of the Multinomial Probability Distribution, I think the Wilson Interval for 95% Confidence Intervals on parameter estimates can be written as follows:

$$\left( \hat{\theta} - \frac{z_{\alpha/2} \sqrt{\hat{\theta}(1-\hat{\theta})/n}}{1+z_{\alpha/2}^2/n}, \hat{\theta} + \frac{z_{\alpha/2} \sqrt{\hat{\theta}(1-\hat{\theta})/n}}{1+z_{\alpha/2}^2/n} \right)$$

However, I am still learning about the details of this.

This brings me to my question: What are the advantages and disadvantages of using any of these approaches for calculating the Confidence Interval for parameter estimates in the Multinomial Distribution?

It seems like many of these methods are approximations - but I am willing to guess that perhaps some of these approximate methods might have better properties than others. As an example:

• Perhaps some of these methods might take longer to calculate in terms of computational power for more complex functions and larger sample sizes?
• Perhaps some of these methods might be less suitable for smaller sample sizes?
• Perhaps some of these methods are known to "chronically" overestimate or underestimate the confidence intervals?
• Perhaps some of these methods are simply "weaker" - i.e. the guarantee of the true parameter estimate lying between predicted ranges is not "as strong a guarantee"?

In any case, I would be interested in hearing about opinions on this matter - and in general, learning about which approaches might be generally more suitable for evaluating the Confidence Intervals on parameter estimates for the Multinomial Distribution.

Thanks!

Note: Or perhaps all these differences in real life applications might be negligible and they are all equally suitable?

• Your binomial proportion confidence interval Wikipedia link gives various different ways of finding a confidence interval for the probability of each face. The minima will add up to less than $1$ and the maxima to greater than $1$, but that is to be expected, as the probabilities are not independent given that they are constrained to sum to $1$ Commented Dec 18, 2022 at 0:24
• @ Henry: Thank you for pointing this out! Commented Dec 18, 2022 at 18:44
• Perhaps this would be better suited for stats.stackexchange.com? (The question is too old to propose migrating it there.) Commented Feb 19, 2023 at 7:19

You must be careful here: if you are interested in inference on a particular value of the die roll (e.g., you only care about interval estimates for the proportion of threes you observe), then the binomial proportion confidence intervals are appropriate. However, if you are interested in simultaneous inference on more than one multinomial parameter, then strictly speaking, the pairwise negatively correlated outcomes imply it is inappropriate to compute the individual binomial confidence intervals (using whichever method you please) because doing so implicitly assumes that they are independent.

Sison and Glaz (1995) describe a method to compute the simultaneous confidence intervals for multinomial proportions. May and Johnson (2000) implements this method in SAS/IML, and there is an implementation described in the R documentation--this is probably what you will want to use.

It is important to understand that the consequences of using a binomial model when simultaneous inference is desired, may include the possibility that your resulting intervals will not have the intended coverage probability, even if you choose to use the Clopper-Pearson interval.

The short answer is that you should probably use one of the methods on the Wikipedia page about binomial confidence intervals. While the overall distribution you're working with is multinomial, the individual $$p_i$$ you're trying to estimate can be treated as binomial; for each roll, you either see the face $$i$$, or you don't.

The paper by Brown, Cai and DasGupta has a good comparison of the different confidence intervals available; three sensible choices are the Wilson, Agresti-Couli and Clopper-Pearson methods.

These methods differ in how conservative they are. The Clopper-Pearson (exact) method guarantees that the probability of the true parameter lying within your confidence interval is at least 95%, something the other methods don't. However, the price for this is that the interval is sometimes larger than necessary; the actual probability of containing the truth can be higher than 95%.

The Wilson method is often recommended, and has a coverage that is usually close to 95%, while being smaller than Clopper-Pearson. However, when $$p$$ is close to $$0$$ or $$1$$, the Wilson interval may be inaccurate. Alternatively, the Agresti-Couli method lies between the other two methods. It's a little conservative, but less than Clopper-Pearson, and doesn't have issues with $$p$$ near $$0$$ or $$1$$; it's also quite easy to calculate.

The delta method is not so much a confidence interval, as a tool a mathematician might use to develop one. In this case, there are lots of confidence intervals already available; you don't need to use the delta method yourself.

Bootstrap methods are a more computationally-intensive approach, that can be useful when the underlying distribution is uncertain, or otherwise hard to work with. In the specific case of a binomial distribution, neither of those are true; you'd just be spending a lot of CPU time to get a worse approximation.

Edit: if you want a confidence set containing all the $$p_i$$ simultaneously, not just a set for each $$p_i$$ separately, then you should use a different approach. Standard techniques for this include the Pearson chi-squared test, G-test, and exact multinomial test, though the sets they provide are not a product of intervals. For a method that does, and discussion of why the separate intervals can sometimes be misleading, see heropup's answer.