Multinomial distribution: probability that one outcome is greater than another

Question

Consider a multinomial distribution with three outcomes. Let $x_i$ denote the number of occurences of the $i^{th}$ outcome, and the $i^{th}$ outcome occurs with probability $p_i$, $i=1,2,3$. Let $n$ be the number of total trials. Then we have

$$(X_1,X_2,X_3)\sim Multi(n;p_1,p_2,1-p_1-p_2).$$

The questions are:

1) How can I calculate the probability of $x_1\geq x_2$ for a given $n$? Is there a concise closed form representation for $P[X_1\geq X_2|n]?$

2) Do we have some asymptotic property of this probability? Intuitively, if $p_1>p_2$, the probability should converge to 1. On the other hand, if $p_1<p_2$, it should converge to zero.. Can we see this from answering the first question?

Answer 1

There is a closed form for this probability, but it is not particularly concise. To obtain the formula, start with the probability of the event of interest, conditional on the event $X_1+X_2=r$. It can easily be shown that:

$$X_1 | X_1 + X_2 = r, n, \mathbf{p} \sim \text{Bin} \Big( r, \frac{p_1}{p_1+p_2} \Big).$$

Under the stated conditioning event, the event $X_1 \geqslant X_2$ is equivalent to the event $X_1 \geqslant r/2$. We therefore have the conditional probability:

$$\begin{equation} \begin{aligned} \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) &= \mathbb{P}(X_1 \geqslant r-X_1 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= \mathbb{P}(X_1 \geqslant r/2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= \mathbb{P}(X_1 \geqslant \lceil r/2 \rceil | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= \frac{1}{(p_1+p_2)^r} \sum_{x = \lceil r/2 \rceil}^r {r \choose x} \cdot p_1^x \cdot p_2^{r-x}. \\[6pt] \end{aligned} \end{equation}$$

In practice, the upper tail of the binomial distribution is usually computed using the regularised incomplete beta function, rather than by summing the mass function. In this latter form, we have:

$$\begin{equation} \begin{aligned} \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) &= 1-\mathbb{P}(X_1 \leqslant \lceil r/2 \rceil - 1 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= 1-I_{p_2/(p_1+p_2)}(r - \lceil r/2 \rceil + 1, \lceil r/2 \rceil). \\[6pt] \end{aligned} \end{equation}$$

Applying the law of total probability gives the marginal probability of interest:

$$\begin{equation} \begin{aligned} \mathbb{P}(X_1 \geqslant X_2 | n, \mathbf{p}) &= \sum_{r=0}^n \mathbb{P}(X_1 + X_2 = r | n) \cdot \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= \sum_{r=0}^n \text{Bin}(r|n, p_1+p_2) \cdot \mathbb{P}(X_1 \geqslant X_2 | X_1 + X_2 = r, n, \mathbf{p}) \\[6pt] &= \sum_{r=0}^n {n \choose r} \cdot (1-p_1+p_2)^{n-r} \sum_{x = \lceil r/2 \rceil}^r {r \choose x} \cdot p_1^x \cdot p_2^{r-x} \\[6pt] &= \sum_{r=0}^n {n \choose r} \cdot (p_1+p_2)^{r} \cdot (1-p_1+p_2)^{n-r} \\[6pt] &\quad \quad \quad \times (1-I_{p_2/(p_1+p_2)}(r - \lceil r/2 \rceil + 1, \lceil r/2 \rceil)). \\[6pt] \end{aligned} \end{equation}$$

Either of the last two lines give a closed form solution to the problem. (The first is a sum over relevant values of the mass function of the binomial. The second uses the regularised incomplete beta function for the distribution. The latter is closer to what you would use when implementing the function in statistical computing.) Either expression is somewhat cumbersome for large $n$.

---

Computing the probability: We can construct a function in R to compute this probability using the standard statistical distribution functions. (Since this computation involves adding very small probabilities we work in log-space for the computation, and then convert back to a standard probability at the end.)

PROB <- function(n, p, log = FALSE) {

#Normalise probability vector
p <- p/sum(p);

#Computeprobability
LOGTERMS <- rep(0, n+1);
LOGTERMS[1] <- dbinom(0, size = n, prob = p[1]+p[2], log = TRUE);
CLOGPROB <- rep(0, n+1);
for (r in 1:n) {
  CLOGPROB[r+1] <- pbinom(ceiling(r/2)-1, size = r, prob = p[1]/(p[1]+p[2]), 
                          lower.tail = FALSE, log.p = TRUE);
  LOGTERMS[r+1] <- dbinom(r, size = n, prob = p[1]+p[2], log = TRUE) + CLOGPROB[r+1]; }
LOGPROB <- matrixStats::logSumExp(LOGTERMS);

if (log) { LOGPROB } else { exp(LOGPROB) } }

This function is sufficiently efficient to compute the probability of interest up to quite large values of $n$. Basic tests show that it will compute for $n = 10^7$ in a few seconds. Below are some examples of use of this function for various values of $n$. As can be seen, even when the probability gets too small to see on the standard scale, we can still compute it on the log-probability scale.

p <- c(0.20, 0.25, 0.55);

PROB(1, p);
[1] 0.75

PROB(2, p);
[1] 0.6625

PROB(10, p);
[1] 0.4987803

PROB(100, p);
[1] 0.2503806

PROB(10^5, p);
[1] 1.174474e-123

PROB(10^6, p);
[1] 0

PROB(10^6, p, log = TRUE);
[1] -2795.471

PROB(10^7, p, log = TRUE);
[1] -27909.27

This function can be used to obtain the probability of interest up to large values of $n$. For values that are so large that they are computationally infeasible with this algorithm, we can then fall back on the normal approximation to the multinomial distribution.

Answer 2

Let $Z_i$ be independent draws from a categorical variable take takes values 1,2 and 3 with probability $p_1,p_2$ and $1-p_1-p_2$. Then, $X_k = \sum_{i=1}^n 1\{Z_i = k\}$ has the desired multinomial distribution. Let $Y_i = \{Z_i = 1\} - 1\{Z_i = 2\}$ and let $W_n = \sum_{i=1}^n Y_i$. The desired probability is $P(W_n \ge 0)$.

Note that $Y_i$ takes values $1, -1$ and zero with probabilities $p_1, p_2$ and $1-p_1-p_2$. $W_n$ is a nonsymmetric lazy random walk. I am going to only consider the asymptotics. We have $E[Y_i] = p_1 - p_2$ and $var(Y_i) = p_1 + p_2 - (p_1-p_2)^2 =: v^2$. By CLT, $$ U_n := \frac{W_n - n(p_1-p_2)}{\sqrt{n} v} \to N(0,1) $$ Thus, $$P(W_n \ge 0) = P( U_n \ge -\sqrt{n}(p_1-p_2)/v) $$ which is going to go to $1$ or $0$ depneding on whether $p_1 > p_2$ or $p_1 < p_2$, respectively. A good approximation for large $n$ is $$ P(W_n \ge 0) \approx \Phi(\sqrt{n}(p_1-p_2)/v) $$ where $\Phi$ is the CDF of the normal distribution.