Confidence Intervals for Dice Rolls?

Question

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?

Answer 1

See also link in @jbowman comment, method 5.

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.

Answer 2

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:

• Wikipedia has a large article about the Binomial proportion confidence interval
• In the comments jbowman referred to a blogpost with a similar content as the article from Wikipedia.

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.