Low probability of a multinomial for expected values

Question

given a population with 4 groups for which their frecuencies are:

A = 0.46
B = 0.075
C = 0.035
D = 0.43

The expected number of individuals for each group in a sample of 10 individuals is:

A = 0.46 * 10 = 4.6 ~ 5
B = 0.075 * 10 = 0.75 ~ 1
C = 0.035 * 10 = 0. 35 ~ 0
D = 0.43 * 10 = 4.3 ~ 4

So, why do I get a low probability of finding 5 individuals from A, 1 from B, 0 from C and 4 from D when extracting a sample of 10 individuals?

dmultinom(x = c(5,1,0,4),prob = c(0.46,0.075,0.035,0.43))
[1] 0.06654184

Answer 1

Having the mean result isn't necessarily the one with the highest probability. Also, the probability of it won't be as high as you'd expect in most cases.

For example, think about a fair coin toss 100 times. The probability of having 50 heads is ${100 \choose 50}0.5^{100}\approx 0.08$. It's not that high as you'd expect. It gets even much smaller as $n$ increases. This is directly related to the number of cases available.

In your problem, there are a total of ${10+3\choose 3}$ different situations.

Answer 2

There are a couple reasons off of the top of my head I can think of why the probability is significantly smaller than your intuition:

1. Playing fast and loose with rounding.

While it seems innocent, and even with a small number of trials the "rounding error" compound very quickly. This is because the multinomial formula,

$$P(X_1 = k_1, X_2 = k_2,\cdots, X_m = k_m) \frac{N!}{k_1!k_2!\cdots k_m!} p_1^{k_1}p_2^{k_2}\cdots p_m^{k_m}$$

involves both factorials and exponentials.

2. Even the most probable multinomial distribution outcome pale in comparison to the number of total combinations.

Even modest probabilities and a low number of trials have a large number of potential outcomes and it becomes increasing difficult to draw computationally or see experimental frequencies of exact outcomes.

Following this problem, and this result might be unexpected, the most probable outcome is $A=5,B=0, C= 0, D= 5$. And can easily be shown by having R run through all the possible combinations:

brute <- c()
d <- 1
ind <- c(0,0,0,0)
maxP <- 0
for (i in 0:10) {
  for (j in 0:10) {
    for (k in 0:10) {
      for (n in 0:10) {
        if (sum(c(i, j, k, n)) == 10) {
          brute[d] <-
            dmultinom(x = c(i, j, k, n),
                      prob = c(0.46, 0.075, 0.035, 0.43))
          if (brute[d] > maxP) {
            ind[1] <- i
            ind[2] <- j
            ind[3] <- k
            ind[4] <- n
            maxP <- brute[d]

          }
          d <- d + 1
        }
      }
    }
  }
}
print(ind)
print(maxP)

And results in:

[1] 5 0 0 5
[2] 0.07630131

Even the most probable outcome happens infrequently.