# Roll an N sided die K times. Let S be the side that appeared most often. What is the expected number of times S appeared?

For example, consider a 6 sided die rolled 10 times. Based on the following monte-carlo simulation, I get that the side that appears most will appear 3.44 times on average.

n = 6
k = 10
samples = 10000
results = []

for _ in range(samples):
counts = {s:0 for s in range(n)}
for _ in range(k):
s = randint(0, n-1)
counts[s] += 1

results.append(max(counts.values()))

print sum(results)/float(len(results))


But I can't figure out how to get this in a closed form for any particular N and K.

• Here is an attempt at formalizing your question: Roll an $N$-sided die $K$ times. The probability that a face with value $n\in\{1,2,\dots,N\}$ is hit $k_n$ times is $\binom{K}{k_n}N^{-k_n}(1-N^{-1})^{K-k_n}.$ It is required that $\sum_{n=1}^N k_n=K$ and we seek $\mathrm{E}[\mathrm{max}(k_n)].$ – Bobson Dugnutt Jan 25 '17 at 23:29
• I believe you are looking for this post: mathoverflow.net/questions/146418/… – Jan Jan 26 '17 at 0:00
• If you try for $N=2$ and $k=2j+1$ odd, I get an expected value of $$\frac{2j+1}{2}\left(1+\frac{\binom{2j}j}{2^{2j}}\right)$$ which is approximately $\frac{k}{2}+\sqrt{\frac{k}{2\pi}}$ – Thomas Andrews Jan 26 '17 at 0:03
• Sorry, I misread the question! – bof Jan 26 '17 at 2:08

Not an answer, but it's worth trying some example with small $N$.

Let $M$ be your number.

For $n=2$, and any $k/2 < m \leq k$ you have $P(M=m)=\frac{\binom{k}{m}}{2^{k-1}}$. If $k$ even, then $P(M=k/2)=\frac{\binom{k}{k/2}}{2^k}$.

So, for $k=2j+1$ odd, you have:

\begin{align}E(M)&=\frac{1}{2^{2j}}\sum_{m=j+1}^{2j+1} m\binom{2j+1}{m}\\ &=\frac{2j+1}{2^{2j}}\sum_{m=j+1}^{2j+1}\binom{2j}{m-1}\\ &=\frac{2j+1}{2^{2j}}\sum_{m=j}^{2j}\binom{2j}{m}\\ &=\frac{2j+1}{2^{2j}}\left(\frac{1}{2}2^{2j}+\frac{1}{2}\binom{2j}{j}\right)\\ &=\frac{k}{2}\left(1+\frac{\binom{k-1}{(k-1)/2}}{2^{k-1}}\right) \end{align}

For $N=2$ and $k=2j$, you get: \begin{align}E(M)&=j\frac{\binom{2j}{j}}{2^{2j}}+\frac{1}{2^{2j-1}}\sum_{m=j+1}^{2j}m\binom{2j}{m}\\ &=j\frac{\binom{2j}{j}}{2^{2j}}+\frac{2j}{2^{2j-1}}\sum_{m=j+1}^{2j}\binom{2j-1}{m-1}\\ &=j\frac{\binom{2j}{j}}{2^{2j}}+\frac{j}{2^{2j-2}}\cdot\frac{2^{2j-1}}{2}\\ &=\frac{k}{2}\left(1+\frac{\binom{k}{k/2}}{2^k}\right) \end{align}

It's gonna get messier when trying it for $N=3$.

Here is a closed form, we will need more sophisticated methods for the asymptotics.

Following the notation introduced at this MSE link we suppose that the die has $m$ faces and is rolled $n$ times. Rolling the die with the most occured value being $q$ and instances of this size being marked yields the species

$$\mathfrak{S}_{=m} (\mathfrak{P}_{=0}(\mathcal{Z}) + \mathfrak{P}_{=1}(\mathcal{Z}) + \cdots + \mathcal{V}\mathfrak{P}_{=q}(\mathcal{Z})).$$

This has generating function

