Probability of having a bounded ratio of two types of balls in each of 'S' bins after random partitioning of a fixed number of balls

Question

Let's say I have a bag with $A$ red balls, $B$ blue balls, and a total number of balls $N = A + B$. With uniform probability, and sampling without replacement from the $N$ balls, I fill an integer number of bins, $S$, with exactly $L$ balls each, where $N = S*L$. As such, the bag of all $N$ balls should be empty by the end of the procedure.

What is the probability of having at most $k \leq B$ blue balls in every one of the $S$ bins?

Edit - I would be very much interested in approximate solutions! This problem should essentially come down to something akin to making sure a fixed number of points randomly placed on a matrix are at least some fixed distance apart (causing them to end up in different bins), so surely there have been problems like this tackled in the literature?

Edit 2 - As per Peter Shor's recommendation, I'm interested in examples where $S$ and $L$ are of fixed size, but where $S$ is large, at least $10^3$ to $10^4$ or so, and $L$ is typically significantly smaller than $S$ (around $10$ to $10^2$ or so).

---

Another way of asking this question might be as follows -

Let's say that I have a matrix composed of $S$ blocks of $L$ cells (pick whatever geometry you'd like), where $S > 10^3$ and $L$ is approximately $10^1$ ~ $10^2$ or so. Some fixed number of cells, $A$, store the value '0', while the remaining $B$ cells store the value '1'. To be clear, there are a total of $A + B = S*L = N$ cells in the matrix.

What is the probability that all $S$ blocks contain at most $k$ cells with the storage value '1' (i.e. at most $k$ of the $B$ cells)?

Answer 1

For some parameters, you will not see much of a difference between the exact value and assuming that the balls are independent, or that the bins are independent. The independent approximations are much easier to calculate exactly.

You can estimate the probability using inclusion-exclusion. Whether the computable estimates are useful depends on the parameters.

The probability that there is at least one bin with more than $k$ blue balls equals the sum over the subsets $\beta$ of the bins

$$\sum_\beta (-1)^{|\beta|-1} P(\text{all bins in $\beta$ have more than $k$ blue balls})$$

$$=\sum_i (-1)^{i-1} {S \choose i} P(\text{the first $i$ bins have more than $k$ blue balls})$$

Each partial sum up to an odd $i$ is an upper bound, while each partial sum up to an even $i$ is a lower bound. The terms increase and then decrease roughly like a Poisson distribution, or the terms in the Taylor series for $\exp(x)$.

You can evaluate the probability that the first $i$ bins each have more than than $k$ blue balls as a sum over the possible number of blue balls in each. I don't know whether this can be simplified significantly, but when $(L-k)$ and $i$ are not too large you can calculate these explicitly:

$${S L \choose B}^{-1}\sum_{k \lt b_1,...,b_i \le L} {L \choose b_1} {L \choose b_2} ... {L\choose b_i} {L(S-i) \choose B-(b_1 + ... + b_i)}.$$

The terms become $0$ when $(k+1)i \gt B$, but that might not be useful.

Better is to collect like terms based on the $L-k-1+i \choose i$ possible multisets $\{b_1,...b_i\}:$

$${S L \choose B}^{-1}\sum_\lambda {i \choose c_{k+1}(\lambda)~ c_{k+2}(\lambda) ... c_L(\lambda)} {L \choose k+1}^{c_{k+1}(\lambda)} ... {L \choose L}^{c_{L}(\lambda)} {L(S-i) \choose B-|\lambda|}.$$

Here $c_n(\lambda)$ is the count of parts of size $n$ in $\lambda$, and $\lambda$ has $i$ parts, all of size from $k+1$ to $L$.

For example, suppose $S=1000$, $L=10$, $k=4$, and $B=1000$.

If you approximate the colors of the balls as independent Bernoulli random variables with parameter $1/10$, then the probability none of the $1000$ bins will have $5$ or more blue balls would be $(1-16349374/10^{10})^{1000} =0.19470389363...$.

If you approximate the bins as independent, then the probability that none of the $1000$ bins will have $5$ or more blue balls would be $(1-311442378665580894806964843/191843012418970806869358022430)^{1000}$ $= (1-0.00162342)^{1000} = 0.196962429... $

If I have calculated correctly, inclusion-exclusion gives the following sum:

$1$ $- 1.6234231038105990908$ $+ 1.2919877435018249272$ $- 0.67199565557524499763$ $+ 0.25695414151721206859$ $-0.077037323936817040460$ $+ 0.018861575339778788209$ $-0.0038785183878317956045$ $+0.00068369841392455375101$ $-0.00010494404713497539259$ $+0.000014199866369591258663$ $-1.7106130082560862260*10^{-6}$ $+ 1.8497238308747411999*10^{-7}$ $-...$ $=0.19206027...$

The chances that there are at most $k$ blue balls in each bin is the complement.

Here is some Mathematica code for calculating this:

subsetToMultiset[sub_] := Table[sub[[i]] - i, {i, Length[sub]}]
Clear[lambdas];
lambdas[i_, k_, l_] := lambdas[i, k, l] = 
   k + 1 + Map[subsetToMultiset, Subsets[Table[n, {n, l + i - k - 1}], {i}]]
counts[lambda_, i_, k_, l_, s_, b_] := 
   Product[Binomial[l, j]^Count[lambda, j], {j, k + 1, l}] 
   Multinomial @@ Table[Count[lambda, j], {j, k + 1, l}] 
   Binomial[l (s - i), b - Total[lambda]]
ieTerm[i_, k_, l_, s_,   b_] := 
   (-1)^i Binomial[s, i] 
   Sum[counts[lambdas[i, k, l][[ind]], i, k, l, s, b], 
       {ind, Length[lambdas[i, k, l]]}]
   /Binomial[s l, b]
Table[{i, N[ieTerm[i, 4, 10, 1000, 1000], 20]}, {i, 0, 12}]
Total[%][[2]]

This takes $7$ seconds on my computer.