---
title: Negative Binomial Order Statistics
url: https://www.emergentmind.com/topics/negative-binomial-order-statistics
type: topic
---

# Negative Binomial Order Statistics

Searching arXiv for the cited papers and closely related terminology.
Negative binomial order statistics are the order statistics obtained from independent and identically distributed negative binomial random variables, most commonly the \(r\)-th smallest value \(X_{(r)}\) among \(D\) iid draws from a negative binomial parent law. In recent probabilistic modeling, they are used to represent observed counts as a minimum, median, maximum, or other order statistic of latent iid counts, especially to model conditional underdispersion that is difficult to capture with a Poisson likelihood [2507.09032]. A closely related but distinct line of work studies the negative-binomial distribution of order \(k\), which is a waiting-time distribution for runs in Bernoulli sequences rather than an order statistic of iid negative binomial draws; its mode theory has been developed in detail [1702.02183]. A further generalization replaces iid Bernoulli trials by a \(q\)-Bernoulli scheme with geometrically varying success probabilities, producing several \(q\)-negative binomial distributions of order \(k\) [2210.03617].

## 1. Terminological scope and neighboring usages

The phrase “negative binomial order statistics” is used most directly for order statistics of iid negative binomial variables, but adjacent literature uses “order” in a different sense to denote waiting-time laws associated with runs of successes. The distinction is structural rather than cosmetic: the first construction transforms a parent count distribution by ranking iid samples, whereas the second defines a count through the occurrence time of combinatorial patterns in Bernoulli trials [2507.09032].

| Usage | Defining object | Representative source |
|---|---|---|
| Negative binomial order statistic | \(X_{(r)}\), the \(r\)-th smallest of \(D\) iid negative binomial draws | 2507.09032 |
| Negative-binomial distribution of order \(k\), type I | Waiting time on support \(kr,kr+1,\dots\) | 1702.02183 |
| \(q\)-negative binomial distributions of order \(k\) | Waiting times under geometrically varying success probabilities | 2210.03617 |

A common terminological confusion is to identify the “negative-binomial distribution of order \(k\)” with an order statistic of a negative binomial sample. The sources distinguish them sharply. The order-statistic construction starts from iid \(NB_{(\alpha,p)}\) variables and then ranks them. The order-\(k\) constructions start from binary trial sequences and count the trial at which the \(r\)-th run event occurs. This suggests that “order” refers to ranking in one literature and to run length or pattern structure in the other.

## 2. Distributional theory for iid negative binomial order statistics

In the formulation developed for count modeling, \(NB_{(\alpha,p)}\) denotes the usual “count-of-failures” negative binomial with shape parameter \(\alpha>0\) and failure probability \(p\in(0,1)\), with
\[
P(Z=k)=\bar p^\alpha (1-\bar p)^k \binom{\alpha+k-1}{k},
\qquad \bar p=1-p.
\]
Fix \(D\in\mathbb N\) and \(r\in\{1,\dots,D\}\), let \(Z_1,\dots,Z_D\) be iid \(NB_{(\alpha,p)}\), and write
\[
X_{(1)}\le \cdots \le X_{(D)}
\]
for their order statistics. The object of interest is
\[
X_{(r)}=\text{the \(r\)-th-smallest of } Z_1,\dots,Z_D.
\]

Its distribution admits the standard binomial-sum representation
\[
F^{(r,D)}_{NB}(k)
=P[X_{(r)}\le k]
=\sum_{t=r}^D \binom{D}{t}[F_{NB}(k)]^t[1-F_{NB}(k)]^{D-t},
\]
where the parent CDF has the incomplete-beta representation
\[
F_{NB}(k)
=\sum_{x=0}^k \binom{\alpha+x-1}{x}\bar p^\alpha (1-\bar p)^x
=I_p(k+1,\alpha)
=\frac{B(p;k+1,\alpha)}{B(k+1,\alpha)}.
\]
The pmf then follows by differencing:
\[
P[X_{(r)}=k]=F^{(r,D)}_{NB}(k)-F^{(r,D)}_{NB}(k-1).
\]

The same framework gives a joint density for the full vector of order statistics on \((k_1\le\cdots\le k_D)\):
\[
P(X_{(1)}=k_1,\dots,X_{(D)}=k_D)=D!\prod_{d=1}^D P(Z=k_d).
\]
Within the modeling program of discrete order statistics, these formulas make the negative binomial order statistic accessible through the parent negative binomial CDF and its incomplete-beta form [2507.09032].

