生物统计(3)-常见的离散型分布

前言

各种概率分布算是统计学中的一个非常重要的知识点，我在网上看到过一篇文章，题目是《常用连续型分布介绍及R语言实现》（网址为：http://blog./r-density/）。这一篇文章只是补充了一下常见的离散型分布介绍及R实现方式。

随机变量取一切可能值的概率的规律称为概率分布(probability distribution)，简称为分布，通常使用表格，图形，公式来表示这些分布，一个概率分布总是一和个总体(population)相关联的，这里的总体也称为样本空间(sample space)，它是由某种试验的所有可能结果构成的集合，这些结果的获取服从某种概率规律，因此， 一个总体（样本空间）是由一个取值范围及与之相连的概率所组成。于是，给出了这个概率分布就等于知道了总体，这些用数学语言表示的概率分布有一些理论参数，称为总体参数(population parameter)，对于一个样本来说，它有自己的一些参数，例如样本均值、样本方差、样本标准差、样本分位数等。相对于总体，它也有总体均值、总体方差、总体标准差，总体分位数等参数。只是对于不同的总体，这些参数可能不同，并且对于某些总体，某些参数可能就不存在。

离散型分布

离散变量很常见，例如在一批产品的质量检查中可能发现的次品数目、一日内通过某电话总体的通话次数、在若干次试验中成功的次数等，这些变量取离散值，而且按照某种概率来取每个值，取每个值的概率应该大于或等于0、小于或等于1，并且取各种值的概率之和应该等于1，用符号来表示就是，假定离散随机变量X的取值范围为x1,x2,…xn，并且取可能值xi的概率记为P(X=xi)，那么这些概率应该同时满足下面的关系：

mark

满足这样关系的概率P(xi)的总和就是该离散变量的概率分布。离散变量的取值不一定是有穷的，例如Poisson分布的取值范围就是所有的非负整数，因而有无穷多可能的值。

除了分布函数P(X=xi)外，还有累积分布函数(cumulative distribution function, CDF)的概念，离散随机变量的累积分布函数的定义如下所示：

mark

