损失次数模型-负二项分布

——非寿险精算的基本理论

1、定义

一个成功概率为ppp的伯努利试验，不断重复，直至失败rrr次。此时成功的次数为一个随机变量，用XXX表示。称XXX服从负二项分布，记作X∼NB(r,p)X \sim NB(r,p)X∼NB(r,p)。

2、概率密度函数

P(X=k)=Ck+r−1kpk(1−p)rP(X=k)=C_{k+r-1}^kp^k(1-p)^rP(X=k)=Ck+r−1kpk(1−p)r

这里涉及到高中的排列组合知识，忘记的的同学可能要去复习一下了，哈哈哈~

公式释义：

rrr：一次实验中，失败次数；

kkk：一次实验中，失败rrr次时，成功的次数；

k+r−1k+r-1k+r−1：k+rk+rk+r是一次实验中总的实验次数，成功的次数加上失败的次数。因为目标失败次数为rrr次，当失败次数达到目标失败次数时，则停止实验，因此最后一次实验一定是失败的，是一个确定事件。所以随机事件只有k+r−1k+r-1k+r−1次；

ppp：一次实验成功的概率；

1−p1-p1−p：一次实验失败的概率；

P(X=k)P(X=k)P(X=k)：一次实验中，失败次数为rrr，成功次数为kkk时的概率；

例子

在一组抛硬币实验中，每次硬币正面向上的概率为p=0.5p=0.5p=0.5，则反面向上的概率为(1−p)(1-p)(1−p)，当反面向上的次数为r=5r=5r=5时停止实验，则在这组实验中正面次数出现k=5k=5k=5次的概率是多少？

<直接带入公式求解>

P(X=5)=C5+5−150.55(1−0.5)5=0.1230P(X=5)=C_{5+5-1}^50.5^5(1-0.5)^5=0.1230P(X=5)=C5+5−150.55(1−0.5)5=0.1230

补充一个公式：Cnm=n!m!(n−m)!C_n^m=\frac{n!}{m!(n-m)!}Cnm=m!(n−m)!n!

<公式理解>

上面的案例可以转化一下，在10次抛硬币实验中，反面向上的次数为5次，并且最后一次为反面向上。则，这一事件发生的概率是多少？

首先，5次反面向上的概率为(1−p)5(1-p)^5(1−p)5，剩余5次正面向上的概率为p5p^5p5。

其次，5次反面向上有4次发生的顺序是随机的，可能在第1，2，3，4次抛硬币时就出现了反面，在第10次才出现第五次反面，也可能在第6，7，8，9，10次抛硬币时出现了反面。但，第10次抛硬币时一定出现的是反面。所以，前四次反面出现的时机，有可能是第1至9次抛硬币实验中的任意4次。反面是前9次出现4次，则正面就是前9次出现5次，这种组合方式共有C95C_9^5C95种。每一种的概率都为p5(1−p)5p^5(1-p)^5p5(1−p)5，所以C95C_9^5C95种组合的概率就为C95p5(1−p)5C_9^5p^5(1-p)^5C95p5(1−p)5。

3、均值和方差

E(X)=rp1−pE(X)=\frac{rp}{1-p}E(X)=1−prp

D(X)=rp(1−p)2D(X)=\frac{rp}{(1-p)^2}D(X)=(1−p)2rp

3.1、均值E(X)E(X)E(X)的推导

E(X)=∑xp(x)E(X)=\sum xp(x)E(X)=∑xp(x)

这个公式在前面介绍泊松分布的时候提及了一下，但并没有展开说明，有兴趣的可以再去回顾一下概率论知识。

=∑xCx+r−1xpx(1−p)r=\sum xC_{x+r-1}^xp^x(1-p)^r=∑xCx+r−1xpx(1−p)r

把P(X)=Cx+r−1xpx(1−p)rP(X)=C_{x+r-1}^xp^x(1-p)^rP(X)=Cx+r−1xpx(1−p)r带入上式便可得到，上述公式

=∑x(x+r−1)!x!(r−1)!px(1−p)r=\sum x\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^r=∑xx!(r−1)!(x+r−1)!px(1−p)r

这是因为Cnm=n!m!(n−m)!C_n^m=\frac{n!}{m!(n-m)!}Cnm=m!(n−m)!n!，所以Cx+r−1x=(x+r−1)!x!(r−1)!C_{x+r-1}^x=\frac{(x+r-1)!}{x!(r-1)!}Cx+r−1x=x!(r−1)!(x+r−1)!