## 3. Moments, dispersion, and the underdispersion mechanism

Closed-form elementary expressions for the first two moments are not available in finite terms, but both moments admit summation formulas in terms of the CDF of \(X_{(r)}\). Defining
\[
S_1
=E[X_{(r)}]
=\sum_{k=0}^\infty P[X_{(r)}>k]
=\sum_{k=0}^\infty [1-F^{(r,D)}_{NB}(k)],
\]
and
\[
S_2
=E[X_{(r)}^2]
=\sum_{k=0}^\infty (2k+1)P[X_{(r)}>k],
\]
one has
\[
\operatorname{Var}[X_{(r)}]=S_2-S_1^2.
\]
In practice, the infinite sums are truncated at \(k\) large enough that \(F^{(r,D)}_{NB}(k)\approx 1\) [2507.09032].

The principal dispersion summary is the index of dispersion
\[
\mathcal D[X_{(r)}]=\frac{\operatorname{Var}[X_{(r)}]}{E[X_{(r)}]}.
\]
As \(\alpha\to\infty\), which is the Poisson limit of the negative binomial, the index converges to
\[
\lim_{\alpha\to\infty}\mathcal D[X_{(r)}]
=\frac{1}{p}\operatorname{Var}_{U\sim N(0,1)^{(r,D)}}[U].
\]
For fixed \(p\) and large \(\alpha\), the limiting constant depends only on \((r,D)\) and scales like \(1/p\). The source states that by changing \(p\) one can achieve arbitrarily large overdispersion, and by choosing \(D\) one can dial in a lower bound on underdispersion. For the median, with \(r\simeq D/2\) and odd \(D\), \(\operatorname{Var}_{N(0,1)^{(r,D)}}\) is nonincreasing in \(D\), so larger \(D\) yields more underdispersion [2507.09032].

This mechanism is central to the modeling motivation. The Poisson distribution is the default choice of likelihood for probabilistic models of count data, but its equidispersion constraint can make predictive uncertainty artificially inflated. Order-statistic models alter dispersion without abandoning a standard parent count law. A plausible implication is that they supply a structured route from ordinary count likelihoods to more regular conditional behavior by compressing latent variability through ranking rather than by modifying only the parent variance function.

## 4. Latent-variable construction, inference, and computational structure

To embed negative binomial order statistics in hierarchical models, latent parent draws are introduced. For each observed count \(Y_i\equiv X_{(r_i)}\), one posits latent variables \(Z_{i,1},\dots,Z_{i,D_i}\) and uses the complete-conditional factor
\[
P(Z_{i,1:D_i}\mid X_{(r_i)}=Y_i,\alpha_i,p_i)
\propto
1\{Y_i=\text{ the }r_i^{\text{th}}\text{-smallest of }Z_{i,1:D_i}\}
\cdot \prod_{d=1}^{D_i} P_{NB}(Z_{i,d};\alpha_i,p_i).
\]
The resulting augmentation is designed to be modular with existing tools tailored to the parent distribution [2507.09032].