统计教科书上的各种概率分布表，通常都是累积分布函数的值，在计算时，这往往比原始的概率分布要方便，例如我要计算P(a<>

mark

或者是通过相减来计算，如下所示：

mark

显然，通过相减计算要容易得多，只需要找到F(b)和F(a)即可。

对于离散分布随机变量X来说，它的总体均值（通常用μ表示）定义如下所示：

mark

它又称为X的数学期望，这里的符号E(x)就是X的期望(expectation)的意思。

变量X的总体方差通常用δ^2表示，它的公式如下所示：

mark

这里的符号Var(x)就是X的方差(variance)的意思，总体标准差定义为总体方差的平方根，如下所示：

mark

显然，这里的均值、方差等概率和样本均值、样本方差很相似，只要把这里的xi换成样本点，而把P(xi)换成1/n（在样本方差经常使用1/(n-1)，在样本量大的时候1/n和1/(n-1)差别不大），就和样本量为n的情况下的样本均值和样本方差一致了，需要注意的是，总体均值和总体方差有可能是无穷或者不存在的（总体的数据很难得到，只能根据样本来推测总体），但样本均值和样本方差总是存在的。

下面看一些常见的离散型分布

二项分布

先看一个案例描述。

某项考虑有10道题，假定某考生答对每道题的概率是0.6，并假定各个题目的答案独立，而及格的标准则是至少答对6道题，那么至少答对6道题的概率就是事件“通过考试“的概率，这时，如果把答对题目的数目记为X，那么X是一个定量变量，而通过考虑的或者考试及格的概率就可以用以下公式来表示，如下所示：

mark

而没有一道题答对的概率为P(X=0)，这就是说10道题都答错，根据题目中的描述，我们知道答错题的概率是0.4，10道全错的概率就是0.4的10次方，即0.0001048576。

这个案例中描述的变量就是一个二项分布随机变量，每回答一道题就是一个试验，每次试验成功的概率是固定的，这里是0.6，这10次试验（答10道题）是互相独立的，在重复的独立试验中，如查每个试验仅可能为两种可能结果之一，而且得到相应结果的概率在每次试验中保持不变，则称这些重复独立试验为Bernoulli试验，一般用p和q=(1-p)分别表示每次试验得到的两种结果的概率，这两种结果通常被抽象地称为”成功“和”失败“，这当然不是实际意义上的成功和失败，仅仅是两种结果的代号而已。

人们一般会关心，如果进行n次Bernoulli试验，每次成功的概率是p，那么某中成功次数X=k次的概率是多少？

我们接着用下面的抽球案例来理解这个问题。

抽球实验

Bernoulli试验：一个罐子中，含有1个红球，3个蓝球，每次随机取出1个球，观察其颜色，再放回，再抽取下一个球（这是有放回的抽样），每次抽得的红球的概率为p=1/4，蓝球的概率是q=1-p=3/4，用R表示红球的事件，用B表示蓝球的事件，而用P(RB)表示第一次得到的是红球，第二次得到的蓝球的概率等，用X表示在n次试验中观测到红球的次数，下面以n=2和n=3为例说明一下这个概率。

其中下面的这个公式

mark

它表示从n个元素中选取k个元素的组合数。

现在看n=2的情况。

n=2表示是我们抽取2次，这里面又分3种情况。

第1种情况：2次抽取，抽到0个红球，也就是两次都是蓝球，这种情况的概率如下所示：

mark

第2种：2次抽取，刚好有一次观测到红球的概率（这种情况可以理解为，第1次红球，第2次蓝球或第1次蓝球，第2次红球），如下所示：

mark

第3种：2次抽取，均是红球的概率：

mark

现在看n=3的情况。

n=3表示是我们抽取3次，这里面又分4种情况。

第1种情况：抽取3次，只有0个红球，概率如下：

mark

第2种情况：抽取3次，只抽到1个红球，概率如下所示：

mark

第3种情况：抽取3次，抽到2个红球，概率如下所示：

mark

第4种情况：抽取3次，抽到3个红球，概率如下所示：

mark

从这个问题的描述中我们可以知道下面的信息：

mark

我们把在成功概率为p的n次Bernoulli试验中成功次数X称为具有参数(n,p)的二项分布(binomial distribution)的随机变量，可以用B(n,p)、Bin(n,p)或Binomial(n,p)等表示这样的二项分布，由于参数中的n和p可能取不同的值，因此二项分布是一族分布，称为二项分布族，在n=1时的二项分布又称为Bernoulli分布，它的概率分布为：P(X=1)=p，P(X=0)=1-p。

我们可以利用排列组合和二项式展开的性质证明：二项分布Binomial(n,p)变量的总体均值为μ=np，方差为δ^2=np(1-p)，标准差为δ=np(1-np)的平方根，即如下所示：

mark

下图就是二项分布Binomial(5,p)的图，这里p取0.1，0.2，，，一直到0.9，可以看出，当p=0.5时，分布是对称的，当p为0.1或0.9时，分布概率偏向左边或右边，如下所示：

k <>0.1,0.9,0.1)
par(mfrow=c(3,3))
for(i in 1:9)
  barplot(dbinom(0:5,5,k[i]),xlab='x',
          ylab='Probability',
          ylim=c(0,0.6),
          main=substitute(Binomial(5,b),
                          list(b=k[i])))


mark

案例分析

假定在一项外语考试中，如果每题有四个选项，只有1个是对的，一共有5道题，这些题目相互独立，如果答题者没有任何知识，全凭随机猜测，那么答对一道题的概率应该是0.25，此时答题就像是一系列Bernoulli试验，每次试验成功的概率是0.25，失败的概率为1-0.25=0.75，如果用X表示答对的题目，则X可以看成是参数为n=5，p=0.25的二项分布，下表就给出了对于k=0,1,2,3,4,5，X的概率分布函数P(X=k)和累积分布函数P(X≦ x)，如下所示：

mark

此时，有以下公式：

mark

