Most likely counts of outcomes of a multinomial after N trials (mode of multinomial) Given a categorical random variable $X = \{x_1,x_2,...,x_k\}$ with corresponding probabilities of each outcome $P = \{p_1, p_2,...,p_k\}$ (i.e. a multinomial):

*

*What are the most likely counts of outcomes aftern $N$ trials?
Let's call the vector of counts after $N$ trials a state: $C_N = \{c_1^N, c_2^N,...,c_k^N\}$. What is the most likely observable state after $N$ trials? In an observable state counts must be integers.
If you just multiply $x_i \times p_i$ you get fractional counts, that doesn't work. Is there a better strategy than computing $x_i \times p_i$ values and then trying all roundings where $\sum_i c_i^N = N$?

*What is the next most likely state (the one that's a little less likely than the best one)?
An example:
Outcome probabilities: $[p_1=0.60, p_2=0.30, p_3=0.10]$.
Number of trials: $N=7$.
Most likely state for $N=7$ trials: $[c_1^7=4.20, c_2^7=2.10, c_3^7=0.70]$.
If you just round it to: $C_{round}=[4, 2, 1]$, you get $P(C_{round})=P([4, 2, 1])=0.1225$, where $P(x)$ is a probability function.
However, $P([5, 2, 0])=0.1470$. So just rounding doesn't work. I guess because of the multinomial coefficient getting in the way of computing probabilities.
 A: You can create an algorithm to compute the mode relatively easily in this problem.  To facilitate this analysis, define the function:
$$H_{i,\ell} \equiv \log(\ell) - \log(p_i).$$
The log-mass for the multinomial distribution can be written as:
$$\begin{align}
\log \text{Mu}(\mathbf{x}|N,\mathbf{p})
&= \log(N!) + \sum_{i=1}^k [x_i \log(p_i) - \log(x_i!)] \\[6pt]
&= \log(N!) + \sum_{i=1}^k [x_i \log(p_i) - \log(x_i) - \log(x_i-1) - \cdots - \log(2)] \\[6pt]
&= \log(N!) - \sum_{i=1}^k \sum_{\ell = 1}^{x_i} H_{i,\ell}. \\[6pt]
\end{align}$$
From this form we can see that the mass function is maximised when we allocate values into categories in a way that minimises the quantity $H_{i,\ell}$ corresponding to each allocation.  In the case where the values in $N \mathbf{p}$ are all integers (i.e., the probabilities are all multiples of $\tfrac{1}{N}$) the mode occurs at the point $\hat{\mathbf{x}} = N \mathbf{p}$ --- at this point we have $H_{i,x_i} = \log(N)$ for each $i=1,....k$.  In the more general case we have $\hat{\mathbf{x}} \approx N \mathbf{p}$, but the exact value of the mode requires some search around this approximating point.  Points that are next most likely, etc., can be determined by looking at the values of the $H$ function here.
Here is a simple algorithm in R to compute the mode, where we allocate each of the $N$ values to categories one at a time by minimising the $H$ values (with tie-breakers being given by the lowest category number).  The function takes in parameters size and prob and outputs a vector that is a mode of the distribution.
mode.multinom <- function(size, prob) {
  
  #Set initial values
  k <- length(prob)
  H <- -log(prob)
  X <- integer(k)
  prob <- prob/sum(prob)
  
  #Compute mode allocation
  for (i in 1:size) {
    Z    <- which.min(H)
    X[Z] <- X[Z]+1
    H[Z] <- log(X[Z]+1) - log(prob[Z]) }
  
  #Give output
  X }

Here is an example where we use this algorithm to compute the mode
#Set parameters
N    <- 76
PROB <- c(0.01, 0.13, 0.03, 0.10, 0.22, 0.05, 0.07, 0.32, 0.02, 0.05)

#Check the approximate mode
N*PROB
[1]  0.76  9.88  2.28  7.60 16.72  3.80  5.32 24.32  1.52  3.80

#Compute the mode
MODE <- mode.multinom(size = N, prob = PROB)
MODE
[1]  0 10  2  8 17  4  5 25  1  4

#Compute the log-mass at the mode
dmultinom(MODE, size = N, prob = PROB, log = TRUE)
[1] -14.06952

A: You could start with a solution rounded down. For instance in the example you start with [4,2,0]. Then you iteratively keep adding the remainder (from rounding down) in the places that are making the best improvement.
This is my intuition. I still have to prove that this works.
