Expected value of the largest item in a multinomial distribution In this question, I use the notation on wikipedia.
Suppose $X=(X_1,\ldots, X_k)$ follows a multinomial distribution with parameter $n,\mathbf p$, where $n$ is the sample size or the number of trials, and $\mathbf p$ has $k$ components, each of which has equal value $1/k$, so each outcome has the same probability.
Let $M=\max\{X_1,\ldots, X_k\}$, which is another random variable. How can I find the expected value $\mathbb E(M)$? Are there any know results about the value of $E(M)$?
This is a related post.
 A: $\newcommand{\p}{\mathbf p}$The paper "Exact Tail Probabilities and Percentiles of
the Multinomial Maximum" by Anirban DasGupta provides an asymptotic answer. This paper contains the following result (I changed the notation slightly):

Let $X_n\sim\text{Mult}(n, \p_n)$ with $\p_n = (1/K_n, \dots, 1/K_n)$ and suppose that $n, K_n\to\infty$ with $\frac{n}{K\log K}\to\infty$. Define $\mu = \frac nK$ and
$$
w = \frac{\log K - \frac 12 \log\log K}{\mu}.
$$
Let $\epsilon$ be the unique positive root of
$$
(1+\epsilon)\log(1+\epsilon) - \epsilon = w.
$$
Then
$$
P\left(\frac{X_\max - \mu(1+\epsilon)}{\sqrt{\frac{n}{2K\log K}}} + \frac 12 \log 4\pi \leq z\right) \to e^{-e^{-z}}
$$

which is the CDF of a standard Gumbel distribution.
This means that
$$
(\gamma - \frac 12\log K)\sqrt{\frac{n}{2K\log K}} + \mu(1+\epsilon)
$$
approximates the expected value of $X_\max$, where $\gamma \approx 0.577$ is the Euler–Mascheroni constant.
Simulating to check:

The vertical red line in the second histogram is the approximate mean.
Code:
set.seed(111)
K <- 1000
pr <- rep(1/K,K)
n <- ceiling(K^2 * log(K))  # I'm using this to keep the required asymptotics correct
nsim <- 10000

mu <- n/K
w <- (log(K) - .5 * log(log(K))) / mu
eps <- uniroot(function(e) {(1+e) * log(1+e) - e - w}, lower=0, upper=1000)$root

Xmax <- apply(rmultinom(nsim, n, pr), 2, max)

par(mfcol=c(2,1))
hist((Xmax - mu * (1 + eps)) / sqrt(n / ( 2 * K * log(K))) + .5 * log(4 * pi), 50, freq=F, main = "Transformed distribution with Gumbel pdf", xlab = "x")
curve(exp(-(x + exp(-x))), from=-50, to=50, n=1000, add=T, col="red", lwd=2)
hist(Xmax, 50, main = "Simulated maxima of Mult(n, 1/K,...,1/K)", xlab="x")
abline(v = (-digamma(1) - .5 * log(4 * pi)) * sqrt(n / (2 * K * log(K))) + mu * (1 + eps), col=2,lwd=2)