下图的左图显示的是概率分布图，下图的右图是累积分布函数图，如下所示：

par(mfrow=c(1,2))
x <>0:5,5,0.25)
barplot(x,xlab='X',
        ylab='Probability Function',
        ylim=c(0,0.4),
        main='Binomial(5,0.25)')
plot(0:5,pbinom(0:5,5,0.25),
     ylim = c(0,1),
     type='s',
     xlab='X',
     ylab='Cumulative Distribution Function',
     main='Binomial(5,0.25)')
points(0:5,pbinom(0:5,5,0.25))

mark

多项分布

二项分布是基于每次的实验有两个可能结果的重复独立Bernoulli试验，如果把二项分布进行拓展，例如每次实验不是有2个结果，而有k个结果，并且针对每个结果的概率分别为(p1,p2,…pk)，并且这些概率都是固定的，其实可以简单地理解为掷骰子，这样在n次实验中，得到的k个可能的结果数目为X1,X2,…,Xk的概率分布为多项分布（multinomial distribution）。

多项分布可以用M(n,p1,…pk)表示，多项分布为多元分布，也就是说我们感兴趣的随机变量有k个(X1,X2,…,Xk)，这是一个k维随机向量，如果用P(X1=n1,X2=n2,…,Xk=nk)=p(n1,..,nk)代表多项分布k的可能结果在n次试验中分别出现p(n1,..,nk)次的概率，而p1,…,pk分别为一次试验时各种结果出现的概率，于是多项分布可以用下面的公式来表示：

mark

其中下面就是多项式系数：

mark

它就是多项式展开，如下所示：

mark

多面项分布的一个典型就是掷骰子，此时k=6，而对于一个正常的骰子来说，每掷一次，它的6个面朝上的概率是一样的，即p1=p2=… =p6=1/6，显然，多项分布是一个分布族，用不同的参数来区分。这里不再列出具体的使用案例。

超几何分布

超几何分布是一种离散概率分布，它描述的是不放回抽样的问题。例如在一批产品里，一共有N件产品，其中有M件次品，那么当我们任何取n件产品，其中恰有X件次品的概率P就可以按归照下面的公式进行计算：

mark

其实超几何分布就是不放回的Bernoulli试验，它也是GO分析的统计学基础。此处我们回忆一下二项分布，如下所示：

二项分布：一个罐子里有m个白球，n个黑球，如果每次随机取出1个，记录下颜色后放回，重复这样的试验就形成了一系列独立的Bernoulli试验，每次摸到白球的概率是p=m/(m=n)，在取了k次后，取到白球的总次数X就服从参数k和p的二项分布。

如果我们不放回白球，我们在取了k后后，取到白球的总次数X就服从超几何分布(hypergeometric distribution)。

R语言中与超几何分布有关的函数为phyper，dhyper，rhyper，qhyper。

其中dhyper的参数为dhyper(x, m, n, k)，这几个参数的含义如下所示：
x：从一个含有白球和黑球的盒子里不放回取出的白球数目。
m：白球的数量
n：黑球的数量
k：抽取的次数

超几何分布的计算

现在计算一下，有100个白球，400个黑球，不放回取50个球，有10个白球的概率，如下所示：

> dhyper(10,100,400,50,log = FALSE)
[1] 0.1474


现在计算一下，有100个白球，400个黑球，不放回取50个球，小于等于10个白球的概率，如下所示：

> phyper(10, 100, 400, 50)
[1] 0.5851

超几何分布的概率密度曲线

现在绘制超几何分布的概率密度曲线，如下所示：

参数说明：其中N表示的是总的球的数目，即N=m+n，这里设定的数目是500个，k表示抽取的次数，这里分别绘制了白球为100个抽取50次，白球为200个抽取60次，白球为300个抽取70次的情况，如下所示：

set.seed(1000)
x<>0,50)
y1<>100, 400, 50)
y2<>200, 300, 60)
y3<>300, 200, 70)
plot(x,y1)
plot(x,y1,col='red',ylab='density',
     xlim=c(0,55),xlab='',
     main='Hypergeometric distribution')
