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.

in each binof blue balls to red balls is bounded between them? My apologies, but the wording of your question seems ambiguous to me, or at least unclear. $\endgroup$3more comments