然后，x!=x(x−1)!x!=x(x-1)!x!=x(x−1)!，(r−1)!=r(r−1)!r(r-1)!=\frac{r(r-1)!}{r}(r−1)!=rr(r−1)!，px=ppx−1p^x=pp^{x-1}px=ppx−1，(1−p)r=(1−p)r+11−p(1-p)^r=\frac{(1-p)^{r+1}}{1-p}(1−p)r=1−p(1−p)r+1，带入上式

=∑x(x+r−1)!x(x−1)!r(r−1)!rppx−1(1−p)r+11−p=\sum x\frac{(x+r-1)!}{x(x-1)!\frac{r(r-1)!}{r}}pp^{x-1}\frac{(1-p)^{r+1}}{1-p}=∑xx(x−1)!rr(r−1)!(x+r−1)!ppx−11−p(1−p)r+1

将p和11−pp和\frac{1}{1-p}p和1−p1提到∑\sum∑前面，xxx和xxx约掉，并且(x+r−1)!x(x−1)!r(r−1)!r\frac{(x+r-1)!}{x(x-1)!\frac{r(r-1)!}{r}}x(x−1)!rr(r−1)!(x+r−1)!分子分母同时乘以rrr，得到

=p1−p∑(x+r−1)!r(x−1)!r(r−1)!rrpx−1(1−p)r+1=\frac{p}{1-p}\sum \frac{(x+r-1)!r}{(x-1)!\frac{r(r-1)!}{r}r}p^{x-1}{(1-p)^{r+1}}=1−pp∑(x−1)!rr(r−1)!r(x+r−1)!rpx−1(1−p)r+1

将(x+r−1)!r(x−1)!r(r−1)!rr\frac{(x+r-1)!r}{(x-1)!\frac{r(r-1)!}{r}r}(x−1)!rr(r−1)!r(x+r−1)!r分子的rrr提到∑\sum∑前面，分母中的rrr约掉，得到

=rp1−p∑(x+r−1)!(x−1)!r(r−1)!px−1(1−p)r+1=\frac{rp}{1-p}\sum \frac{(x+r-1)!}{(x-1)!r(r-1)!}p^{x-1}{(1-p)^{r+1}}=1−prp∑(x−1)!r(r−1)!(x+r−1)!px−1(1−p)r+1

因r(r−1)!=r!=((r+1)−1)!，(x+r−1)!=((x−1)+(r+1)−1)!r(r-1)!=r!=((r+1)-1)!，(x+r-1)!=((x-1)+(r+1)-1)!r(r−1)!=r!=((r+1)−1)!，(x+r−1)!=((x−1)+(r+1)−1)!，所以上式可写成

=rp1−p∑((x−1)+(r+1)−1)!(x−1)!((r+1)−1)!px−1(1−p)r+1=\frac{rp}{1-p}\sum \frac{((x-1)+(r+1)-1)!}{(x-1)!((r+1)-1)!}p^{x-1}{(1-p)^{r+1}}=1−prp∑(x−1)!((r+1)−1)!((x−1)+(r+1)−1)!px−1(1−p)r+1

另，x−1=m,r+1=nx-1=m,r+1=nx−1=m,r+1=n，则上式变为

=rp1−p∑(m+n−1)!m!(n−1)!pm(1−p)n=\frac{rp}{1-p}\sum \frac{(m+n-1)!}{m!(n-1)!}p^{m}{(1-p)^{n}}=1−prp∑m!(n−1)!(m+n−1)!pm(1−p)n

∑(m+n−1)!m!(n−1)!pm(1−p)n\sum \frac{(m+n-1)!}{m!(n-1)!}p^{m}{(1-p)^{n}}∑m!(n−1)!(m+n−1)!pm(1−p)n，这个公式是不是和p(x)=(x+r−1)!x!(r−1)!px(1−p)rp(x)=\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^rp(x)=x!(r−1)!(x+r−1)!px(1−p)r一摸一样，因此它是一个负二项分布。在所有情况下，所有事件的概率和为1，既∑p(x)=1\sum p(x)=1∑p(x)=1，所以∑(m+n−1)!m!(n−1)!pm(1−p)n=1\sum \frac{(m+n-1)!}{m!(n-1)!}p^{m}{(1-p)^{n}}=1∑m!(n−1)!(m+n−1)!pm(1−p)n=1