lines(x,y1,col='red')
points(x,y2,col='blue')
lines(x,y2,col='blue')
points(x,y3,col='orange')
lines(x,y3,col='orange')
legend('topright',legend=paste('N=',c(500,500,500),
                               'k=',c(50,60,70)),
       lwd=1,col=c('red','blue','orange'))


mark

超几何分布的概率积累分布曲线

现在绘制超几何分布的概率积累分布曲线，参数说明：其中N表示的是总的球的数目，即N=m+n，这里设定的数目是500个，k表示抽取的次数，这里分别绘制了白球为100个抽取50次，白球为200个抽取60次，白球为300个抽取70次的情况，如下所示：

set.seed(1000)
x<>0,50)
y1<>100, 400, 50)
y2<>200, 300, 60)
y3<>300, 200, 70)
plot(x,y1)
plot(x,y1,col='red',ylab='density',
     xlim=c(0,55),xlab='',
     main='Hypergeometric Cumulative distribution')
lines(x,y1,col='red')
points(x,y2,col='blue')
lines(x,y2,col='blue')
points(x,y3,col='orange')
lines(x,y3,col='orange')
legend('bottomright',legend=paste('N=',c(500,500,500),
                               'k=',c(50,60,70)),
       lwd=1,col=c('red','blue','orange'))

mark

案例分析

假定在一个20个人的研究生班中，有12个女生（用m表示，即m=12），有8个男生（用n表示，即n=8），因此，女生的比例是60%，如果随机抽取10个学生进行论文主审查，那么X=8篇女生的论文被抽中的概率是多少呢？

按照超几何分布Hyper(12,8,10)来计算，容易得出，其概率为0.0750，如果按照二项分布Binomial(10,0.6)来计算，则得出的概率是0.1209，显然这两种计算的差异很大，但是如果有120个女生和80个男生，那么随机抽取10个学生的论文，其中有8篇是女生的论文的概率就基本上差不多了，用超几何分布过计算，此概率是0.1183，用二项分布来计算，是0.1209。

现在计算一下这两个概率。

第1个概率（女生12个，男生8个），如下所示：

> m=12;n=8;k=10;x=8
> dhyper(x,m,n,k)
[1] 0.07501786
> dbinom(x,k,m/(m+n))
[1] 0.1209324


dhyper的参数前面了解了，现在看一下dbinom的参数，它的参数有3个，如下所示：

dbinom(x, size, prob, log = FALSE)

其中x表示目标事件的次数，我们的目标事件是”论文中是女生论文“，目标次数是8篇，此处x就等于8，k表示抽样的次数，这里是取了10个同学，prob表示目标事件的概率，我们的目标事件是”论文中是女生论文“，因此这个prob此处为0.6，也就是m/(m+n)。

现在扩大一下样本，女生120个，男生80个，抽10份论文，女生论文是8篇的概率，如下所示：

> m=120;n=80;k=10;x=8
> dhyper(x,m,n,k)
[1] 0.1182677
> dbinom(x,k,m/(m+n))
[1] 0.1209324


在这个案例中，我们研究的是具体的某个事件的概率，也就是抽到女生论文是8篇的概率，如果我们要研究至少投到8篇是女生论文的概率，那么就需要求出8，9，10篇是女生论文的概率之和（女生人数为120人，男生80个），如下所示：

> dhyper(8,120,80,10)+dhyper(9,120,80,10)+dhyper(10,120,80,10)
[1] 0.1606976

这里使用到了一个phyper函数，它的用法如下所示：

phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)


q就是目标事件的次数，这里是8，因此使用这个函数计算出来的结果其实就是，抽取10篇，其中小于等于8篇是女生论文的概率，因此这个概率并不是我们想要的，我们想要的概率是至少8篇是女生论文的概率，也就是8，9，10篇。因此，此处我们计算出小于等于7篇是女生论文的概率，再用1减去就是我们要想的概率，如下所示：

> 1-phyper(7,m,n,k)
[1] 0.1606976

