Suppose we have $M_q$ balls of color $A_q$ where $1\le q\le N$ and $k$
balls are drawn without replacement and we ask about the probability
that each color is represented. We put $M=\sum_q M_q.$
Using the same concept as in this MSE
link we introduce
the generating function
$$[z^k] \prod_{q=1}^N (1+zA_q)^{M_q}.$$
Now when we set a subset of the $A_q$ to zero we are left with those
terms that are missing these $A_q$ and possibly additional
terms. Therefore we can apply inclusion-exclusion to compute the terms
where none of the $A_q$ is missing. We subtract those where one or
more $A_q$ is missing, then add those where two or more are missing
and so on. Finally we set the remaining $A_q$ to one to obtain the
count. Putting $A=[N]$ we get
$$[z^k] \sum_{S\subseteq A} (-1)^{|S|}
(1+z)^{M-\sum_{q\in S} M_q}
\\ = \sum_{S\subseteq A} (-1)^{|S|}
{M-\sum_{q\in S} M_q\choose k}.$$
We then get for the probability
$${M\choose k}^{-1} \sum_{S\subseteq A} (-1)^{|S|}
{M-\sum_{q\in S} M_q\choose k}.$$
We have the following Maple program to verify these numbers. The total
enumeration routine implements the problem before optimization -- this
was done deliberately in order to adhere to the problem definition.
with(combinat);
CHOOSE :=
proc(Ml, k)
local choice, src, Mq, tp, count, res, cols;
src := []; count := 0;
for tp to nops(Ml) do
Mq := Ml[tp];
src :=
[op(src), seq([tp, q+count], q=1..Mq)];
count := count + Mq;
od;
res := [];
for choice in choose(src, k) do
cols := [seq(choice[q][1], q=1..k)];
res := [op(res), cols];
od;
res;
end;
PROB :=
proc(Ml, k)
option remember;
local choice, count, N, M;
count := 0; N := nops(Ml);
for choice in CHOOSE(Ml, k) do
if nops(convert(choice, `multiset`)) = N then
count := count + 1;
fi;
od;
M := add(q, q in Ml);
count/binomial(M, k);
end;
X :=
proc(Ml, k)
local res, S, N, rest, M;
N := nops(Ml); M := add(q, q in Ml);
res := 0;
for S in powerset([seq(q, q=1..N)]) do
rest := M - add(Ml[q], q in S);
res := res + (-1)^nops(S)*binomial(rest, k);
od;
res/binomial(M,k);
end;
X2 :=
proc(Ml, k)
local res, S, N, rest, M, FF;
FF := (n, q) -> mul(n-p, p=0..q-1);
N := nops(Ml); M := add(q, q in Ml);
res := 0;
for S in powerset([seq(q, q=1..N)]) do
rest := M - add(Ml[q], q in S);
res := res + (-1)^nops(S)*FF(rest, k);
od;
res/FF(M, k);
end;
Following an idea by @Fimpellizieri we can re-write the formula as
$$\frac{1}{M^{\underline{k}}}
\sum_{S\subseteq A} (-1)^{|S|}
\left(M-\sum_{q\in S} M_q\right)^{\underline{k}}.$$