# Multinomial distribution approach as negative binomial?

I was asked to solve a problem, following the multinomial distribution (if you think is the correct approach, which I actually think it is). I paste the whole problem just in case I am approaching, but it's the last item which really bothers me.

Formulation The ABO blood group system allows classifying individuals into four mutually exclusive blood groups: A, B, AB, and O. In a population, the frequencies of these groups are 0.46 (A), 0.075 (B), 0.035 (AB), and 0.43 (O).

• Suppose the blood of 8 individuals is analyzed. How many individuals do we expect to find in each group?

• In the same sample as the previous case, what would be the probability of obtaining two individuals from each group?

• And how many individuals it would be required to analyze to obtain at least two individuals per each group? (Meaning 2 or more in each group)

$$P(X_1 = x_1, X_2 = x_2, ..., X_k = x_k) = \frac{n!}{x_1!x_2!...x_k!} \cdot p_1^{x_1} \cdot p_2^{x_2} \cdot \ldots \cdot p_k^{x_k}$$

As for me the two first questions follow like this:

# Prob for each group
prob <- c(0.46, 0.075, 0.035, 0.43)

# Sample size
n <- 8

# Applying multinomial distribution
prob_sol <- rmultinom(1, n, prob)

# Probabilities of obtaining two per group

prob_two_group <- dmultinom(c(2, 2, 2, 2), n, prob)



But for the last part, when trying to approach as reaching a specific minimum number of successes as negative binomial.

I have seen a function dnegmn from the package MGLM, which is called multinomial negative. I don't know for sure if this is the correct function, and how the arguments operate

• I don't get your comment. The probabilities exposed before are not sufficient to reach a solution? Commented Mar 24 at 16:21
• The last question doesn't need a probability calculation to be answered. For example, if there are 1000 people in the population, then 460 of them are A, 75 are B, 35 are AB, and 430 are O. Therefore, you need to sample at least 1000 - 33 = 967 to guarantee two per group. Commented Mar 24 at 16:21
• But if your population is infinite, there's no absolute guarantee at any sample size that you will have at least two per group. Commented Mar 24 at 16:22
• Because there's no certainty you will succeed in obtaining two individuals from each group, you need to formulate your last question as a probability. For example, you could find the smallest sample needed to include at least two from each group with 95% (or whatever) probability. I therefore suspect that either your problem is not well formulated or you haven't accurately described it here. Could you please verify your statement of the problem is correct?
– whuber
Commented Mar 24 at 16:45
• Sorry for the confussion. The question is not well formulated. My guess, the question, as whuber proposed, tries to ascertain the smallest sample necessary, with a predetermined level of confidence. Can you help me with that? Commented Mar 24 at 18:31

Assuming the adapted case (how many individuals for a given probability), then this sounds a lot like an occupancy problem or coupon collectors problem.

The general solution is not easy to compute, but in this small case an exact distribution can be computed

• Exact computation: you could compute a Markov chain for the 80+1 states where the occupancy of at least one is below 2, and the end state where all are at least 2.

• Estimate: A simple way would be to estimate by considering only the rare groups B and AB which are the bottleneck and approximate them with a Poisson distribution or beta-binomial.

The two methods give very similar results

### create a transition matrix
T = matrix(rep(0,81*81),81)
for (A in 0:2) {
for (B in 0:2) {
for (AB in 0:2) {
for (U in 0:2) {
row = 1+A+B*3+AB*9+U*27
T[row,row+(A<2)]    = T[row,row+(A<2)] + 0.46
T[row,row+3*(B<2)]  = T[row,row+3*(B<2)] + 0.075
T[row,row+9*(AB<2)] = T[row,row+9*(AB<2)] + 0.035
T[row,row+27*(U<2)] = T[row,row+27*(U<2)] + 0.43
}}}}

## begin state is position 0
state = t( c(1,rep(0,80)) )
### variable to keep track of end state
p_end = c(0)

### compute 150 transitions
t = 150
for (i in 1:t) {
state = state %*% T
p_end = c(p_end,state[81])
}

ts = 0:t

### approximate with independent binomial distributions
p_approx = (1-pbinom(1,ts,0.075))*(1-pbinom(1,ts,0.035))

### plotting
plot(ts,p_approx, main = "compare exact (points) with approximation (line)", xlab = "number of individuals", ylab = "probability of occupancy for each group 2 or more", type = "l", col = 1)

points(ts[ts %% 5 == 0],p_end[ts %% 5 == 0], pch = 21, bg = 0, cex = 0.8)


Other approach

• For larger problems (more states and larger number of occupancy) the distribution may be approximated with a Gumbel distribution, which relates to the distribution of the minimum of several variables (and the problem is equivalent to the minimum of the variables being above some value).

Or for intermediate problems the minimum distribution can be approximated by assuming independence and multiplying the survival functions (one minus the cumulative distribution function). In the code above this is simplified by regarding only the blood groups AB and B.

• Thank you for the answer. It's a lot to take in, and I haven't heard most of the theory background. I have several doubts, why this = 1+A+B*3+AB*9+U*27 , and the final result is extrapolated from the graph, should it be close to 150 (p = 0.95). It was mentioned before "there's no certainty you will succeed in obtaining two individuals from each group" Commented Mar 25 at 14:01
• @JavierHernando The graph is just explaining how it works. There is no need to extrapolate from the graph. The code can be extended to compute different values, but you haven't explained what you are looking for exactly. Commented Mar 25 at 16:23
• The idea behind 1+A+B*3+AB*9+U*27 is to map decimal numbers to ternary numbers. Here we only regard the number of cases in the A, B AB and O groups, when they have values $0$, $1$ or $\geq 2$. Commented Mar 25 at 16:26