MLE of Negative Binomial distributions of different sizes

Question

There are two teams that are competing in a series of matches. These matches are a best-of-$x$ format, so after either team wins $\lfloor{\dfrac{x}{2}}\rfloor+1$ the series is over. This is, in essence, a negative binomial distribution.

These two teams play a few different matches, with different $x$ values, and I want to know what the MLE is of $p$, where $p$ is the probability that Team $1$ defeats Team $2$.

Here is what I attempted:

Our likelyhood function is as follows: $$ L(k_i, x_i; p) = \prod_{i=1}^n c_ip^{k_i}(1-p_i)^{x_i-k_i} $$ where $x_i$ notes that match $i$ was a Best-of-$x$ and $k_i$ is the number of matches won by Team $1$.

If Team $1$ won match $i$, then $c_i = \binom{x-1}{k-1}$, else $c_i = \binom{x-1}{k}$. Take the $\log$ of the likelyhood, we find:

$$ \log{L} = \sum_{i=1}^n c_i+\log(p)\sum_{i=1}^nk_i + \log(1-p)\sum_{i=1}^n(x_i-k_i) $$

differentiating this with respect to $p$, and setting that to $0$ to find the maximum, we get that:

$$ 0 = \dfrac{\sum_{i=1}^nk_i}{p} - \dfrac{\sum_{i=1}^n(x_i-k_i)}{1-p} $$ $$ 0 = (1-p)\sum_{i=1}^nk_i - p\sum_{i=1}^n(x_i-k_i) $$ $$ \sum_{i=1}^nk_i = p\left(\sum_{i=1}^nk_i + \sum_{i=1}^n(x_i-k_i) \right) = p\sum_{i=1}^nx_i $$

which leads us to the final result

$$ p = \dfrac{\sum_{i=1}^nk_i}{\sum_{i=1}^nx_i} = \dfrac{wins}{games} $$

From here I was a little skeptical of this result, so I simulated sets of matches for various probabilities of victory between teams to find the numerical MLE.

For example, I did the scenario in which there is are $2$ Best-of-$1$s, and one Best-of-$5$. Team $1$ loses the best of $1$s, but wins the Best-of-$5$ $3$-$1$. Numerically, I found the MLE to be $p\approx .65$, but using the MLE I just found I would get that

$$ p = \dfrac{0+0+3}{1+1+4} = \dfrac{1}{2} $$

Answer 1

The likelihood that we want to maximize for your example is

$$(1-p)\cdot (1-p)\cdot 3 \cdot p^3\cdot (1-p)=3[p(1-p)]^3$$

It indeed has maximum value at $p=0.5$.

Perhaps you might like to check your simulation. If the success probability is more than $0.5$, we should expect the team to win more frequently.

Answer 2

The terminology is confusing and not sufficiently precisely stated, which I would like to address by establishing some definitions and notation.

Let us call a game a single instance of competition between the two teams, which I will denote by $A$ and $B$, in which the outcome of the game is either $1$ if team $A$ wins, and $0$ if team $A$ loses. In other words, the outcome we are counting in a single game is the number of points won by team $A$.

A match comprises a sequence of games, of which there is a predetermined number of games played called the match length, and the outcome of the match is in favor of the team with more winning games in the match. The outcome of the match is expressed by the total points won by $A$.

A series comprises multiple matches, in which the match lengths may vary from match to match. The winner of the series is determined in a similar fashion to match decisions; i.e., whichever team wins more matches wins the series.

Now let us define some notation. Let $m$ represent the number of matches played, and for each $i = 1, 2, \ldots, m$, let $n_i$ represent the match length of match $i$. Moreover, let $X_i$ represent the random number of games/points won by team $A$ in match $i$.

Next, suppose that the outcomes of individual games are independent and identically distributed, in which the probability of team $A$ winning a game/point is $p$, and the probability of team $B$ winning a game/point is $q = 1-p$. Then the number of wins $X_i$ for team $A$ in match $i$ is a binomial random variable; $$X_i \sim \operatorname{Binomial}(n_i, p), \quad i = 1, 2, \ldots, m.$$ Furthermore, each $X_i$ is assumed to be independent, although not identically distributed since $n_i$ may not be equal for each $i$.

The goal is to estimate $p$ based on the observed data $$\boldsymbol x = (x_1, x_2, \ldots, x_n)$$ and match lengths $\boldsymbol n = (n_1, n_2, \ldots, n_m)$, where each $x_i$ is the outcome of match $i$. In such a case, note that the joint distribution of $\boldsymbol x$ is $$f_{\boldsymbol X}(\boldsymbol x \mid p, \boldsymbol n) = \prod_{i=1}^m \binom{n_i}{x_i} p^{x_i} (1-p)^{n_i - x_i} = \left(\prod_{i=1}^m \binom{n_i}{x_i}\right) p^{\sum x_i} (1-p)^{\sum n_i - \sum x_i} = h(\boldsymbol x) g(T(\boldsymbol x) \mid p),$$ where $$h(\boldsymbol x) = \prod_{i=1}^m \binom{n_i}{x_i}, \quad g(T \mid p) = p^T (1-p)^{\sum n_i - T}, \quad T(\boldsymbol x) = \sum x_i,$$ by the Factorization Theorem. Consequently $$T(\boldsymbol x) = \sum x_i$$ is a sufficient statistic for $p$, and it is a straightforward exercise to show that $$T \sim \operatorname{Binomial}(\Sigma n_i, p).$$ Consequently, the MLE of $p$ is expressible as a function of $T$ and is the same functional form as the MLE for the binomial distribution.

