# Combinatory probability, group composed by different balls of different colors

Let's say you have a group of $$M$$ balls of different colors in a box. For example, 20 balls are red, 15 are blue, 10 are green, 5 are grey, 5 are yellow and 5 violet, for a total of $$M=60$$ balls. You pick $$1 \leqslant n \leqslant M$$ of them without replacement. The order of the colors does not count, so for instance, if $$n=2$$ and you pick red then gray, it is the same as picking gray then red.

How do I calculate the probability of all the possible outcomes for $$n$$ elements? Is there a general formula for this problem? In particular, if $$M$$ and $$n$$ are large then the number of possible combinations is large, so how can I find the most probable combinations?

• "Probability of all the possible outcomes" seems either too inclusive or too vague. Specifically, if 4 balls drawn without replacement, then probability of getting exactly 2 red, 1 green and 1 yellow is $\frac{{20\choose 2}{10\choose 1}{5\choose 1}{25\choose 0}}{{60\choose 4}}.$ May 18, 2020 at 21:22
• See Wikipedia on 'multinomial distribution' for more. May 18, 2020 at 21:30

Yes, there is a general formula. Consider an urn containing $$M$$ balls, where $$M_1$$ balls have the color $$c_1$$, $$M_2$$ balls have the color $$c_2$$,..., $$M_r$$ balls have the color $$c_r$$, and $$M_1+\cdots+M_r=M$$. If you draw a sample of size $$n without replacement the sample space is $$\Omega=\{\omega\colon \omega=(a_1,\dots,a_n), a_i\ne a_j, i\ne j\}$$and $$|\Omega|=(M)_n=\frac{M!}{(M-n)!}$$. Consider an event $$B_{n_1,\dots,n_r}$$ in which $$n_1$$ balls have color $$c_1$$,..., $$n_r$$ balls have color $$c_r$$, where $$n_1+\cdots+n_r=n$$. The $$c_1$$ balls can get $$C_n^{n_1}=\binom{n}{n_1}=\frac{n!}{(n-n_1)!n_1}$$ sets of $$n_1$$ indexes in $$(a_1,\dots,a_n)$$, the $$c_2$$ balls can get $$C_{n-n1}^{n_2}$$ sets of $$n_2$$ indexes, etc., and you can choose $$(M_i)_{n_i}=\frac{M_i!}{(M_i-n_i)}$$ balls that have color $$c_i$$. The general number of events is:\begin{align*}|B|&=\frac{n!}{(n-n_1)!n_1!}\frac{(n-n_1)!}{(n-n_1-n_2)!n_2!}\cdots\frac{(n-n_1-\dots-n_{r-1})!}{(n-n_1-\cdots-n_r)!n_r!}\prod_{i=1}^r (M_i)_{n_i}\\&=\frac{n!}{(n-n_1)!n_1!}\frac{(n-n_1)!}{(n-n_1-n_2)!n_2!}\cdots\frac{n_r!}{0!n_r!}\prod_{i=1}^r (M_i)_{n_i}\\&=\frac{n!}{n_1!\cdots n_r!}\frac{M_1!}{(M_1-n_1)!}\cdots\frac{M_r!}{(M_r-n_r)!}\\&=n!C_{M_1}^{n_1}\cdots C_{M_r}^{n_r}\end{align*}and$$P(B)=\frac{|B|}{|\Omega|}=\frac{n!C_{M_1}^{n_1}\cdots C_{M_r}^{n_r}}{(M)_n}=\frac{C_{M_1}^{n_1}\cdots C_{M_r}^{n_r}}{C_M^n}$$ The set of probabilities $$\{P(B_{n_1,\dots,n_r})\}$$ is called the multivariate hypergeometric distribution. See Shiryaev, Probability, 1996, or Probability 1, 2016, Chapter 1, §2.

If you are using R to compute $$P(B)$$, you can install the extraDistr package:


> K <- 10 # sample size
> x <- subset(expand.grid(red=0:20, blue=0:15, green=0:10, gray=0:5, yellow=0:5, violet=0:5), red+blue+green+gray+yellow+violet==K)
> dim(x)
[1] 2625    6
red blue green gray yellow violet
11   10    0     0    0      0      0
31    9    1     0    0      0      0
51    8    2     0    0      0      0
71    7    3     0    0      0      0
91    6    4     0    0      0      0
111   5    5     0    0      0      0
> tail(x)
red blue green gray yellow violet
739201   0    0     0    2      3      5
753986   1    0     0    0      4      5
754006   0    1     0    0      4      5
754321   0    0     1    0      4      5
757681   0    0     0    1      4      5
776161   0    0     0    0      5      5
> p <- dmvhyper(x, n=c(20,15,10,5,5,5), k=K)
> max(p)
[1] 0.008930581
> x[which.max(p),]
red blue green gray yellow violet
159646   3    2     2    1      1      1
> dmvhyper(x[which.max(p),], n=c(20,15,10,5,5,5), k=K)
[1] 0.008930581
> choose(20,3)*choose(15,2)*choose(10,2)*5*5*5/choose(60,K)
[1] 0.008930581

• Usually, comments are not for this but your answer is absolutely beautiful and impressive. I was able to get a result but yours is amazing and efficient. May 20, 2020 at 13:32
• You are welcome. Happy to be useful! May 20, 2020 at 13:35

After the great input of @BruceET I was able to come with a decent approximation of the question. Unfortunately is uses a brute force approach. The idea is the following. Simulate the process of extracting to see which outcome is the most probable (or which are the "n" most probable results). In the following case, I extracted 12 elements and I did a statistics using 10 million extractions.

library(parallel)
library(plyr)

# name of the elements of the set
acolors=c("B","G", "O", "R", "Y", "Gr")

# generating a real set with a certain composition
alist=c(rep("B",20), rep("G", 15), rep("O", 10), rep("R",5), rep("Y",5), rep("Gr",5) )

# number of extraction to simulate
pulls=10000000

# parallel version of the extraction
all_res=mclapply(1:pulls, function(x, alist, acolors){
ares=NULL
asamp=list(table(sample(alist, 12)))
for(ac in acolors){
if(is.na(asamp[[1]][ac])){
ares=c(ares,0)
}
else{
ares=c(ares,asamp[[1]][ac])
}
}
return(ares)
}, alist=alist, acolors=acolors, mc.cores=8)

# tdata store the result of each extraction.  Each column has a given name
# corresponding to "acolors"
# a line can look like 4,2,3,1,1,1, that mean 4 from color B, 2 from color G and so on...
tdata=as.data.frame(do.call(rbind, all_res))

colnames(tdata)=acolors

## now we do the statistics of the result, we count how many times a given line is duplicated
stat_res=as.data.frame(ddply(tdata,.(B, G, O, R, Y, Gr),nrow))

## we sort the data frame from the most probable to the least probable
stat_res=stat_res[order(stat_res$V1, decreasing = TRUE ),] ## we calculate the frequency of each line stat_res$$frequency=stat_res$$V1/sum(stat_res$V1)


With the results is then possible to use the binomial formula given by @BruceET and calculate the probability for the most frequent results.