# Confidence Intervals for 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?

• "Or perhaps all these differences in real life applications might be negligible and they are all equally suitable?" I agree with this statement. I find a lot of statistics is actually a lot simpler than it first appears and many people like to over complicate the subject by decorating it with entirely unnecessarily methods and vocabulary. At the end of the day this is a subject suited for applications, so do whatever method is reliable and the simplest. Just think about this, how often is the answer "just take the average of the numbers", to a statistical problem? Commented Dec 30, 2022 at 3:06
• @ Nicolas Bourbaki: thank you so much for your reply! Commented Dec 30, 2022 at 5:14
• This article towardsdatascience.com/…. may expand your horizons and also answer some of your questions. It looks at five methods with sample size 100 over the range of true probabilities, with plots of performance, commentary, and R code! Commented Jan 3, 2023 at 3:19
• @ jbowman: thank you for your suggestion! Commented Jan 3, 2023 at 3:42
• Are you looking for a confidence interval, seperately for each number on the dice, or are you looking for a confidence region for the different numbers together, which have a correlation? In the first case the problem is just equivalent to finding confidence intervals for the parameter $p$ of a binomial distribution. Commented Jan 3, 2023 at 10:14

Are you willing to entertain a Bayesian approach? If so, you could specify a uniform prior on the five-dimensional surface $$\sum\limits_{x=1}^{6} p_x = 1$$, sample from the posterior distribution, and identify the $$100 (1-\alpha)$$% region with the highest posterior density.

Algorithm sketch

loop:

• Gibbs sample parameter vector $$\theta_i$$ from prior surface.
• Compute likelihood $$w_i \propto \mathcal{L}(\theta_i | X)$$.
• Accumulate sums
• weight $$\sum w_i$$
• weighted vectors $$\sum w_i \theta_i$$
• weighted squared vectors $$\sum w_i \theta_i \theta_i^\textrm{T}$$

output:

• multivariate posterior sample mean and variance

library(dplyr)

# 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)
x <- tibble(value=observations) %>% count(value) %>% rename(x=n) %>% pull(x)
print(x)
## [1]  8 18 32 22 11  9
theta = rep(1 / 6, 6) # starting point for Gibbs sampler

sum_wt <- 0
sum_wt_x <- rep(0,6)
sum_wt_xxt <- matrix(0,nrow=6,ncol=6)

for (i in 1:1e5) {

# Gibbs sampler - reallocate between two elements of vector
idx <- sample(1:6, 2)
w <- sum(theta[idx])
u <- runif(1)
theta[idx] <- c(u * w, (1 - u) * w)

# weight sampled theta proportional to likelihood (ignore constant)
wt <- exp((log(theta) %*% x))[1]
sum_wt <- sum_wt + wt
sum_wt_x <- sum_wt_x + wt*theta
sum_wt_xxt <- sum_wt_xxt + wt* (theta %*% t(theta))
}

mu <- sum_wt_x / sum_wt[1]
sigma <- (sum_wt_xxt/sum_wt[1] - mu %*% t(mu))

print(mu)
## [1] 0.08614837 0.18522926 0.31384024 0.21021317 0.11096186 0.09360710
print(sigma,digits=5)
##             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
## [1,]  6.5268e-04 -0.00011527 -0.00022615 -0.00015234 -9.4170e-05 -6.4751e-05
## [2,] -1.1527e-04  0.00155075 -0.00063112 -0.00044187 -2.0724e-04 -1.5525e-04
## [3,] -2.2615e-04 -0.00063112  0.00201363 -0.00062085 -2.1323e-04 -3.2228e-04
## [4,] -1.5234e-04 -0.00044187 -0.00062085  0.00160344 -2.2775e-04 -1.6063e-04
## [5,] -9.4170e-05 -0.00020724 -0.00021323 -0.00022775  8.0827e-04 -6.5879e-05
## [6,] -6.4751e-05 -0.00015525 -0.00032228 -0.00016063 -6.5879e-05  7.6879e-04


mu is the posterior mean, sigma is the posterior covariance matrix. Note off diagonal entries of sigma are negative. This is expected. A positive error in one element of mu is offset by negative errors in the other elements of mu.

• The posterior is a Dirichlet distribution for which the mean and covariance matrix can be computed with simple closed formulas. A problem is that the mean and covariance matrix do not yet give you a highest density region and neither a confidence interval. Commented Jan 9, 2023 at 6:48
• One thing that one might do with sampling is compute the distribution of the density (using the known posterior) to find a cut-off density for the highest density region. Then use that to compare with the density at the point for a given hypothesis. (I guess that it will be difficult to describe the region parametrically, it will be a region within the surface defined by $\sum \alpha_i \log(p) = \text{const}$ and $\sum p = 1$) Commented Jan 9, 2023 at 6:59
• @SextusEmpiricus thanks for pointing out the existence of the simple closed formulas. As to your second point, one could accumulate the posterior density above and below the likelihood of the comparison point, and thereby determine if the comparison point lies in the desired HPD region. Commented Jan 9, 2023 at 13:18
• Skilling, John. "Nested sampling for general Bayesian computation." Bayesian Analysis 1.4 (2006): 833-859. Commented Jan 9, 2023 at 13:25
• In addition, for a particular comparison point one might also compute a hypothesis test. The biggest problem is expressing the entire region. Commented Jan 9, 2023 at 13:38

When you speak about confidence intervals rather than confidence regions then the problem is equivalent to confidence intervals for a binomial distribution. There are several sources that describe this situation:

In the case that you are looking for a confidence region, then the situation is different (more complicated). In the comments there is a reference about the Sison–Glaz method described in Sison and Glaz's "Simultaneous Confidence Intervals and Sample Size Determination for Multinomial Proportions". And a software implementation is in the R-package MultinomCI.