=rp1−p=\frac{rp}{1-p}=1−prp

3.2、方差D(X)D(X)D(X)的推导

根据之前泊松分布里面的文章，我们知道D(X)=E(X2)−(E(X))2D(X)=E(X^2)-(E(X))^2D(X)=E(X2)−(E(X))2，现在E(X)E(X)E(X)已经求出，因此只需推导E(X2)E(X^2)E(X2)即可。

E(X2)=∑x2p(x)E(X^2)=\sum x^2p(x)E(X2)=∑x2p(x)

这个公式我也曾在泊松分布这篇文章里面提到过，有兴趣的可以去查看。

把p(x)=(x+r−1)!x!(r−1)!px(1−p)rp(x)=\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^rp(x)=x!(r−1)!(x+r−1)!px(1−p)r带入上式

=∑x2(x+r−1)!x!(r−1)!px(1−p)r=\sum x^2\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^r=∑x2x!(r−1)!(x+r−1)!px(1−p)r

把px=ppx−1,(1−p)r=(1−p)r+11−p,(x+r−1)!x!(r−1)!分子和分母同时乘以r，p^x=pp^{x-1},(1-p)^r=\frac{(1-p)^{r+1}}{1-p},\frac{(x+r-1)!}{x!(r-1)!}分子和分母同时乘以r，px=ppx−1,(1−p)r=1−p(1−p)r+1,x!(r−1)!(x+r−1)!分子和分母同时乘以r，带入上式

=∑x2r(x+r−1)!x!r(r−1)!ppx−1(1−p)r+11−p=\sum x^2\frac{r(x+r-1)!}{x!r(r-1)!}pp^{x-1}\frac{(1-p)^{r+1}}{1-p}=∑x2x!r(r−1)!r(x+r−1)!ppx−11−p(1−p)r+1

把r,p,11−pr,p,\frac{1}{1-p}r,p,1−p1提到∑\sum∑前面得

=rp1−p∑x2(x+r−1)!x!r(r−1)!px−1(1−p)r+1=\frac{rp}{1-p}\sum x^2\frac{(x+r-1)!}{x!r(r-1)!}p^{x-1}(1-p)^{r+1}=1−prp∑x2x!r(r−1)!(x+r−1)!px−1(1−p)r+1

因x2=xx,x!=x(x−1)!,r(r−1)!=r!x^2=xx,x!=x(x-1)!,r(r-1)!=r!x2=xx,x!=x(x−1)!,r(r−1)!=r!所以

=rp1−p∑xx(x+r−1)!x(x−1)!r!px−1(1−p)r+1=\frac{rp}{1-p}\sum xx\frac{(x+r-1)!}{x(x-1)!r!}p^{x-1}(1-p)^{r+1}=1−prp∑xxx(x−1)!r!(x+r−1)!px−1(1−p)r+1

因xxx分子和分母可以约掉，所以

=rp1−p∑x(x+r−1)!(x−1)!r!px−1(1−p)r+1=\frac{rp}{1-p}\sum x\frac{(x+r-1)!}{(x-1)!r!}p^{x-1}(1-p)^{r+1}=1−prp∑x(x−1)!r!(x+r−1)!px−1(1−p)r+1

令y=x−1,n=r+1y=x-1,n=r+1y=x−1,n=r+1则

=rp1−p∑(y+1)(y+n−1)!y!(n−1)!py(1−p)n=\frac{rp}{1-p}\sum (y+1)\frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n}=1−prp∑(y+1)y!(n−1)!(y+n−1)!py(1−p)n

把∑(y+1)(y+n−1)!y!(n−1)!py(1−p)n\sum (y+1)\frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n}∑(y+1)y!(n−1)!(y+n−1)!py(1−p)n拆开

=rp1−p(∑y(y+n−1)!y!(n−1)!py(1−p)n+∑(y+n−1)!y!(n−1)!py(1−p)n)=\frac{rp}{1-p}(\sum y\frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n}+\sum \frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n})=1−prp(∑yy!(n−1)!(y+n−1)!py(1−p)n+∑y!(n−1)!(y+n−1)!py(1−p)n)

