# Distribution of the size of the largest set(s) of faces in roll of multiple dice?

Over the holidays I played a dice game with the family (a gift one got - don't recall the name) that I'll generalize:

Each player has a set of $D$ identical dice with face values $1$ to $F$. The players all roll all their dice, and whatever face shows up the most in their hand, they put all dice with that face aside and roll again. Any dice in the new roll matching the face value of those put aside is added to the aside pile, and the process is repeated until all of some player's dice have been put aside.

For example, say my first roll (using 10 regular dice here) was $\{1,2,1,1,5,6,3,2,4,3\}$: I'd put the three one's aside. If my next roll was $\{5,4,6,3,2,5,1\}$, I'd add the single one to the aside pile, and continue.

First player to get all dice aside wins.

To the question:

On the flight home, I thought about the distribution of $S$, the size of the largest set of same faces in the first roll. In the above example, it was $3$. Had the roll been say $\{1,2,3,1,4,1,5,5,5,4\}$ it would also be $3$ (I can pick either the three ones or the three fives).

For sizes $S>D/2$, it's simply the probability of the binomial trial times $F$, e.g. for $8$ regular six-sided dice the probability of the maximum set size in the first roll being say $5$ is just $6 *C(8, 5) (1/6)^5 (1 - 1/6)^3 = 875/34992$

I'm stuck on the cases for $S\le D/2$ - I think this will require inclusion-exclusion, but I've had no success so far.

Edit: I was able to get what I'm after using ideas from the paper "The exact distribution of the maximum, minimum and the range of Multinomial/Dirichlet and Multivariate Hypergeometric frequencies" by Corrado.

Using the ideas within, I am able to get exact PMF for 100's of dice of 20+ faces quite quickly.

• Suppose you start with $10$ dice, two come up $6$ and no number comes up three times. You put those aside and roll the other $8$. Suppose you get $3$ aces the next roll and no $6$s. Do you keep the two sixes or the three ones? – Ross Millikan Jan 8 '17 at 20:25
• This is a kind of quotient space of a multinomial distribution. You must put together all multinomial probabilities that hold for some $S$. Idk if there is a more easy and direct way to go. – Masacroso Jan 8 '17 at 20:40
• @RossMillikan - The way we played (the standard rules of the game) is once you've set aside dice from the first roll, you're "committed" to that number, so even if you get more of some other number in a later roll than the size of the "aside" pile, you can't switch. The rules had alternate play types, I think one was like your description (Yahtzee like). But I hope my question was not poorly worded - I'm not after some end-game probability result (where things change along the way), but just the simple result for the very first roll. – AppleMyEye Jan 8 '17 at 20:54
• In your most recent edit, for $F=6$ and $D= 8,$ you have $P(S = 5) \approx 0.050.$ The simulation program in my Answer gives the same result (within the margin of simulation error) for the same parameters. – BruceET Jan 9 '17 at 20:20
• Nice article... Unfortunately I dont have a good background on probability/combinatorics to understand completely the article. Anyway I think it would be interesting the comparison between the estimated error through a simulation and the real error, maybe this is a silly question for someone that had a better background on statistics and modeling. – Masacroso Jan 10 '17 at 13:24

Simulation. I don't have anything to contribute to a combinatorial discussion that leads to an analytic solution. But here is a simulation in R statistical software that counts values of $S$ in repeated experiments, each with $D$ independent rolls of a die with $F$ faces. I begin with a demonstration of the simulation method, and then give results of many experiments to simulate the distribution of $S.$ Knowing approximately the correct results may be useful to assess the accuracy of an analytic solution.

Demonstration of concept. The R function sample can be used to simulate rolling a die, and the function rle (for 'run length encoding') can be used as the basis of a method to find $S$ for each experiment. Here are results of experiment yielding various values of $S.$

f=6; d=10;  die=1:f
roll = sample(die, d, repl=T)
srt.roll = sort(roll); srt.roll
## 1 1 2 2 3 3 4 4 5 6
len = rle(srt.roll)$lengths; len ## 2 2 2 2 1 1 # 'run lengths' of sorted rolls s = max(len); s ## 2 # s = 2 (There two 1's, two 2's, etc.) ... ## 1 1 2 2 2 3 4 5 5 6 len = rle(srt.roll)$lengths;  len
## 2 3 1 1 2 1
s = max(len); s
## 3                 # There were three 2's

...
## 2 4 4 4 4 4 5 5 6 6
len = rle(srt.roll)$lengths; len ## 1 5 2 2 s = max(len); s ## 5 # there were five 4's  I hope this is enough to show that I am counting$S$as you intended. Simulation for ten fair six-sided dice. With a million experiments of ten rolls each, the simulated probabilities should be accurate to about three decimal places. In the run shown, possible values of$S$are$2, 3, \dots, 10.$Three is by far the most likely value of$S.$Ten (getting the same face on all ten rolls) is extremely unlikely [$(1/6)^9$], and did not happen to occur during this simulation run. m = 10^6; s = numeric(m); f = 6; d = 10; die = 1:f for (i in 1:m) { srt.roll = sort(sample(die, d, repl=T)) s[i] = max(rle(srt.roll)$lengths)  }
table(s)/m
#  s
##        2        3        4        5        6        7        8        9
## 0.067348 0.529659 0.310519 0.077897 0.012943 0.001536 0.000097 0.000001


Note: If you are familiar with R, you can try out additional runs with different parameter values, to compare with any analytical results you get. With m = 10^6 the program takes a couple of minutes to run, mainly because of the sorting; you might find m = 10^5 gives enough accuracy. If you're not familiar with R, you could easily get it for free from www.r-project.org (Windows, Mac, Linux supported), and just type the program into the Session window.

• This seems the faster and useful approach... The analytic approach seems to rely in some approximations because no closed forms exists for these distributions. – Masacroso Jan 9 '17 at 13:44
• +1. I use R among others. The problem with a sim. is of course when probabilities get small, the number of required iterations for a decent accuracy grows from prohibitive to intractable. See my edit to the question, I found a paper that had ideas I used to code up a routine that gives exact results for any case that's realistic for my needs, and quite quickly. – AppleMyEye Jan 10 '17 at 5:21

This is just an outline for an analytic solution. If we call $S$ as the random variable that count the maximum length of the different groups in a throw of dice then

$$\Pr[S=s]=\frac1{D^n}\sum_{\|\alpha\|_\infty =s}\binom{n}{\alpha}K_\alpha$$

where $n$ is the number of dice thrown, $D$ is the number of faces of each dice, and $\alpha\in\Bbb N_{\ge 0}^D$ and $\|\alpha\|_1=n$ (it is a multi-index).

The main problem seems to be the coefficients $K_\alpha$, that count the variations of each setup $\alpha$, this depends on $n$, $D$ and of course $\alpha$ itself, the big problem is count the multiplicities inside each $\alpha$.

All of this make me thing that the best approach is a simulation, with estimated error and variance, as the previous answer did.

• +1, thanks - see my edit in question, solution found for exact results. – AppleMyEye Jan 10 '17 at 5:19