# Probability of empty urns for undistinguishable balls in distinguishable urns

Distinguishable balls in distinguishable urns is a well known problem in probability theory and can be easily simulated through Montecarlo simulations. Here we provide our code for n balls and m urns:

n = 20;
m = 10;
mvec = Range[1, m];
nrun=10000;
monteCarlo=Table[RandomChoice[mvec, n], {i, 1, nrun}];


The probability of finding empty urns in this case is given by the normalized histogram of the following vector:

pE1=Map[Length,Table[Complement[mvec, Union[RandomChoice[mvec, n]]], {mrun}]]


If we now want to now count the same probability, but for UNDISTINGUISHABLE balls, then the only solution we found to this counting problem is:

allConfig = ParallelTable[Split[Sort[monteCarlo[[i]]]], {i, 1, nrun}];
config = Union[allConfig];
emptysites = ParallelTable[Complement[mvec, Flatten[config[[i]]]], {i, 1, Length[config]}];


And the probability of finding empty urns in this case is given by the normalized histogram of the following vector:

   pE2 = ParallelMap[Length, emptysites];


The problem is that even for not too large m,n the number of possible configuration is huge (Binomial[n + m - 1, m - 1]), and for an nrun that is very high our code is very slow. Anyone can find a fastest solution?

• This is quite difficult to follow. What is shown in the example? Which are the balls and which are the urns? What you might be looking for is something along the lines of draws = 1000; monteCarloUD = Map[Differences, Table[Join[{0}, Sort[RandomChoice[Range[m + n - 1], n - 1]], {m + n - 1}], {draws}]]; Histogram[Map[Count[#, 0] &, monteCarloUD]] – Daniel Lichtblau Nov 20 '18 at 15:44
• There are explicit formulas for the probabilities of the number of empty urns for both models that you describe. Here's one such reference: pdfs.semanticscholar.org/14f5/…. – JimB Nov 20 '18 at 15:53

For distinguishable balls there is a formula for the probability of the number of empty urns:

prDistinguishable[nEmptyUrns_, nBalls_, nUrns_] :=
If[0 <= nEmptyUrns < nUrns,
StirlingS2[nBalls, nUrns - nEmptyUrns] FactorialPower[nUrns, nUrns - nEmptyUrns]/
nUrns^nBalls, 0]


(See Mingshu 2011, for example.) We can create a table of probabilities:

nBalls = 20;
nUrns = 10;
(probabilities =
Table[{k, prDistinguishable[k, nBalls, nUrns]}, {k, 0, nUrns - 1}]) // N // TableForm


$$\left( \begin{array}{cc} 0. & 0.214737 \\ 1. & 0.435865 \\ 2. & 0.275261 \\ 3. & 0.0673962 \\ 4. & 0.00651079 \\ 5. & 0.00022656 \\ 6. & \text{2.2796986414104\grave{ }*{}^{\wedge}-6} \\ 7. & \text{4.1803664112\grave{ }*{}^{\wedge}-9} \\ 8. & \text{4.718583\grave{ }*{}^{\wedge}-13} \\ 9. & \text{1.\grave{ }*{}^{\wedge}-19} \\ \end{array} \right)$$

For indistinguishable balls one can find the probability of the number of empty urns:

prIndistinguishable[nEmptyUrns_, nBalls_, nUrns_] :=
If[IntegerQ[nEmptyUrns] && IntegerQ[nBalls] && IntegerQ[nUrns] && 0 <= nEmptyUrns < nUrns,
Binomial[nBalls - 1, nBalls - (nUrns - nEmptyUrns)] Binomial[nUrns, nUrns - nEmptyUrns]/
Binomial[nBalls + nUrns - 1, nUrns - 1], 0]


from which we can create a table of probabilities:

(probabilities =
Table[{k, prIndistinguishable[k, nBalls, nUrns]}, {k, 0, nUrns - 1}]) // N // TableForm


$$\left( \begin{array}{cc} 0. & 0.00922396 \\ 1. & 0.0754688 \\ 2. & 0.226406 \\ 3. & 0.325096 \\ 4. & 0.243822 \\ 5. & 0.0975289 \\ 6. & 0.0203185 \\ 7. & 0.00204893 \\ 8. & 0.0000853719 \\ 9. & \text{9.985017481269356\grave{ }*{}^{\wedge}-7} \\ \end{array} \right)$$

There's no need to perform simulations to obtain the probabilities. If you want to samples from the distribution, then sampling from the resulting probability distribution would be more straight forward (and definitely speedier). For example:

ProbabilityDistribution[(Binomial[10, 10 - \[FormalX]] Binomial[19, 10 + \[FormalX]])/10015005, {\[FormalX], 0, 9, 1}]
x = RandomVariate[indistinguishable, 100000];


I don't think your current code for the indistinguishable case produces the correct distribution as the resulting distribution is nearly identical to the distinguishable case.