Sampling the latent \(Z\)'s is carried out by an exact \(O(D_i)\) algorithm based on three-way indicators
\[
C_{i,d}\in\{<Y_i,=Y_i,>Y_i\},
\]
together with sufficient statistics
\[
N_d=(\#\text{'s }<Y_i,\ \#=Y_i,\ \#>Y_i)
\]
up to step \(d\). The procedure sequentially draws \(C_{i,d}\) from a three-category distribution whose weights are ratios of partial order-statistic CDFs. Given \(C_{i,1:d}\), one samples \(Z_{i,d}\) from the negative binomial truncated to the appropriate region, and the loop may stop early if break-conditions on \(N_d\) are satisfied.

Once \(Z_{i,1:D_i}\) has been completed, the latent draws are iid \(NB(\alpha_i,p_i)\), so standard negative-binomial augmentations become available. In particular, the “CRT + Poisson” augmentation gives
\[
Z_{i,d}\sim NB(\alpha_i,p_i)
\iff
\left\{
L_{i,d}\mid Z_{i,d}\sim CRT(Z_{i,d},\alpha_i),\quad
Z_{i,d}\mid L_{i,d}\sim LogSeries(p_i)
\right\},
\]
with marginal
\[
L_{i,d}\sim Pois\!\left(\alpha_i\cdot \log\!\frac{1}{1-p_i}\right).
\]
Summing over \(d=1,\dots,D_i\) yields
\[
L_i=\sum_d L_{i,d}\sim Pois\!\left(D_i\alpha_i \ln\!\frac{1}{1-p_i}\right),
\qquad
Z_i=\sum_d Z_{i,d}\sim NB(D_i\alpha_i,p_i).
\]
The corresponding Gibbs updates are
\[
p_i\mid\{Z_{i,\cdot},\alpha_i\}\sim Beta\!\left(e+\sum_d Z_{i,d},\,f+D_i\alpha_i\right),
\]
\[
\alpha_i\mid\{L_i\}\sim Gamma\!\left(a+\sum_d L_{i,d},\,b+D_i\ln\!\frac{1}{1-p_i}\right).
\]
If \(D_i\) itself is random, with a prior such as shifted-Binomial or “OddBinomial,” its complete-conditional is proportional to
\[
P(D_i=d)\cdot P_{NB^{(r_i,d)}}(Y_i;\alpha_i,p_i).
\]

Several practical heuristics are emphasized. “Median” order statistics, with \(r\approx (D+1)/2\), tend to be most symmetric and concentrate around the negative binomial mean; they can be viewed as a discrete analog of location-scale families. Min and Max models yield heavy one-sided skew, and MaxOS admits a shortcut in sparse data because \(Y_i=0\) forces all \(Z\)'s to be \(0\). Worst-case sampling cost is \(O(D_i)\) per data point, but break-conditions often halt the loop early. Because inference re-uses standard conjugacies for negative binomial or Poisson models, regression on \(\alpha_i\) or \(p_i\), random effects, factor models, and spatio-temporal structure can be added by incorporating the \(Z\)- and \(L\)-augmentations. The framework has been illustrated on commercial flight times, COVID-19 case counts, Finnish bird abundance, and RNA sequencing data [2507.09032].

## 5. Negative-binomial distributions of order \(k\): type I and mode theory

The negative-binomial distribution of order \(k\), type I, is a natural generalization of the ordinary negative binomial, recovered by setting \(k=1\), and of the geometric distribution of order \(k\), recovered by setting \(r=1\). Here \(k\) and \(r\) are fixed positive integers, \(0<p<1\), and \(q:=1-p\). A discrete random variable \(X\) has this distribution if its support is
\[
n=kr,kr+1,kr+2,\dots
\]
and its pmf is
\[
P\{X=n\}=P_n
=
p^r
\sum_{\substack{n_1,\dots,n_k\ge 0\\ n_1+2n_2+\cdots+kn_k=n-kr}}
\frac{(n_1+\cdots+n_k+r-1)!}{n_1!\cdots n_k!(r-1)!}
\,q^{\,n_1+\cdots+n_k}.
\]
By convention \(P_n=0\) for \(n<kr\), and
\[
P_{kr}=p^r.
\]
This law is a waiting-time distribution rather than an order statistic of iid negative binomial samples [1702.02183].

Its analysis is based on the recurrence, valid for \(n>kr\),
\[
(n-kr)P_n
=
q\sum_{j=1}^k \bigl[n-kr+j(r-1)\bigr]p^{j-1}P_{n-j}.
\]
Together with the initial conditions \(P_n=0\) for \(n<kr\) and \(P_{kr}=p^r\), this recurrence completely determines the sequence \(\{P_n\}\). Let \(m_k(r,p)\) denote any mode, that is, any value \(n\) at which \(P_n\) attains its global maximum. Substituting \(n=m\) into the recurrence yields the upper bound
\[
m_k(r,p)\le
kr+\left\lfloor
\frac{(r-1)\bigl[1-p^k(1+kq)\bigr]}{q\,p^k}
\right\rfloor.
\]

When \(r=1\), the upper bound collapses to \(k\), and since necessarily \(m_k(1,p)\ge kr=k\), one obtains the exact identity
\[
m_k(1,p)=k.
\]
Thus the geometric distribution of order \(k\) always has its maximum exactly at \(n=k\) [1702.02183].

A complementary lower bound is obtained by studying successive differences \(\Delta_n:=P_n-P_{n-1}\). For \(k,r\ge 2\) and \(0<p\le (r-1)/r\), the lower bound takes the form
\[
m_k(r,p)\ge kr+\lfloor v_2\rfloor,
\]
where \(v_2\) is the larger real root of a quadratic \(f(v)=0\) whose leading coefficient is \(-qp^k<0\). The source states that the upper and lower bounds “pinch” the mode between two explicit integers and often coincide, yielding an exact formula. In the special case \(p=\tfrac12\), the mode can indeed be written in closed form:
\[
m_k\!\left(r,\tfrac12\right)
=
kr+(r-1)\bigl(2^{k+1}-k-2\bigr)+\delta_{k,r},
\]
where
\[
\delta_{k,r}=
\begin{cases}
0,\pm 1, & (k,r)=(2,2),\\
0, & \text{otherwise.}
\end{cases}
\]
The exceptional \(\pm 1\) reflects a flat region at \(n=6,7,8\) when \((k,r)=(2,2)\).

Distributions of order \(k\) arise naturally whenever one counts overlapping runs or patterns in sequences of Bernoulli trials. The source notes applications in designing acceptance-sampling plans, computing confidence intervals in runs-based tests, and sensitive-survey methodology where runs-of-successes are censored. In that setting, locating the mode is operationally important because it identifies the most-likely count and can guide thresholds, control-limits, or penalty levels [1702.02183].

## 6. \(q\)-negative binomial distributions of order \(k\) and generalized run schemes

A further branch of the subject studies variations of negative-binomial-of-order-\(k\) waiting-time laws under a \(q\)-Bernoulli trial model. In this scheme, the probability of success in the \((n+1)\)-st trial depends on the number of failures already observed:
\[
P(X_{n+1}=1\mid F_n=i)=\theta q^i,
\qquad
P(X_{n+1}=0\mid F_n=i)=1-\theta q^i,
\]
with \(0<\theta<1\) and \(0\le q<1\). As \(q\to 1\), the model collapses to iid Bernoulli\((\theta)\), and one recovers the classical negative-binomial-of-order-\(k\) distributions [2210.03617].

The literature distinguishes four main types, together with an \(\ell\)-overlapping generalization.

| Model | Interpretation | Support |
|---|---|---|
| \(W_{r,k}^{(1)}\) | \(r\)-th non-overlapping run of \(k\) successes | \(\{kr,kr+1,\dots\}\) |
| \(W_{r,k}^{(2)}\) | \(r\)-th run of length at least \(k\) | \(\{r(k+1)-1,r(k+1),\dots\}\) |
| \(W_{r,k}^{(3)}\) | \(r\)-th overlapping run of length \(k\) | \(\{k+r-1,k+r,\dots\}\) |
| \(W_{r,k}^{(4)}\) | \(r\)-th run of exactly \(k\) successes | \(\{r(k+1)-1,r(k+1),\dots\}\) |
| \(W_{r,k,\ell}\) | \(r\)-th \(\ell\)-overlapping run | \(\{\ell+r(k-\ell),\ell+r(k-\ell)+1,\dots\}\) |

Each type has a combinatorial kernel: \(A_q^k\) for non-overlapping runs, \(B_q^k\) for runs of length at least \(k\), \(C_q^k\) for overlapping runs, \(D_q^k\) for runs of exactly \(k\), and \(E_q^{k,\ell}\) for the \(\ell\)-overlapping generalization. These kernels admit “add one more cell” recurrences of the generic form
\[
\text{kernel}(r,s,t)=\sum_{j=0}^{\dots} q^{(r-1)j}\,\text{kernel}(r-1,s-j,\dots).
\]
The associated pmfs are expressed as combinatorial sums whose weights combine factors of the form
\[
\prod_{j=1}^{i}(1-\theta q^{j-1})
\]
with additional powers of \(q\) determined by weighted sums of cell occupancies.

The source identifies specific authors for the various types: Yalçın-Eryilmaz (2014) for Type 1, Makri-Psillakis (2016) for Type 2, Yalçın (2013) for Type 3, Oh-Jang (2022) for Type 4, and Kinacı-Coşkun-Karakaya-Akdoğan (2016) for the \(\ell\)-overlapping case. Ordinary generating functions in \(u\) may also be written by summing \(u^nP(W=n)\); by duality with \(q\)-binomial-of-order-\(k\) counting distributions, these have rational forms in \(u\) and \(q\), although the source describes the combinatorial sums as the most compact closed forms. As \(q\to1\), the \(q\)-factors disappear and the models revert exactly to the classical negative binomial distributions of order \(k\) [2210.03617].

The coexistence of these frameworks clarifies the modern landscape. Negative binomial order statistics in the strict sense are rank transforms of iid negative binomial counts and are motivated by tractable modeling of dispersion, especially underdispersion. Negative-binomial distributions of order \(k\) and their \(q\)-extensions are run-count waiting-time laws. The shared terminology reflects a common concern with structured count distributions, but the underlying stochastic constructions are fundamentally different.

Source: https://www.emergentmind.com/topics/negative-binomial-order-statistics