Probability that exactly y of n rolls of an r-sided die are unique

Question

Consider an $r$-sided die that is rolled $n$ times. What is the probability that of the $n$ rolls exactly $y$ of the rolls are unique?

For example, consider $n = 2$ and $r = 3$. The possibilities are

0 0
0 1
0 2
1 0
1 1
1 2
2 0
2 1
2 2

In this case, $P(y = 0) = 1/3$, $P(y = 1) = 0$, and $P(y = 2) = 2/3.$

I feel as though this should be a fairly simple problem to solve, but for the life of me I haven't been able to figure it out (I need this information to bound the behavior of an algorithm I'm working on).

Answer 1

There is an efficient, simple, $O(n)$ solution.

In expanding the polynomial

$$f_{n,r} = \left(x_1+x_2+\cdots+x_r\right)^n = \sum_{i_1,i_2,\ldots,i_r} \binom{n}{i_1,i_2,\ldots,i_r} x_1^{i_1}x_2^{i_2}\cdots x_r^{i_r},$$

for each of the $\binom{r}{y}$ subsets of $y$ of the variables there will be a term like this one

$$\binom{n}{1,1,\ldots,1,i_{y+1},\ldots, i_r}\left(x_1x_2\cdots x_y\ x_{y+1}^{i_{y+1}}\cdots x_r^{i_r}\right)$$

whose coefficient gives the number of times at least $y$ of the variables appear just once. This coefficient can be found by differentiating $f_{n,r}$ with respect to each of those $y$ variables, setting the values of those variables to $0$, and setting the values of the remaining $r-y$ variables to $1$, because

$$\frac{\partial^y}{\partial x_1\partial x_2\cdots \partial x_y}\left(x_1x_2\cdots x_y\ x_{y+1}^{i_{y+1}}\cdots x_r^{i_r}\right) = x_{y+1}^{i_{y+1}}\cdots x_r^{i_r}$$

evaluates to $1$ and all other terms have at least one of the first $y$ variables as a factor, whence they evaluate to $0$.

Computing this derivative for the original expression of $f_{n,r}$ yields (using the falling factorial notation for the coefficient)

$$\eqalign{&\frac{\partial^y}{\partial x_1\partial x_2\cdots \partial x_y} \left(x_1+x_2+\cdots+x_r\right)^n \\&= n(n-1)\cdots(n-y+1)\left(x_1+x_2+\cdots+x_r\right)^{n-y} \\ &= n_{(y)}\left(x_1+x_2+\cdots+x_r\right)^{n-y}. }$$

When $y$ of the $x_i$ equal $0$ and the remaining $r-y$ equal $1$, the right hand side evaluates to

$$ n_{(y)}(r-y)^{n-y}.$$

Multiplying by $\binom{r}{y}$ to account for all possible combinations of $y$ variables and applying the Principle of Inclusion Exclusion ("PIE") produces the number of times $y$ variables appear exactly once, which is

$$\binom{r}{y}\sum_{j=y}^{\min(r,n)} (-1)^{j-y} (r-j)^{n-j} n_{(j)}\binom{r-y}{j-y}.$$

Dividing this by $r^n$ gives the associated probabilities. The computational effort is $O(\min(r,n)-y)$.

Nothing comes for free! As in most applications of the PIE, this is an alternating sum of terms that can vary radically in size, with the final result being much smaller than the largest terms. There can be catastrophic loss of precision, so high-precision (or, better yet, exact rational) arithmetic is needed. With that available, the implementation is remarkably short. Here it is in Mathematica:

p[n_, k_] := n^k; p[n_, 0] := 1;
f[n_, d_, k_] := Binomial[d,k]  
            Sum[(-1)^(j-k) Binomial[d-k,j-k] FactorialPower[n,j] p[d-j,n-j],{j,k,Min[d,n]}]

As an example, let's plot the full distribution for a particular $n$ and $r$:

With[{n = 100, r = 100}, DiscretePlot[f[n, r, y]/r^n, {y, 0, Min[n, r]}]]

---

As an example, consider the case $n=4$, $r=3$, and values of $y$ from $0$ through $3$. The expansion of $f_{4,3}$ is

$$x_1^4+x_2^4+x_3^4 \\\color{blue}{+4 x_2 x_1^3+4 x_3 x_1^3+4 x_1 x_2^3+4 x_3x_2^3+4 x_1 x_3^3+4 x_2x_3^3}\\+6 x_2^2 x_1^2+6 x_2^2 x_3^2+6 x_3^2 x_1^2\\\color{red}{+12 x_1 x_2 x_3^2 +12 x_1 x_2^2 x_3 +12 x_1^2 x_2 x_3}.$$

