R统计学(03): 超几何分布
这篇博客将介绍另一个离散分布：超几何分布(Hypergeometric distribution)，其自变量X定义为从N个有限物品中抽出n个物品，成功抽出指定种类的物品的个数。
对于二项分布和多项分布，它们均基于伯努利试验，每次试验结果的发生概率是不变的，而超几何分布试验结果的概率会随着每一次试验的发生而改变。以随机抽样为例，二项分布试验和几何分布试验是有放回抽样(总体数量保持不变)，因此每次试验结果的发生概率是保持不变的；而超几何分布试验则是在有限总体中进行无放回抽样(总体数量不断减少)，所以每次试验结果发生的概率将发生变化。
1. 概率分布
超几何分布是一种重要的离散型概率分布，其概率质量函数可以这样定义：假设有限总体包含N个样本，其中质量合格的为m个，则剩余的N-m个为不合格样本，如果从该有限总体中抽取出n个样本，其中有k个是质量合格的概率为：
其中表示从N个总体样本中抽取n个样本的方法数目；
表示从m个质量合格样本中抽取k个样本的方法数目；
表示从N-m个质量不合格样本中抽取n-k个样本的方法数目。
由上式可知，超几何分布由样本总量N、质量合格的样本数m和抽取数目n决定，记为X~H(N, m, n)。
2. 性质
- k的取值范围为{max(0, n+m-N), ..., min(n, m) }，其期望值和方差分别为：
具体可参考https://en.wikipedia.org/wiki/Hypergeometric_distribution。
-
如果抽取数目n=1(即从有限总体中只抽取一个样本)，那么超几何分布将还原成伯努利分布。
-
如果样本总体N为无穷大(也即将有限总体换成无限总体)，此时是否放回抽中的样本对于质量合格样本所占比例没有影响，超几何分布也可视为二项分布。在实际应用中，只要样本总体的数目是抽取数目的10倍以上(即N>10n)，就可用二项分布来近似描述超几何分布，通过两种概率质量函数计算得到的概率几乎相同。对于这个性质，我们将在后面用一个例子来说明。
3. R中的相关函数
R中也有四个函数可用于超几何分布，分别是：
-
dhyper(x, m, n, k)：返回抽中x个质量合格样本的概率 -
phyper(q, m, n, k)：返回累积概率 -
qhyper(p, m, n, k)：返回相应分位点x，详情见下面的例子 -
rhyper(nn, m, n, k)：返回每组抽中质量合格样本的个数
这四个函数都有m、n和k参数，分别对应于超几何分布中的质量合格数目m、不合格数目N-m和抽取数目n。下面通过一个例子来了解如何使用它们：
假设某服装店举行十一促销抽奖活动，抽奖箱中总共有30个乒乓球，其中只有3个乒乓球上写有“中奖”两字。结账时，每个顾客抽出2个乒乓球。如果抽中一个中奖字样，商品总价打7折；如果抽中两个，商品总价打5折；如果没抽中就不打折。
第一个问题：抽到0个，1个和2个带有“中奖”字样乒乓球的概率分别是多少？此时要用到dhyper(x, m, n, k)函数，其中x参数指定抽中的数目，函数返回相应概率，比如：
> m <- 3 ##带有“中奖”字样乒乓球的数目
> n <- 30-m ##没有带有“中奖”字样乒乓球的数目
> k <- 2 ##抽取的数目
> dhyper(0:2, m, n, k)
[1] 0.806896552 0.186206897 0.006896552
从计算结果可知，无法中奖的概率竟然高达80.69%，而中奖的概率仅有19.31%，可见该服装店是不希望顾客中奖的。
我们上面提到，如果总体数目是抽取数目的十倍以上，可以用二项分布的概率质量函数近似超几何分布的概率质量函数。在本案例中，总体数目为30，是抽取数目2的15倍，尝试用二项分布的概率质量函数来计算顾客中奖的概率，中奖乒乓球的比例是10%(即p=0.1)，随机抽取2个乒乓球，二项分布概率质量函数的计算结果如下：
> p <- 0.1
> dbinom(0:2, 2, p)
[1] 0.81 0.18 0.01
结果表明，两套概率结果的数值非常接近，验证了上述的第三个性质。
第二个问题：至多抽到1个带有“中奖”字样乒乓球的概率是多少？此时要用到phyper(q, m, n, k)函数，其中q参数指定至多抽中的数目（这里为1），函数返回相应累积概率，比如：
> m <- 3 ##带有“中奖”字样乒乓球的数目
> n <- 30-m ##没有带有“中奖”字样乒乓球的数目
> k <- 2 ##抽取的数目
> phyper(1, m, n, k)
[1] 0.9931034
结果表明，至多抽到1个带有“中奖”字样乒乓球的概率高达99.3%
第三个问题：90%概率下我们至多能抽到几个带有“中奖”字样乒乓球的概率是多少？此时要用到qhyper(p, m, n, k)函数，其中p参数指定概率(这里是0.9)，函数返回相应分位点x(即F(x)≥0.9对应的最小x值)，比如：
> m <- 3 ##带有“中奖”字样乒乓球的数目
> n <- 30-m ##没有带有“中奖”字样乒乓球的数目
> k <- 2 ##抽取的数目
> qhyper(0.9, m, n, k)
[1] 1
结果表明，90%概率下我们至多能抽中1个
最后一个问题：重复10000组抽奖，每组抽中的个数是多少？这时就要用到rhyper(nn, m, n, k)函数，其中nn参数指定抽奖组数(这里为10000)，函数返回每组抽中的个数，比如：
> m <- 3 ##带有“中奖”字样乒乓球的数目
> n <- 30-m ##没有带有“中奖”字样乒乓球的数目
> k <- 2 ##抽取的数目
> set.seed(123)
> ns <- rhyper(10000, m, n, k)
> table(ns)
ns
0 1 2
8117 1821 62
> mean(ns) ##抽中数目的平均值
[1] 0.1945
> k*m/(m+n) ##理论值
[1] 0.2
> var(ns) ##抽中数目的方差
[1] 0.1690867
> k*m/(m+n)*n/(m+n)*(m+n-k)/(m+n-1) ##理论值
[1] 0.1737931
模拟1万位顾客抽奖，高达8116位顾客没有中奖，与理论上80.69%不中奖很接近。此外均值和方差也与理论值接近。
感谢您的阅读！想了解更多有关技巧，请关注我的微信公众号“R语言和Python学堂”，我将定期更新相关文章。
