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?
 A: 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.