现在用dhyper函数来验证下，分别计算出8，9，10篇是女生论文的概率，如下所示：

> dhyper(8,m,n,k)
[1] 0.1182677
> dhyper(9,m,n,k)
[1] 0.03726013
> dhyper(10,m,n,k)
[1] 0.005169843
> dhyper(8,m,n,k)+dhyper(9,m,n,k)+dhyper(10,m,n,k)
[1] 0.1606976


结果是一样的。

二项分布与超几何分布的联系与区别

在抽样调查的实践中，一般不会重复调查同一个人，这相当于不放回抽样。严格来说应该用超几何分布来描述，但对于很大的总体，通常都考虑二项分布模型。这是因为，在一个很大的一群人中进行随机调查，就相当于从有大量黑白球的罐子中摸球，无论是放回抽样，还是不放回抽样，抽样时待抽取的黑白球所点的比例不会有显著的区别，由于二项分布计算简单，所以都用二项分布来近似。但是如果仅仅在一个人数不多的单位中进行调查，就相当于上面的不放回摸球，而且由于球的总数不多，每次抽样后，待抽取的黑白球所占的比例变化会很大，因此就应该考虑用超几何分布。

Poisson分布

在二战中，德国的导弹不断袭击英国，为了研究导弹袭击的规律，英国的统计学家R.D. Clarke把一个地区划分为许多小块，使得每个小块被击中的概率极小，在袭击是随机的假定下，各个小块被击中的概率在整个区域假定是一个常数，这样，一个小块被南路的次数就大致相当于在大量重复的成功概率很小的独立Bernoulli试验中成功的次数，他于是就导出了Poisson分布，而实际观测的导弹袭击的频数和用Poisson分布预测的频率相当吻合，在许多事件中，例如某种疾病的发生数目，一个自动生产线出现的废品数，某种事故发生的数目等，只要出现一次的概率很小，一般被看成近似地服从Poisson分布，Poisson分布的总体参数是其均值，而且其总体均值和总体方差相等，这里用λ表示该总体均值，Poisson分布常用P(λ)或Poisson(λ)表示，此分布P(λ)的概率分布为：对于x=0,1,…，以及λ>0，有：

mark

下图就是参数为3、6、10的Poission分布图（只标出了20之内的部分），Poisson分布是一个分布族，其成员用不同的λ来区分，若干个独立的Poisson分布变量的和也服从Poisson分布，其参数等于那些Poisson分布的参数之和。

x<>0:20
z<>3),dpois(x,6),dpois(x,10))
matplot(x,z,
        pch=c(19,8,17),
        cex=1.5,
        xlab='x',ylab='P(X=x)',
        col=1,
        main='Poisson Distribution',
        type='b')
legend(15,0.23,
       c('P(3)','P(6)','P(10)'),
       cex=1.5,
       pch=c(19,8,17),
       col=1)

mark

正如前面的导弹袭击案例所说，当Bernoulli试验的次数很大，而每次成功概率很小时，可以用Poisson分布来近似二项分布，数学上，如果二项分布Binomial(n,p)的p趋于0，而np趋于常数λ时，该二项分布则趋于Poisson分布P(λ)，由于当n很大，以及p很小时，二项分布计算很麻烦，因此在计算机出现之前，人们在计算概率时，往往用P(np)来对二项分布Binomial(n,p)做近似。

看一个案例，现在使用不同的n（例如取值10，100，1000，10000，100000）和p（例如取值0.1，0.01，0.001，0.0001，0.00001）分别利用二项分布和Poisson分布近似来计算X=2的概率，如下所示：

> n<>10,100,1000,10000,100000)
> p<>0.1,0.01,0.001,0.0001,0.00001)
> cbind(dbinom(2,n,p),dpois(2,n*p))
          [,1]      [,2]
[1,] 0.1937102 0.1839397
[2,] 0.1848648 0.1839397
[3,] 0.1840317 0.1839397
[4,] 0.1839489 0.1839397
[5,] 0.1839406 0.1839397


参考资料

1. 统计学：从数据到结论.吴喜之.第四版