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?
 A: 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 

A: 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.
(It is too late now, I will add details tomorrow)
