本文是一篇手记，记录了我对负二项分布及其应用的理解。
目录如下：
- 1.理解”负”的含义
- 2.理解定义1：伯努利过程视角
- 3.理解其与泊松分布的极限关系
- 4.理解期望/方差/overdispersion之间的关系
- 5.理解定义2：gamma-poisson mixture distribution
- 6.counting modeling中的应用
- 7.理解其与gamma分布的极限关系
1.理解”负”的含义
知乎回答有提到：https://www.zhihu.com/question/24253978?sort=created
文档中有提到：http://www.johndcook.com/negativebinomial.pdf
负二项级数：https://brilliant.org/wiki/negative-binomial-theorem/
负是指负二项级数，但实际上，“负”除了告诉我们负二项分布的由来，对理解其直观意义并无帮助。
根据组合数定义：
\[\begin{aligned}\left(\begin{array}{c}x+r-1 \\x\end{array}\right) &=\frac{(x+r-1)(x+r-2) \cdots r}{x !} \\&=(-1)^{x} \frac{(-r-(x-1))(-r-(x-2)) \cdots(-r)}{x !} \\&=(-1)^{x} \frac{(-r)(-r-1) \cdots(-r-(x-1))}{x !} \\&=(-1)^{x}\left(\begin{array}{c}-r \\x\end{array}\right)\end{aligned}\]把组合数定义带入下面的式子（(1-q)^-r是负二项级数）：
\[\begin{aligned}1 &=p^{r} p^{-r} \\&=p^{r}(1-q)^{-r} \\&=p^{r} \sum_{x=0}^{\infty}\left(\begin{array}{c}-r \\x\end{array}\right)(-q)^{x} \\&=p^{r} \sum_{x=0}^{\infty}(-1)^{x}\left(\begin{array}{c}-r \\x\end{array}\right) q^{x} \\&=\sum_{x=0}^{\infty}\left(\begin{array}{c}x+r-1 \\x\end{array}\right) p^{r} q^{x}\end{aligned}\]2.理解定义1：伯努利过程视角
https://www.johndcook.com/blog/computing_binomial_coefficients/
定义1完全是基于负二项级数得到的，从伯努利过程的视角出发，也能自然的能理解负二项分布；
实际上，负二项分布描述的是第r次成功前失败的次数，记为k，那么有：
\[\text{Pr}(X=k) = \left(\begin{array}{c}k+r-1 \\r-1\end{array}\right) \cdot(1-p)^{k} p^{r}\]2.1.理解其作为sum of i.i.d geom
一旦理解了定义1，很自然的就可以理解为什么sum of i.i.d geom 是 负二项分布：
因为几何分布刻画的是在多次伯努利试验中，第一次成功前所经历的失败次数。
2.2.理解其期望/方差/MGF的推导
推导链接参考:
https://statisticalmodeling.wordpress.com/2011/07/28/the-negative-binomial-distribution/
理解了sum of i.i.d geom 是 NB后，NB的很多数字特征是可以通过geom推导出来的。
包括：
期望/方差/MGF
假设由r个i.i.d的几何分布相加，得到了NB(r,p)，那么期望与方差可以写作：
\[E(X)=\frac{r(1-p)}{p} \quad \operatorname{Var}(X)=\frac{r(1-p)}{p^{2}}\]MGF：
\[M(t)=\frac{p^{r}}{\left(1-(1-p) e^{t}\right)^{r}}\]只需要从几何分布的期望/方差/MGF出发就能得到上面的结果。
2.3.理解定义1的拓展
由于gamma函数在整数时恰好等于整数-1的阶乘，即：
\[\Gamma(a)=(a-1)!\]因此，可以把组合数的定义式利用gamma函数对定义1进行拓展，这样使得r得取值范围在整个正实数范围。即：
# failures before r-th success, denote by k：
\[f(k ; r, p) \equiv \operatorname{Pr}(X=k)=\frac{\Gamma(k+r)}{k ! \Gamma(r)}p^{r} (1-p)^{k} \quad \text { for } k=0,1,2, \dots\]3.理解其与泊松分布的极限关系
https://en.wikipedia.org/wiki/Negative_binomial_distribution#Poisson_distribution
令X~NB(r,p)，若 stopping parameter r趋于无穷，且每次成功的概率p趋于0，且保持期望不变（E(X)=r(1-p)/p），并把期望记作λ，那么λ就等于：
将p带入定义2的表达式，得失败次数为k的概率为：
\[P(X=k ; r, p)=\frac{\Gamma(k+r)}{k ! \cdot \Gamma(r)} p^{r}(1-p)^{k}\\=\frac{1}{k !} \cdot \frac{\Gamma(r+k)\lambda^{k}}{\Gamma(r)(r+\lambda)^{k}} \cdot (\frac{r}{r+\lambda})^r\\=\frac{\lambda^{k}}{k !} \cdot \frac{\Gamma(r+k)}{\Gamma(r)(r+\lambda)^{k}} \cdot \frac{1}{\left(1+\frac{\lambda}{r}\right)^{r}}\]现在，考虑r趋于无穷，第二项将会趋于1（分子分母同阶）;第三项根据compound interest limit 趋于e^{-λ}，所以得到：
\[\lim _{r \rightarrow \infty} P(X=k ; r, p)=\frac{\lambda^{k}}{k !} \cdot 1 \cdot \frac{1}{e^{\lambda}}\]即泊松分布的PMF。
换句话说，上面的r参数实际上控制了NB与Poisson的deviation；这个性质使得负二项式分布可以作为泊松的稳健替代，即：当r很大得时候就会逼近泊松分布，当r很小的时候，其方差比泊松分布大。【这一点可以从NB的mean与var中都能看出来，具体见下面】
4.理解期望/方差/overdispersion之间的关系
根据前面第2节的结果，期望与方差分别为：
\[\mu=E(X)=r \frac{1-p}{p}\\\sigma^2=Var(X)=r\frac{1-p}{p^2}\]重新改写一下方差，并以期望的形式表示（一点代数技巧）：
\[\sigma^{2}=\mu+\frac{1}{r} \mu^{2} > \mu\]我们发现，在NB中，方差总是大于期望的；而当r（stopping parameter）趋于无穷大的时候，方差与期望相等。前面我们已经推导了当r趋于无穷的时候，NB就成了泊松分布。这里再次印证了这个结论。
我们不妨把1/r称为dispersion parameter，它能够帮助我们检验数据的overdispersion情况(利用Wald test检验原假设1/r=0是否成立)。
这一点性质使得我们可以把NB当作一个更general得模型来用，在实际count modeling的时候比poisson适合面更广，因为它能够handle overdispersion！
5.理解定义2：gamma-poisson mixture distribution
https://statisticalmodeling.wordpress.com/2011/07/28/the-negative-binomial-distribution/
实际上，负二项分布也叫gamma-poisson mixture distribution；其由来如下：
令Y l λ ~ Pois(λ)，λ ~ gamma(ro,bo)，有X的边际分布服从NB，X ~ NB(r0,b0/(1+b0))；
推导如下：
\[\begin{aligned} P(Y=y) &=\int_{0}^{\infty} P(Y=y | \lambda) f_{0}(\lambda) d \lambda \\ &=\int_{0}^{\infty} \frac{e^{-\lambda }\lambda ^{y}}{y !} \frac{\left(b_{0} \lambda\right)^{r_{0}} e^{-b_{0} \lambda}}{\Gamma\left(r_{0}\right)} \frac{d \lambda}{\lambda}\\ &=\frac{ b_{0}^{r_{0}}}{y ! \Gamma\left(r_{0}\right)} \int_{0}^{\infty} e^{-\left(b_{0}+1\right) \lambda} \lambda^{r_{0}+y} \frac{d \lambda}{\lambda} \\ &=\frac{\Gamma\left(r_{0}+y\right)}{y ! \Gamma\left(r_{0}\right)} \frac{ b_{0}^{r_{0}}}{\left(b_{0}+1\right)^{r_{0}+y}} \int_{0}^{\infty} \frac{1}{\Gamma\left(r_{0}+y\right)} e^{-\left(b_{0}+1\right) \lambda}\left(\left(b_{0}+1\right) \lambda\right)^{r_{0}+y} \frac{d \lambda}{\lambda} \\ &=\frac{\Gamma(y+r_0)}{y ! \Gamma(r_0)}\left(\frac{1}{b_{0}+1}\right)^{y}\left(\frac{b_{0}}{b_{0}+1}\right)^{r_{0}} \end{aligned}\]这正好是负二项分布的PMF。
从这个定义中，我们可以猜一下：
the unconditional NB的方差直观上要比conditional poisson的大：由于NB分布中存在一个λ参数额外的不确定性，因此the unconditional NB的方差直观上要比conditional poisson的大。
这实际上是所有mixture distribution的特点，证明留到方差分解一节。
【注意这和第4节的结论并不一样：第4节的意思是说对于NB来说，其variance总是大于mean，而这个程度由参数r来控制；本节的这个猜测是从mixture分布的特点出发，认为the unconditional NB方差要大于conditional poisson的方差】
gamma-poisson mixture distribution的定义允许我们探究很多有意思的事实，因此，在进入方差分解一节以前，我们不妨先看看NB与gamma期望的关系。
5.1理解NB与gamma的期望关系：利用全期望公式
gamma-pois mixture model的定义:
若X l λ ~ Pois(λ)，λ ~ gamma(ro,bo)，有X的边际分布服从NB，X ~ NB(r0,b0/(1+b0))；
如果以定义1中NB(r0,p)的视角来看，可以定义1与定义2中p和b0的关系：
\[p=\frac{b_0}{1+b_0}\]根据定义1，NB(r0,p)的期望可以写作： \(\mu=r_0 \frac{1-p}{p}\)
带入p=b_0/(1+b_0)，可以得到： \(\mu=\frac{r_0}{b_0}\)
ro/bo恰好是gamma(ro,bo)的期望，这说明在gamma-pois mixture model的视角下，NB与gamma的期望是一致的。这是巧合吗？当然不是，如果你还记得全期望公式：
\[E(E(X|\lambda))=E(X) \quad(全期望公式)\\E(X|\lambda)=\lambda \quad(泊松分布的期望)\\\]所以：
\[E(X) = E(\lambda)=r_0/b_0\]可以很自然的发现，在gamma-pois mixture model的视角下，NB与gamma的期望是一致的。
5.2.variance decomposition of mixture distribution
https://statisticalmodeling.wordpress.com/2011/06/16/the-variance-of-a-mixture/
https://statisticalmodeling.wordpress.com/2011/07/28/the-negative-binomial-distribution/
1.gamma-poisson mixture distribution
若X l λ ~ Pois(λ)，λ ~ gamma(ro,bo)，有X的边际分布服从NB，X ~ NB(r0,b0/(1+b0))
2.gamma的方差与期望
本节的证明需要对gamma分布的期望与方差有了解。对于X~gamma(a,b)来说：
E(X)=a/b Var(X)=a/b^2
在gamma-poisson mixture distribution的视角下，我们猜测由于NB分布中存在λ这个额外的不确定性，因此NB的方差直观上看要比conditional泊松的大。我们的猜测是完全有依据的，这个依据可以通过方差分解公式(EVVE法则)来证实：
\[\operatorname{Var}(X)=E[\operatorname{Var}(X | \lambda)]+\operatorname{Var}[E(X | \lambda)]\]第一项：Var(X)是指unconditional NB的方差，根据定义：
\[\sigma^2=Var(X)=r\frac{1-p}{p^2}\]带入p=b_0/(1+b_0)，可以得到：
\[\sigma^2=Var(X)=r\frac{1-p}{p^2}\]第二项 由于Var(X l λ) = λ，所以E(Var(X l λ))=E(λ) = ro/bo = ro(1-p)/p（根据全期望公式得到）
第三项 Var(E(X l λ))=Var(λ)=ro/bo^2 ，带入bo=p/(1-p)，可以得到：
\[Var(E(X|λ))=r_0/b_0^2 = \frac{r_0 (1-p)^2}{p^2}\]第二项和第三项加起来正好等于第一项。
6.counting modeling中的应用
GLM介绍：https://www.youtube.com/watch?v=vpKpFMUMaVw
泊松回归介绍：https://www.youtube.com/watch?v=c73Q8nQ31qc
学习负二项分布的目的最终是为了应用，计数资料建模中最常见的分布是泊松/负二项分布，即：泊松/负二项回归。
在正式介绍这俩方法前，不妨先回顾一下线性回归以及GLM的概念。
6.1 回忆GLM
本节回忆几个基本问题。
- 线性回归中根据方程预测出来的是什么？
假设Y l X ~ N(μ,σ^2)，根据回归方程：E(Y l X=x)=μ = ax+b；我们预测出来的值本质上是条件期望！
- GLM是什么？
广义线性回归与线性回归也是类似的，只不过它做了一点泛化：
(1) 不再假设Y只能服从normal了，而是假设它服从指数族分布，但不会改变我们预测的仍然是条件期望E(Y l X)；
(2) 需要一个link function来连接解释变量与被解释变量，即：E(Y l X)=g(βX)，这个g()可以是指数函数，也可以是logisitc函数，等等；
- 泊松/负二项回归与GLM的关系？
泊松/负二项分布都属于指数族分布，而这俩分布的定义域都是非负整数；因此，可以假设Y l X服从泊松/负二项分布进行计数资料的预测；
具体来说，对于泊松回归，预测的值是E(Yi l X=xi)=λi；而λi恰好是泊松分布的参数，它反映的是在单位观测时间下，稀有事件的发生次数，因此泊松回归更直观；而对于负二项回归来说，在NB(r,p)中，r和p无法直观体现出期望与方差，因此，负二项分布的第一种定义对我们理解其回归的过程没有帮助；我们需要对NB进行重参数化，将r和p以期望和方差的形式来表达负二项分布的参数。
- 4.对于泊松/负二项回归来说，选择什么link function呢？
显然Y l X ≥0，而由于我们不对X做任何分布假设，因此βX的取值范围是整个实数集；这就使得如果直接利用βX预测得到的结果很多可能都会小于0，从而不具备解释性；因此，我们需要一个能把βX的取值范围限制在非负数的link function；很显然，指数函数e^(·)是一个好的选择。因此，我们选择e^(·)作为link function，即：E(Y l X)=e^(βX)。
6.2.reparameterized for counting modeling
对泊松分布来说，其参数λ，而期望=方差=λ；在泊松回归中，我们要预测的是正是λ【E(Yi l X=xi)=λi】。所以，泊松回归理解起来比较直观。
但对于负二项分布来说，其参数是NB(r,p)的r和p，和我们的计数建模有什么关系呢？所以，仅靠前面的2种定义对理解负二项回归是不够的，因为在回归中，相比于NB(r,p)的r和p，我们更关注条件期望E(Y l X)与条件方差Var(Y l X)。
换句话说，NB(r,p)中的参数最好能够用mean和variance重新表达。
回忆定义1的拓展（# failures before r-th success, denote by k）：
\[f(k ; r, p) \equiv \operatorname{Pr}(X=k)=\frac{\Gamma(k+r)}{k ! \Gamma(r)}p^{r} (1-p)^{k} \quad \text { for } k=0,1,2, \dots\]再顺便回忆一下期望与方差：
\[\mu=E(X)=r \frac{1-p}{p}\\\sigma^2=Var(X)=r\frac{1-p}{p^2}\]将r和p以μ和σ的形式表达，得到：
\[\begin{aligned}&p=\frac{\mu}{\sigma^{2}}\\&r=\frac{\mu^{2}}{\sigma^{2}-\mu}\\&\operatorname{Pr}(X=k)=\left(\begin{array}{c}k+\frac{\mu^{2}}{\sigma^{2}-\mu}-1 \\\frac{\mu^{2}}{\sigma^{2}-\mu}-1\end{array}\right)\left(\frac{\sigma^{2}-\mu}{\sigma^{2}}\right)^{k}\left(\frac{\mu}{\sigma^{2}}\right)^{\mu^{2} /\left(\sigma^{2}-\mu\right)}\end{aligned}\]有时候，也不那么关注其方差（例如：在负二项回归中，我们可能只关注条件期望），那么，可以将r和p以r和μ的形式表达，得到：
\[p=\frac{r}{r+\mu}\\Pr(X=k)=\frac{\Gamma(r+k)}{k ! \Gamma(r)}\left(\frac{r}{r+\mu}\right)^{r}\left(\frac{\mu}{r+\mu}\right)^{k} \quad \text { for } k=0,1,2, \dots\]6.3. 负二项回归
https://www.johndcook.com/blog/2009/11/03/negative-binomial-poisson-gamma/
负二项回归推导:http://www.karlin.mff.cuni.cz/~pesta/NMFM404/NB.html
即便对负二项分布进行了重参数化，其PMF还是非常抽象，对于直观理解其意义比较费劲。
现在，请你忘记负二项分布的标准定义，你对负二项分布的印象只保留在它是泊松分布一个很好的近似。此时，负二项分布的参数NB(r,p)已经不再重要，你在回归任务中，更加关注条件期望E(Y l X)与条件方差Var(Y l X)。
因此，我们的目的就是依据解释变量来预测被解释变量的条件期望：
\[\mathrm{E}\left[Y_{i} | X=\mathbf{x}_{i}\right]=e^{\mathbf{x}_{i}^{\top} \boldsymbol{\beta}}\\Y|X\sim NB(r,\frac{r}{r+\mu})\]后面就是利用MLE来求解相应的系数了： \(\prod_{i=1}^{n}\text{Pr}(y_i|x_i)\)
在实践中，负二项回归往往应用更广泛。
因为计数资料的条件方差往往大于条件期望，前面推导的结论发现负二项分布在r很大的时候，其实是泊松分布很好的近似，因此负二项回归应用更广泛。
6.4.基因组学中的应用
计数过程在基因组学中简直再常见不过，生成reads的过程就可以看作一个计数过程。
我们可以把采样自某个基因的reads数量建模为泊松分布，并利用泊松回归估计具体的数量。
然而，对于同一个基因来说，通常var比mean要大（比如进行了3次生物学重复实验，假设测了3只同样疾病状态的小鼠；注意，生物重复不是技术重复；这3次得到某基因的reads counts，一般来讲方差大于其均值），特别是对于高表达的基因。因此，负二项回归要更合适一点。
7.理解其与gamma分布的极限关系
考虑到伯努利过程的几何/二项分布对应了泊松过程的指数/泊松分布；我们自然有理由认为负二项分布(sum of i.i.d geom)相似的对应了gamma分布(sum of i.i.d exponential)；
上一篇: 闲扯概率论中的计数问题
下一篇: 如何快速发布一个python软件包