Expected number of empty bins with clumpy balls?

Question

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!}$

Answer 1

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$$

Answer 2

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*
add(coeftayl(mul(1+w/q, q=1..n-1), w=0, k-1)
    *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);
}