Expected number of empty bins with clumpy balls?

I have $b$ bins and a reservoir of $n$ balls.

The balls have a tendency to "clump", that is, when I try to grab one from the reservoir, the number of balls removed is uniformly random on $[1, m]$, where $m$ is the number of balls remaining in the reservoir for that draw.

The "clump" is added to one of the $b$ bins, the bin selected from all bins with equal probability.

When the reservoir is depleted, how many empty bins are expected?

I don't even know where to start...

Edit: I asked a neighbor (an actuary) this question, found this on a note on my windshield this morning, seems to match the answer by Quasi :

$\frac{b\prod _ {z=\frac{b-1}{b}}^{\frac{b-1}{b}+n-1} z}{n!}$

• Are you expecting a closed form? If so, why? – quasi Sep 10 '17 at 0:33
• Does the source problem give actual values for $n$ and $b$? – quasi Sep 10 '17 at 0:44
• @quasi my math stopped at college algebra, I can do most simple combinatorics (like simple probability, coupon colletor, etc.). Had to wikipedia "closed form" to make sure I remembered it correctly. So, I don't know if a simple form might exist. Say for n=10000 and b=100. – Tommy Sep 10 '17 at 1:07
• This is a family web site you know. – zhw. Sep 10 '17 at 1:33

For positive integers $n,b$, let $f(n,b)$ be the expected number of empty bins at the end of the process.

I'm not sure if there's a closed form for $f(n,b)$, but here's a recursion, implemented in Maple, to compute $f(n,b)$ for given (not too large) values of $n,b$, . . .

Shown above is a sample calculation for $n=10,\;b = 6$, which shows that $$f(10,6) = \frac{751583152441}{208971104256} \approx 3.596588893$$ For $n=10000,\;b=100$, the recursion fails (too many levels of recursion for my version of Maple), but a simulation gives the approximate result $$f(10000,100) \approx 90.67$$

• Thank you! See edit to my original question for pretty simple form. – Tommy Sep 10 '17 at 22:57
• The product formula works perfectly! – quasi Sep 10 '17 at 23:47

Credit goes to @quasi for the first answer and the recurrence. Using generating functions and classifying on the number of rounds $k$ we find the PGF

$$\bbox[5px,border:2px solid #00A000]{ G(u) = \frac{1}{n} \sum_{k=1}^n [w^{k-1}] \prod_{q=1}^{n-1} \left(1+\frac{w}{q}\right) \times \frac{k!}{b^k} [z^k] (u-1+\exp(z))^b}$$

where the coefficient on $[u^p]$ is the probability of obtaining $p$ empty bins in the experiment with the given $n$ and $b.$ Here we have used the labeled combinatorial species $\mathfrak{S}_{=b}(\mathcal{U}+\mathfrak{P}_{\ge 1}(\mathcal{Z}))$ to construct the EGF that was used. The desired expectation is then given by

$$\left. \frac{d}{du} G(u)\right|_{u=1.}$$

The derivative yields

$$\left. \frac{1}{n} \sum_{k=1}^n [w^{k-1}] \prod_{q=1}^{n-1} \left(1+\frac{w}{q}\right) \times \frac{k!}{b^k} [z^k] b (u-1+\exp(z))^{b-1}\right|_{u=1.} \\ = \frac{b}{n} \sum_{k=1}^n [w^{k-1}] \prod_{q=1}^{n-1} \left(1+\frac{w}{q}\right) \times \frac{k!}{b^k} [z^k] \exp((b-1)z) \\ = \frac{b}{n} \sum_{k=1}^n \frac{(b-1)^k}{b^k} [w^{k-1}] \prod_{q=1}^{n-1} \left(1+\frac{w}{q}\right) \\ = \frac{b-1}{n} \sum_{k=1}^n \frac{(b-1)^{k-1}}{b^{k-1}} [w^{k-1}] \prod_{q=1}^{n-1} \left(1+\frac{w}{q}\right) \\ = \frac{b-1}{n} \prod_{q=1}^{n-1} \left(1+\frac{b-1}{bq}\right) \\ = \frac{b-1}{n!} \prod_{q=1}^{n-1} \left(q+\frac{b-1}{b}\right) = \frac{b}{n!} \prod_{q=0}^{n-1} \left(q+\frac{b-1}{b}\right).$$

This is

$$b \times {(b-1)/b+n-1\choose n}$$

or indeed

$$\bbox[5px,border:2px solid #00A000]{ b \times {n-1/b\choose n}}$$

as claimed in edit to question by OP. (This formula will produce zero when there is just one bin and this is the correct value, no empty bins are possible in this case.)

Here is some Maple code to help explore this interesting problem. Using the multiplicative version of the closed form gives

$$\bbox[5px,border:2px solid #00A000]{90.66863442}$$

for $n=10000$ and $b=100.$

with(combinat);

ENUM_GF :=
proc(n, b)
option remember;
local recurse, gf;

gf := 0;

recurse :=
proc(prob, rest, cc)
local rm;

if rest = 0 then
gf := gf +
prob*cc!*coeftayl((u-1+exp(z))^b, z=0, cc)/b^cc;

return;
fi;

for rm to rest do
recurse(prob/rest, rest-rm, cc+1);
od;
end;

recurse(1, n, 0);

gf;
end;

X_GF := (n, b) -> 1/n*
*k!*coeftayl((u-1+exp(z))^b, z=0, k)/b^k,
k=1..n);

ENUM := (n, b) -> subs(u=1, diff(ENUM_GF(n, b), u));

X := (n, b) -> b*binomial(n-1/b, n);
Y := (n, b) -> b*mul(n-1/b-q, q=0..n-1)/n!;


We also have a very basic simulation, which confirmed the closed form on all values that were examined.

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <time.h>

int main(int argc, char **argv)
{
int n = 6 , b = 3, trials = 1000;

if(argc >= 2){
n = atoi(argv[1]);
}

if(argc >= 3){
b = atoi(argv[2]);
}

if(argc >= 4){
trials = atoi(argv[3]);
}

assert(1 <= n);
assert(1 <= b);
assert(1 <= trials);

srand48(time(NULL));
long long data = 0;

for(int tind = 0; tind < trials; tind++){
int rest = n;

int bins[b];

for(int bidx = 0; bidx < b; bidx++)
bins[bidx] = 0;

while(rest > 0){
int rm = 1 + drand48() * (double)(rest);
int bidx = drand48() * (double)(b);

bins[bidx] += rm;

rest -= rm;
}

int empty = 0;
for(int bidx = 0; bidx < b; bidx++)
if(!(bins[bidx]))
empty++;

data += empty;
}

long double expt = (long double)data/(long double)trials;
printf("[n = %d, b = %d, trials = %d]: %Le\n",
n, b, trials, expt);

exit(0);
}


• Thank you for such a detailed answer. Way over my head, indistinguishable from magic for me, but appreciated! – Tommy Sep 11 '17 at 23:56
• Thank you, it was an interesting question and it will perhaps seem less like magic as you proceed in your studies. – Marko Riedel Sep 12 '17 at 1:40