Quantifying the confidence that the most sampled outcome is the most probable outcome

Question

Let us say we have a cheater's six-sided die, which we can assume to be unfair with an unknown probability distribution. We want to know the most likely roll with this die, and so we roll it $N \gg 6$ times and find that the frequency of rolling a 5 is higher than the frequencies of other rolls. How can we quantify our confidence that 5 is actually the most probable roll, i.e. that the underlying probability distribution of the die is peaked at the outcome 5?

It seems clear that if we roll less than 6 times, it is pretty random, while if we roll it a thousand times, we can be pretty confident due to the law of large numbers. Furthermore, if it was actually an almost fair die, with only a 1% higher chance to roll 5 than other outcomes, we would need many more rolls to identify this unfairness than if there was a 50% chance of rolling a 5.

However, I have so far been unsuccessful in quantifying these statements. Is there some way to estimate our confidence of identifying the most probable roll in terms of the number of trial rolls we have made, and the relative frequencies they have resulted in?

Answer 1

Using Bayesian methods, you could start with a conjugate Dirichlet prior for the probabilities of the six sides, update it with your observations, and then find the probability from the Dirichlet posterior that side five has the highest underlying probability of the six sides.

This will be affected slightly by the prior you choose and substantially by the actual observations. It may produce slightly counter-intuitive results for small numbers of observations. To take a simpler example with a biased coin,

• if you start with a uniform prior for the probability of it being heads then toss it once and see heads, the posterior probability of it being biased towards heads would be $0.75$ and towards tails $0.25$;
• if instead you tossed it $200$ times and see heads $101$ times then the posterior probability of it being biased towards heads would be about $0.556$;
• if you tossed it $200$ times and see heads $115$ times then the posterior probability of it being biased towards heads would be about $0.983$.

I do not see a simple way of doing the integration with six-sided dice to find the probability a given face is most probable, but simulation will get close enough. The following uses R and a so-called uniform Dirichlet prior for the biases, supposing you observed $21$ dice throws of $2$ ones, $3$ twos, $4$ threes, $5$ fours, $6$ fives and $1$ six:

library(gtools)

probmostlikely <- function(obs, prior=rep(1, length(obs)), 
  cases=10^6) {
  posterior <- prior + obs
  sims <- rdirichlet(cases, posterior)
  table(apply(sims, 1, function(x) which(x == max(x))[1])) / cases
  }

set.seed(2023)
probmostlikely(c(2, 3, 4, 5, 6, 1))

#        1        2        3        4        5        6 
# 0.027885 0.072102 0.152415 0.279509 0.460498 0.007591 

so suggesting that the die is biased most towards five with posterior probability about $0.46$ (and most towards six with posterior probability just under $0.008$).

Seeing that pattern of observations ten times as often would increase the posterior probability that the die is biased most towards five to just under $0.82$ (and reduces those for one and six to something so small that they never appeared as most likely in a million simulations).

probmostlikely(c(20, 30, 40, 50, 60, 10))

#        2        3        4        5 
# 0.000222 0.013702 0.167825 0.818251 

Answer 2

We can use profile likelihood methods to construct a confidence interval for the maximum probability $\theta = \max_{j=1}^k p_k$. Here $p_1, p_2, \dotsc, p_k$ represent the discrete distribution of dice rolls, where in your example $p_5$ is somewhat larger than the others. The results of $n$ dice rolls is given by the random variable $X=(X_1, \dotsc, X_k)$ where in the dice example $k=6$. The likelihood function is then $$ L(p) = p_1^{X_1} p_2^{X_2} \dotsm p_k^{X_k} $$ and the loglikelihood is $$ \ell(p) =\sum_1^k x_j \log(p_j) $$ The profile likelihood function for $\theta$ as defined above is $$ \ell_P(\theta) = \max_{p \colon \max p_j = \theta} \ell(p) $$

With some simulated data we get the following profile log-likelihood function

where the horizontal lines can be used to read off confidence intervals with confidence levels 0.95, 0.99 respectively. I will add the R code use at the end of the post. A paper using bootstrapping for estimating $\theta$ is Simultaneous confidence intervals for multinomial proportions

But this is only a partial solution, in a comment you say

@Dave No, I want to find the most probable outcome of the die.

I read that as finding the maximum probability (done above), but also which of the sides of the dice correspond to the max probability. The Bayesian approach in the answer by user Henry is a direct answer to that. It is not so clear how to approach that in a frequentist way, maybe bootstrapping could be tried? One old approach is subset selection, choosing a subset of the sides of the coin which contains the side with max probability with a certain confidence level. Papers discussing such methods is A subset selection procedure for multinomial distributions and SELECTING A SUBSET CONTAINING ALL THE MULTINOMIAL CELLS BETTER THAN A STANDARD WITH INVERSE SAMPLING.

R code for the plot above:

library(alabama) 

make_proflik_max <- function(x) {
  stopifnot(all(x >= 0))
  k <- length(x) 
  
  Vectorize(function(theta) {
   par <- rep(1/k, k)  # initial values
   fn <- function(p) -sum(log(p) * x)  
   gr <- function(p) -x/p
   hin <- function(p) {  # each component must be positive
     c(p)
   }
   heq <- function(p) {   # must be zero 
     c(sum(p)-1 , max(p) - theta)
   }
   res <- alabama::auglag(par, fn, hin=hin, heq=heq) 
   -res$value
  }   )
}

set.seed(7*11*13) # My public seed

x <- sample(1:6, 200, replace=TRUE, prob=c(9,9,9,9,10,9)) 
            # 5 is a little more probable
x <- table(x)

proflik_max <- make_proflik_max(x)

plot(proflik_max, from=1/6 + 0.001, to=0.35, xlab=expression(theta))
loglik <- function(p) sum(x * log(p))
maxloglik <- loglik(x/200)
mle_pmax <- max(x/200)  

abline(h=maxloglik - qchisq(0.95,1)/2, col="red")
abline(h=maxloglik - qchisq(0.99,1)/2, col="blue")  
abline(v=mle_pmax)