What the above tells us is that if you know the match scores, then it doesn't matter what the individual match lengths are in so far as estimating $p$. This makes sense because each game is an independent and identically distributed $\operatorname{Bernoulli}(p)$ random variable.

The situation is more complex, however, when instead of knowing the match scores, you only know the match result--i.e., whether team $A$ or team $B$ won match $i$ for each $i = 1, 2, \ldots, m$. (However, as before, you still know the match lengths, you just don't know how many points were won in each match by each team.) In such a case, the likelihood function is less informative, and we must define another random variable to characterize this situation.

Specifically, as described above, team $A$ wins match $i$ if $X_i \ge \lfloor n_i/2 \rfloor + 1$. This of course is problematic when $n_i$ is even, since an even match length admits a tied outcome for that match with probability $\binom{n_i}{n_i/2} (p(1-p))^{n_i/2}$, and it is not clear in the question statement how to handle such a case. So for the time being, let us assume that $n_i$ is odd for all $i$. We then define the random variable $$W_i = \begin{cases}1, & X_i > n_i/2 \newline 0, & x_i < n_i/2. \end{cases}$$ Then $$W_i \sim \operatorname{Bernoulli}(\pi_i), \quad \pi_i = \Pr[X_i > n_i/2] = \sum_{x = \lceil n_i/2 \rceil}^{n_i} \binom{n_i}{x} p^x (1-p)^{n_i-x}.$$ Then our likelihood must be constructed from the sample $\boldsymbol w = (w_1, w_2, \ldots, w_m)$ rather than $\boldsymbol x$, and we have $$\mathcal L(p \mid \boldsymbol w, \boldsymbol n) = \prod_{i=1}^m \pi_i^{w_i} (1-\pi_i)^{1-w_i}.$$ The log-likelihood is $$\ell(p \mid \boldsymbol w, \boldsymbol n) = \sum_{i=1}^m w_i \log \pi_i + (1-w_i) \log (1 - \pi_i)$$ and we seek a critical point satisfying $$0 = \frac{d\ell}{dp} = \sum_{i=1}^m \left( \frac{w_i}{\pi_i} - \frac{1-w_i}{1 - \pi_i} \right)\frac{d\pi_i}{dp}.$$ As you can already see, this is quickly turning intractable. We could continue to calculate the derivative but it seems unlikely that we will be able to analytically solve for a critical point for a general sample.

But maybe we can we do it for a simple case? Suppose $m = 3$, $\boldsymbol n = (5,3,7)$, and $\boldsymbol w = (1,1,0)$. Then we can calculate $$\begin{align*} \pi_1 &= p^3 \left(6 p^2-15 p+10\right) \newline \pi_2 &= (3-2 p) p^2 \newline \pi_3 &= p^4 \left(-20 p^3+70 p^2-84 p+35\right). \end{align*}$$ Our likelihood is $$\begin{align*} \mathcal L(p \mid \boldsymbol w, \boldsymbol n) &= \pi_1 \pi_2 (1-\pi_3) \newline &= -240 p^{15}+1800 p^{14}-5668 p^{13}+9602 p^{12}-9240 p^{11}+4795 p^{10} \newline &\quad -1050 p^9-12 p^8+48 p^7-65 p^6+30 p^5. \end{align*}$$ This confirms our earlier concern about tractability, since the result is a high-degree polynomial in $p$ for which it is not possible to obtain an analytic solution for the maximum. After factoring out the trivial critical points $p \in \{0, 1\}$, we can use a numeric approach (e.g. Newton's method) to locate the unique critical point in $(0,1)$: $$\hat p \approx 0.55616759186866770384.$$ So we can do the calculation in practice, provided that the number of games in a match is not too large nor the number of matches too large, in which case the polynomial $\mathcal L$ has a very high degree.

How would we confirm the above result by simulation? One way is to simulate numerous samples $\boldsymbol x$ each with some random $p$ drawn from a uniform distribution on $(0,1)$. Then select only those samples corresponding to the observed outcome $\boldsymbol w$. Then compute the sample mean of $p$ and compare with the MLE. The problem with this approach is that you can waste a lot of computational effort on simulating outcomes $\boldsymbol x$ that do not give rise to the desired $\boldsymbol w$. I invite the interested reader to devise a more efficient simulation.

---

The above likelihood can be adjusted to consider tied outcomes, in which case $w_i$ is no longer Bernoulli but categorical. Also, the model becomes completely different if we suppose that the match length is not predetermined but is stopped when a certain margin of wins is attained (e.g. a "win-by-2" rule). But the only interpretation of the original question that makes sense to me is to consider match lengths fixed, since it was not specified that a minimum number of games needed to be played in a match, nor a minimum number of matches in a series.