# Expected Balls in Bins of Unequal Capacities

## Question:

There are $n$ bins each with different capacities, $c_1, c_2, .., c_n$.

There are $m$ balls, where $m = p \sum_i c_i$, for $0<p<1$.

Suppose the capacities and number of balls are large so this problem emulates one of a continuous nature.

Balls are sequentially and randomly assigned to bins with uniform probability amongst those that have remaining capacity.

Suppose the number balls allocated to bin $i$ is defined by random variable $B_i$.

Is there a closed form function (of parameters $n, p, c_1, .., c_n$) for $E[B_i]$?

### Some numeric analysis

I have implemented the following Monte Carlo Analysis in python:

def bin_expectation(c, p, samples=1000):
"""
Calculate the expectation of bin allocation based on the uneven capacities of each bin.
Ball distribution is performed sequentially and uniformly over all (unfilled) bins.

Args:
c (array): array of bin capacities
p (float): proportion of total capacity to be distributed, i.e. sum(c)*p = total number of balls.
samples (int): the number of monte carlo trials to generate statistics from

Returns:
Array: the mean of number of allocated balls in each bin.

Method:
i) In each iteration divide the remaining ball allocation into n equal piles.
ii) Generate a uniform RV for each pile which describes the pile proportion allocated to the bin.
iii) If the bin becomes over allocated reset to its max capacity.
iV) repeat until all balls allocated.
"""
t_balls = c.sum()*p               # total number of balls
n = len(c)                        # number of bins
c_arr = np.tile(c, (samples, 1))  # create repeated c arr for each MC sample

a = np.zeros(shape=(samples, n))  # allocation array for each sample

for i in range(40):               # 20-30 iterations is generally enough, 40 is conservative
rv = stats.uniform.rvs(size=n*samples, loc=0, scale=1/n)
rv = np.reshape(rv, newshape=(samples, n))  # generate RVs for allocation
r_balls = t_balls - a.sum(axis=1)  # the remaining number of balls for each sample
a += np.einsum('ij,i->ij', rv, r_balls)  # adjust the allocation
a = np.where(a > c_arr, c_arr, a)  # correct for over allocated bins

return a.mean(axis=0)


If I run the above function on the parameters:

c = np.array([100, 200, 300, 400, 500])
p = 0.5
Exp_a = bin_expectation(c, p, samples=2000)


Then the chart below plots the expected allocation of each bin for different $p$. The graphs look like there might be some nice analytical formula but then it also might be quite horribly algebraically cumbersome?

For any fixed parameters $c_1,\ldots, c_n$ and $m$ you can derive a formula by calculating coefficients of related combinatorial generating functions.

I would not expect this problem to admit a simple closed solution however for all parameter choices (though I may be wrong). Doing asymptotic analysis of the generating functions may then yield some asymptotic statement about the continuous limit (aside: I think this is what you are interested in, as I saw a previous question about the continuous parameter which has now been deleted...)!

First of all I note that you are describing the uniform distribution on all admissable assignments of $m$ balls to the $n$ boxes: where I use admissable to mean that there are at most $c_j$ balls in box $1 \leq j \leq n$.

That is any given configuration has probability $1/N_m$ with $N_m$ the total number of admissable configurations of $m$ balls.

The generating function $G$ of this sequence is the polynomial

$$G(x) = \sum_{m=0}^\infty N_m x^m,$$

For your problem $G(x)$ can be shown to be given by

$$G(x) = \prod_{j=1}^n (1 + x+ \cdots +x^{c_j} )$$

A brief explanation of this formula is that each of the $n$ terms in the product corresponds to a box, and for the $j$-th entry $(1 + x + \cdots + x^{c_j})$ then the term in $x^k$ corersponds to having $k$ balls in box $j$. When you expand the product, the total number of combinations of $m$ balls is the coefficient corresponding to $x^m$.