∑y(y+n−1)!y!(n−1)!py(1−p)n\sum y\frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n}∑yy!(n−1)!(y+n−1)!py(1−p)n这个形式是不是和我们上面提到的E(X)=∑x(x+r−1)!x!(r−1)!px(1−p)r=rp1−pE(X)=\sum x\frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^r=\frac{rp}{1-p}E(X)=∑xx!(r−1)!(x+r−1)!px(1−p)r=1−prp一模一样。所以，∑y(y+n−1)!y!(n−1)!py(1−p)n=E(y)=np1−p\sum y\frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n}=E(y)=\frac{np}{1-p}∑yy!(n−1)!(y+n−1)!py(1−p)n=E(y)=1−pnp，又因为n=r+1n=r+1n=r+1，所以，E(y)=np1−p=(r+1)p1−pE(y)=\frac{np}{1-p}=\frac{(r+1)p}{1-p}E(y)=1−pnp=1−p(r+1)p。再看∑(y+n−1)!y!(n−1)!py(1−p)n)\sum \frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n})∑y!(n−1)!(y+n−1)!py(1−p)n)是不是和∑p(x)=∑(x+r−1)!x!(r−1)!px(1−p)r=1\sum p(x)=\sum \frac{(x+r-1)!}{x!(r-1)!}p^x(1-p)^r=1∑p(x)=∑x!(r−1)!(x+r−1)!px(1−p)r=1一模一样，所以，∑(y+n−1)!y!(n−1)!py(1−p)n)=1\sum \frac{(y+n-1)!}{y!(n-1)!}p^{y}(1-p)^{n})=1∑y!(n−1)!(y+n−1)!py(1−p)n)=1，带入上式

=rp1−p((r+1)p1−p+1)=\frac{rp}{1-p}(\frac{(r+1)p}{1-p}+1)=1−prp(1−p(r+1)p+1)

=r2p2+rp(1−p)2=\frac{r^2p^2+rp}{(1-p)^2}=(1−p)2r2p2+rp

经过上面的层层推导，我们终于把E(X2)=r2p2+rp(1−p)2E(X^2)=\frac{r^2p^2+rp}{(1-p)^2}E(X2)=(1−p)2r2p2+rp推导出来了，再根据E(X)=rp1−pE(X)=\frac{rp}{1-p}E(X)=1−prp，便可以求出D(X)D(X)D(X)的值

D(X)=E(X2)−(E(X))2D(X)=E(X^2)-(E(X))^2D(X)=E(X2)−(E(X))2

=r2p2+rp(1−p)2−(rp1−p)2=\frac{r^2p^2+rp}{(1-p)^2}-(\frac{rp}{1-p})^2=(1−p)2r2p2+rp−(1−prp)2

=rp(1−p)2=\frac{rp}{(1-p)^2}=(1−p)2rp

4、负二项分布在精算中的应用

4.1、定义

假设损失次数NNN服从参数为r和βr和\betar和β的负二项分布，则发生kkk次损失的概率为
P(N=k)=Ck+r−1k(β1+β)k(11+β)rP(N=k)=C_{k+r-1}^k({\frac{\beta}{1+\beta}})^k(\frac{1}{1+\beta})^rP(N=k)=Ck+r−1k(1+ββ)k(1+β1)r

连接上文，这里的β=p1−p\beta=\frac{p}{1-p}β=1−pp

负二项分布的均值和方差为：

E(N)=rβE(N)=r\betaE(N)=rβ

D(N)=rβ(1+β)D(N)=r\beta(1+\beta)D(N)=rβ(1+β)

将β\betaβ替换为ppp之后，我们可以发现和上面推导的结果是一样的。

4.2、性质

（1）方差大于均值。在实际运用中如果方差大于均值，负二项分布比泊松分布可以更好的拟合数据。

（2）负二项分布是一种混合泊松分布，即如果假设泊松分布的参数服从伽马分布，由此得到的混合泊松分布即为负二项分布。
针对这个性质的证明，到时在学完伽马分布的时候再看吧。

5、软件实操

5.1、相同的ppp，不同的rrr

画图代码