$$G(z,v) = \left(\sum_{r=0}^{q-1} \frac{z^r}{r!} + v\frac{z^q}{q!}\right)^m.$$

Subtracting the values where sets of size $q$ did not occur we obtain the generating function

$$H_{q}(z) = \left(\sum_{r=0}^{q} \frac{z^r}{r!}\right)^m - \left(\sum_{r=0}^{q-1} \frac{z^r}{r!}\right)^m.$$

This also follows more or less by inspection.

We then obtain for the desired quantity the closed form

$$\bbox[5px,border:2px solid #00A000]{ \frac{n!}{m^n} [z^n] \sum_{q=1}^n q H_q(z).}$$

Introducing

$$L_{q}(z) = \left(\sum_{r=0}^{q} \frac{z^r}{r!}\right)^m$$

we thus have

$$\frac{n!}{m^n} [z^n] \sum_{q=1}^n q (L_{q}(z) - L_{q-1}(z)).$$

This is

$$\frac{n!}{m^n} [z^n] \left(n L_n(z) - \sum_{q=0}^{n-1} L_q(z)\right).$$

We also have for

$$[z^n] L_q(z) = \sum_{k=0}^{\min(q, n)} \frac{1}{k!} [z^{n-k}] \left(\sum_{r=0}^{q} \frac{z^r}{r!}\right)^{m-1}$$

Furthermore we obtain for $m=1$

$$[z^n] L_q(z) = [[n \le q]] \times \frac{1}{n!}.$$

With these we can implement a recursion, which in fact on being coded proved inferior to Maple's fast polynomial multiplication routines. It is included here because it memoizes coefficients of $L_q(z)$, thereby providing a dramatic speed-up of the plots at the cost of allocating more memory.

All of this yields the following graph where we have scaled the plot by a factor of $n/m.$ This is it for a six-sided die:

  3+  H
+
|  H
+   H
|
+    H
+     HH
|      HH
2+        HH
|         HHH
+            HHHH
+                HHHHHHH
|                       HHHHHHHHHHH
+                                  HHHHHHHHHHHHHHHHHHHHHHHHH
|                                                           HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
-+--+---+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
0              20             40             60             80             100            120


And here is the plot for a twelve-sided die. (I consider it worth observing that we have the exact value for the expectation in the case of $120$ rolls of this die, a case count that has $130$ digits, similar to what appeared in the companion post.)

  8+
|
+
+
|
+
+ H
|
6+
|
+
+
|  H
+
+  H
|
4+   HH
+     H
|      H
+      HH
+        HH
|          HHHH
+              HHHHHH
2+                    HHHHHHHHHHH
|                               HHHHHHHHHHHHHHHHHHHHHHHHH
+                                                        HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
-+--+---+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
0              20             40             60             80             100            120


This was the Maple code.

with(plots);
with(combinat);

ENUM :=
proc(n, m)
option remember;
local rolls, res, ind, counts, least, most;
res := 0;

for ind from m^n to 2*m^n-1 do
rolls := convert(ind, base, m);

counts := map(mel->op(2, mel),
convert(rolls[1..n], multiset));

res := res + max(counts);
od;

res/m^n;
end;

L := (m, rmax) -> add(z^r/r!, r=0..rmax)^m;

X :=
proc(n, m)
option remember;
local H;

H := q -> expand(L(m,q)-L(m,q-1));

n!/m^n*
end;

LCF :=
proc(n,m,q)
option remember;

if n < 0 then return 0 fi;

if m = 1 then
if n <= q then return 1/n! fi;
return 0;
fi;

end;

LVERIF :=
(m, q)  -> add(LCF(n, m, q)*z^n, n=0..q*m);

XX :=
proc(n, m)
option remember;
local res;

res :=

res*n!/m^n;
end;

DICEPLOT :=
proc(nmx, m)
local pts;

pts := [seq([n, XX(n,m)/(n/m)], n=1..nmx)];
pointplot(pts);
end;

• Very nice! (+1) – Markus Scheuer Jan 27 '17 at 9:12