As an example if you start the expansion of the above you see that according to the above $N_0 = 1$ (which is as you'd expect, there is exactly one way to not distribute any balls), $N_1 = n$ (there are $n$ boxes, and you can put one ball in any of them), etc. Note also that if $m > c_1 + \cdots + c_n$ then there is $x^m$ does not appear in $G(x)$, i.e. has a coefficient of $0$, and so $N_m = 0$ as expected.

We can in general get at any coefficient by taking derivatives, that is $$N_m = \frac{1}{m!}\left( \frac{d^m}{dx^m}G\right)(0).$$

The Expected Number of Balls in Box $j$

We will find $B_j$, the expected number of balls in a given box $j$, by conditioning. That is we note

\begin{align*} \mathbf E_m[ B_j] & = \sum_{k=0}^{c_j} k \mathbf P[B_j = k] \\ & = \frac{1}{N_m} \sum_{k=0}^{c_j} k N_m^{(k,j)}, \end{align*} where I use $N_m^{(k,j)}$ to denote the number of admissible configurations which have exactly $k$ balls in box $j$. We respectively denote the generating function of this sequence by $G^{(k,j)}$ and note that this is given by the formula

\begin{align*} G^{(k,j)}(x) & = x^k \prod_{i \neq j} (1 + x + \cdots + x^{c_i} ) \\ &= \frac{x^k}{(1 + x + \cdots + x^{c_j} )} \prod_{i=1}^n (1 + x + \cdots + x^{c_i} ) \\ & = \frac{x^k}{(1 + x + \cdots + x^{c_j} )} G(x) \end{align*} As before, we have that

$$N_m^{(j,k)} = \frac{1}{m!}\left( \frac{d^m}{dx^m} G^{(k,j)}\right)(0),$$ but computing the expectation is then facilitated by the fact that we have

\begin{align*} \sum_{k=1}^{c_j} k N_m^{(k,j)} & = \frac1{m!} \sum_{k=0}^{c_j} k \left( \frac{d^m}{dx^m} G^{(k,j)}\right)(0) \\ & =\frac1{m!} \left( \frac{d^m}{dx^m} \sum_{k=0}^{c_j} k \, G^{(k,j)}\right)(0). \end{align*} On using the identity $$\sum_{k=0}^{c_j} k \,x^k = \frac{x}{(x-1)^2} \big( c_j x^{c_j+1} - (c_j + 1) x^{c_j} + 1\big)$$ we can define

\begin{align*} G^{(j)}(x) & \colon= \sum_{k=0}^{c_j} k \, G^{(k,j)}(x)\\ & = \frac{x}{(x-1)^2} \frac{\big( c_j x^{c_j+1} - (c_j + 1) x^{c_j} + 1\big)}{\big(1 + x + \cdots + x^{c_j} \big)} G(x), \end{align*}

and we have

$$\mathbf E_m[B_i] = \frac{ \left( \frac{d^m}{dx^m} G^{(j)} \right)(0)}{ \left( \frac{d^m}{dx^m} G \right)(0)} \propto \left( \frac{d^m}{dx^m} G^{(j)} \right)(0)$$

So we have derived a formula for which an expression for the expectation can be computed; how easy this is to do, however, is not all that clear. The complexity of the formulae above hint that finding a closed formula is out of scope.

However, finding a statement that holds in a limit may however be possible (note that defining an appropriate limit is a first essential step,since for any fixed $c_j$ the problem is only defined for $m \leq 1 + \cdots + c_n$); a comprehensive reference for this is Flajolet and Sedgewick's Analytic Combinatorics which is available freely online.

• Thanks, Owen, thats really clear and helpful in my situation - this problem does serve a practical use in another area. I appreciate the time. I did delete the other question since after a bit of research this seemed to be a better reformulation and I had also done a bit of numerical analysis that might have been considered useful for the problem. – Attack68 Mar 7 '18 at 7:09
• No problem at all; I was originally going to respond to the original by suggesting an approximation with the discrete case! – owen88 Mar 7 '18 at 12:55