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

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?

• Are you looking for a confidence interval for the probability of a "5"? Or (closely related) for hypothesis testing of the null hypothesis that $P("5")=1/6$? Mar 2 at 12:52
• I am looking for a confidence interval for the probability of each roll, which would then allow me to evaluate how sure I am that the most rolled result is also the most probable result. I am not interested in evaluating whether the die is fair. Thanks for the clarification help. Mar 2 at 12:59
• You could create a confidence interval for the estimate of the probability of each of the six sides but it sounds like you want to test a hypothesis that side a is the side with the most probability against all the other five sides? Mar 2 at 13:09
• @Dave No, I want to find the most probable outcome of the die. Mar 2 at 13:16
• For essentially the same question (with answers) please see stats.stackexchange.com/questions/274211. There the metaphor for a discrete random variable is a population of voters for various candidates -- but it's exactly the same problem.
– whuber
Mar 5 at 17:04

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)))) / 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

• Can you please elaborate a bit on the theory behind your computation? How do you compute or obtain a formula for $P(p_i>p_j\mbox{ for all }j\neq i \mid obs)$? Mar 3 at 8:53
• @cdalitz each of the million sims generates a random value for the probabilities of each side from the posterior Dirichlet distribution. For example, the first one is 0.1201143 0.1533349 0.2085550 0.2718155 0.2220087 0.0241715 (note these add up to $1$ as they should), and the largest of these is the 4th value. It does this a million times and counts how often each face is the most probable from this sample, reporting the frequency for each, giving a close approximation to the probabilities each face is the most probable in the posterior Dirichlet distribution. Mar 3 at 9:08
• Thanks for explaining the simulation. Actually I had to look up the term "conjugate prior" and learn, what parameter set of the Dirichlet distribution corresponds to the uniform distribution and how the prior parameters propagate to the posterior parameters. This helped me understand your answer, but it nevertheless assumes a lot of knowledge in probability theory, which might be beyond the grasp of the original poster. It might be helpful to elaborate somewhat on the underlying ideas in the answer so that more readers can understand it. Mar 3 at 9:38
• @cdalitz I tried to warn of the need for some knowledge in probability theory with my introductory words "Using Bayesian methods". I also tried to explain some of the issues with the simplified example of a biased coin. It would be possible to expand this answer into a full lecture including covering uniformity over a standard $(K-1)$-simplex in multidimensional geometry, as well as conjugate distributions and Monte Carlo methods, but it would have got too long. Mar 3 at 9:51

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)