import numpy as np
import math
import matplotlib.pyplot as pltdef nbpmf(sample,r,p):'''Parameters----------sample : int负二项分布的样本数量.r : int贝努力实验中的失败次数.p : float贝努力实验中的成功的概率.Returns-------k_sample : int样本数量.probability : TYPE概率.'''sample=sampler=rp=pk_sample=np.arange(sample)#会返回一个0至sample-1的连续整数列表probability=[]for k in k_sample:temp = (math.factorial(k+r-1)/(math.factorial(k)*math.factorial(r-1)))*p**k*(1-p)**r#计算负二项分布的概率#math.factorial()是阶乘函数probability.append(temp)return k_sample,probabilityr=5#可以跟换不同的r值查看概率图像
p=0.6
sample=40
k,probability=nbpmf(sample,r,p)#调用自定义的负二项分布函数
plt.subplot(1,1,1)#定义图片的布局，针对一张图片来说可以省略该代码
plt.plot(k,probability,"r",label="r="+str(r))#画曲线图，并定义曲线为红色"r"，并显示图例label（后面是图例名称）
plt.plot(k,probability,"ro")#画散点图，并定义散点为红色
plt.title("Negative Binomial distribution(p="+str(p)+")")#设置图片标题
plt.xlabel("Number of successes (k)")#设置图片x轴说明
plt.ylabel("probability")#设置图片y轴说明
plt.legend(loc=0)#将图例放在合适的位置，让系统自动查找

输出结果

从图片可以看出，当成功的概率ppp一定时，随着失败次数rrr的增大，函数图像逐渐向左偏移，峰值逐渐变小，图像由细高变得宽矮。

5.2、相同的rrr，不同的ppp

代码与上面相同，只是保持rrr不变，而修改了ppp值和图像标头

输出结果

从图片可以看出，其变化形态和5.1、相同的ppp，不同的rrr。是一样的。只是变化幅度有些许不同。

5.3、负二项分布与正态分布

画图代码

import numpy as np
import math
import matplotlib.pyplot as pltdef nbpmf(sample,r,p):'''Parameters----------sample : int负二项分布的样本数量.r : int贝努力实验中的失败次数.p : float贝努力实验中的成功的概率.Returns-------k_sample : int样本数量.probability : TYPE概率.'''sample=sampler=rp=pk_sample=np.arange(sample)#会返回一个0至sample-1的连续整数列表probability=[]for k in k_sample:temp = (math.factorial(k+r-1)/(math.factorial(k)*math.factorial(r-1)))*p**k*(1-p)**r#计算负二项分布的概率#math.factorial()是阶乘函数probability.append(temp)return k_sample,probabilitydef normalpdf(mu,sigma,simple):"""Parameters----------mu : float正态分布的均值.sigma : float正态分布标准差.simple : int生成的样本数量.Returns-------x : floatx轴值.probability : float取值为x轴值时，对应的概率密度函数值，相当于离散变量的概率."""mu = musigma = sigmasimple = simplex = np.arange(sample)probability=np.exp(-(x-mu)**2/(2*sigma**2))/(math.sqrt(2*math.pi)*sigma)#正态分布的概率密度函数return x,probabilityr=90#通过更改r值便可以画出不同参数的图像
p=0.5
plt.subplot(1,1,1)
mean = r*p/(1-p)#负二项分布的均值
sigam = math.sqrt(r*p/((1-p)**2))#负二项分布的标准差
sample=mean*2
k,probability_nb=nbpmf(sample,r,p)
k,probability_no=normalpdf(mean,sigam,sample)
plt.plot(k,probability_no,"r",label="normal")
plt.plot(k,probability_nb,"g--",label="n_binomial")
plt.title("Negative Binomial and Normal distribution" + "(r,p="+str(r)+","+str(p)+")")#设置图片标题
plt.xlabel("Number of successes (k)")#设置图片x轴说明
plt.ylabel("probability")#设置图片y轴说明
plt.legend(loc=0)#将图例放在合适的位置，让系统自动查找

输出结果

通过图片可以看出，在ppp值不变的情况下，随着失败次数rrr的逐渐变大，负二项分布逐渐趋近正态分布。所以在大样本情况下，可以用正态分布近似负二项分布。

—End—

*** 参考资料 ***

1、《非寿险精算学》孟生旺著

2、负二项分布_虚胖一场的博客-CSDN博客_负二项分布

3、概率质量函数(PMF)、概率密度函数(PDF)、累积分布函数(CDF)_Uncertainty!!的博客-CSDN博客_概率质量函数

4、排列组合cn和an公式_xxxx追风的博客-CSDN博客_排列组合

5、负二项分布的数学期望和方差的一种求法 - 百度文库 (baidu.com)

6、Python中如何求阶乘？教你四个方法-群英 (qycn.com)