Consider the calculation for $y=1$.

• The terms with exactly one $x_1$ in them are

  $$\color{blue}{4 x_1 x_2^3+4 x_1 x_3^3}+\color{red}{12 x_1 x_2 x_3^2 +12 x_1 x_2^2 x_3}.$$

  These coefficients sum to $4+4+12+12=32$. Thus we would estimate the total number of terms with just a single one of the $x_i$ would be $3\times 32=96$.

• We're not done yet. The terms with one $x_1$ and one other $x_i$ in them are

  $$\color{red}{12 x_1 x_2 x_3^2 + 12 x_1 x_2^2 x_3}.$$

  This tells us that when we counted the $x_1$ terms previously, we overcounted by $12+12 = 24$. The total overcount therefore is $3\times 24 = 72$.

• We are done now, because there are no terms possible with exactly one instance of three sides.

Consequently, the count for $y=1$ is

$$96 - 72 + 0 = 24.$$

Indeed, this is the sum of coefficients of

$$\color{blue}{4 x_2 x_1^3+4 x_3 x_1^3+4 x_1 x_2^3+4 x_3x_2^3+4 x_1 x_3^3+4 x_2x_3^3}.$$

Answer 2

Define $F(n, k)$ to be the number of ways to allocate $k$ options to $n$ flips such that each option appears either 0 or $\geq 2$ times. Then the probability that you see exactly $y$ unique values when you roll a $k$-sided dice $n$ times is:

$$ Pr(Y=y) = \frac{{k\choose y}{n\choose y}y!F(n-y, k-y)}{k^n} $$

Basically, there are ${k \choose y}$ ways to select the $y$ unique options from all $k$ options, ${n\choose y}$ ways to select the $y$ rolls for these unique options, and $y!$ orderings of the $y$ options within these rolls.

All that remains is to compute $F(n, k)$. There are a few simple cases and then a recursive definition:

\begin{align*} F(0, k) &= 1 &\forall~k\geq 0 \\ F(1, k) &= 0 &\forall~k\geq 0 \\ F(n, 0) &= 0 &\forall~n\geq 1 \\ F(n, k) &= F(n, k-1) + \sum_{i=2}^n {n\choose i}F(n-i, k-1) &\forall~n\geq 2, k\geq 1 \end{align*}

The recursive step selects an arbitrary option and separately considers the number of allocations for which it appears $0, 2, 3, \ldots, n$ times. This formulation enables the calculation of the entire pmf in $O(n^2k)$ runtime, which should be a good deal more efficient than summing over all valid partitions of the multinomial distribution. Here's an R implementation:

uniquePMF <- function(n, k) {
  F <- matrix(0, nrow=n+1, ncol=k+1)
  F[1,] <- 1
  for (.k in 1:k) {
    for (.n in 2:n) {
      F[.n+1,.k+1] <- F[.n+1,.k] + sum(choose(.n, 2:.n)*F[.n-(2:.n)+1,.k])
    }
  }
  out <- sapply(0:min(n, k), function(y) choose(k, y)*choose(n, y)*factorial(y)*F[n-y+1,k-y+1]) / k^n
  names(out) <- 0:min(n, k)
  out
}

This returns your hand-calculated results for the $n=2, k=3$ case:

uniquePMF(2, 3)
#         0         1         2 
# 0.3333333 0.0000000 0.6666667

It can also comfortably handle larger instances (here $n=k=100$):

plot(0:100, uniquePMF(100, 100), xlab="y", ylab="Pr(Y=y)")

Answer 3

This is basically a generalized coupon-collector problem. I don't think you're going to find an easy closed form solution to it. In general, the distribution of dice throws is multinomial:

$$\binom{n}{a_1,a_2,\cdots,a_r}(1/r)^n,$$

where $a_i$ represents the number of times $i$ is rolled, and $a_1+\cdots+a_r=n$. Then the answer is:

$$(1/r)^n\sum\binom{n}{a_1,a_2,\cdots,a_r},$$

where the sum is over all combinations of $a_i$ such that exactly $y$ of them are nonzero. I.E. you are looking at all partitions of $n=a_1+\cdots+a_r$ such that exactly $y$ are nonzero. You can use some symmetry to simplify this further by just summing over all combinations of $n=a_1+\cdots+a_y$ with $a_i>0$ for $i=1,2,\cdots,y$ and $a_i=0$ for $i>y$, and then accounting for their multiplicity via $\binom{n}{a_1}\binom{n-a_1}{a_2}\cdots$.