第3章负超几何随机变量的两种抽样算法
在数理统计中，常把被考察对象的某一个（或多个）指标的全体称为总体。人们总是把总体看成一个具有分布的随机变量，总体中的每一个单元称为个体。把从总体中抽取的部分个体x1,x2,…,xn称为样本。随机变量有若干可能的取值，每个可能的取值有一定的可能性，根据随机变量遵循的规律使用一种方法获得它的具体值，叫做抽样。随机变量的抽样方法有多种，不同的分布采用的抽样方法不尽相同，但各种分布的随机抽样方法都是以服从［0,1］均匀分布的随机变量的抽样为基础产生的。
离散型随机变量的直接抽样是将随机数和阶梯型的随机变量的分布值逐项比较而确定相应的随机事件［66］。为了得到所需要的随机事件，直接抽样法所需要的比较次数也形成了另外一个随机变量，而且这个随机变量的数学期望和原随机变量的数学期望是一致的。这就形成了两个问题： 不确定的比较次数提高了抽样算法实现的复杂度。从抽样过程可以看出，离散型随机变量的直接抽样是一个“理想”的抽样方法，只能对分布律简单的随机变量有效。负超几何分布律中含有组合、阶乘运算，参数个数多，分布律函数复杂，所以用直接抽样法来生成负超几何随机变量不可行［67］，需要借助其他的分布，使用间接抽样法来生成。
本章利用乘抽样法，基于本书2.2节推导的改进的负二项近似，构造一种高效的负超几何随机变量抽样算法，分析抽样算法的效率，得到抽样效率为1.25，即平均生成1.25个［0,1］上均匀分布的样本，得到一个负超几何分布的样本。
此外，利用舍选抽样法，基于本书2.4节推导的伽马近似，构造了负超几何随机变量的一种精确抽样算法，严格证明抽样算法的正确性。
3.1负超几何随机变量的一种高效抽样算法
本章中，f,g,h表示概率分布，f(x)表示分布律或概率密度函数，χ表示集合，X,ξ表示随机变量，x表示随机变量在某一时刻的取值（即随机值），E表示抽样算法的效率，U(0,1)表示区间［0,1］上的均匀分布。
3.1.1乘抽样法
乘抽样法（Multiply Method）的基本思想是将一个复杂分布的抽样，转换成一个已知的、简单分布的抽样，它能克服直接抽样法实现困难的缺点，适用于较复杂的概率分布。抽样思想是通过一个容易生成的概率分布f1和一个取舍准则生成另一个与f1相近的概率分布f，具体来说分布f和f1同为集合χ上的分布且满足f(x)=f1(x)·h(x)，其中f1(x)的抽样方法已知，h(x)是x的非负函数且上确界为a，则可以通过对分布f1的样本使用某种机制舍弃其中的一些样本个体，保留剩下的样本，使得其服从分布f。具体步骤如下［68］：
（1） 生成f1的样本X；
（2） 生成ζ～U(0,1)，如果ζ≤h(X)/a，输出X； 否则转步骤（1）；
（3） 如能以任意精度计算ζ和ζ≤h(X)/a的值, 那么以上输出的X服从f分布，抽样效率E=1/a。
使用乘抽样法时，应注意如下两点：
（1） 为了提高乘抽样法的效率，a的值应该尽可能得小，也就是使分布f1和f分布更为相近，通常情况下选择f分布的近似分布作为f1；
（2） f1的高效抽样算法已知。
3.1.2算法构造和分析
本节首先找与负超几何分布相关的f1(x)和h(x)，然后计算h(x)的最大值。在第2章中，推导了负超几何的一种改进的负二项近似，它由负二项分布NBD(r,p)乘以一个修正因子得到。使用乘抽样法生成服从负超几何分布的随机变量，把NBD(r,p)当作f1(x)，把修正因子当作h(x)，即有
f1(x)=NBD(r,p)(3.1)
h(x)=1+x(x-1)2N-
12Npq(x-r)(x-r-1)p+r(r-1)q
(3.2)
接下来求函数h(x)的上确界。
定理3.1令h(x)=1+x(x-1)2N-12Npq［(x-r)(x-r-1)p+r(r-1)q］，当x取不小于rp+1的最小正整数时，函数h(x)取得最大值1.25，其中r≥1,x≥r,q=1-p。
证明：
由于
h(x)-h(x-1)
=x(x-1)2N-(x-r)(x-r-1)p2Npq-(x-1)(x-2)2N+
(x-r-1)(x-r-2)p2Npq
=x(x-1)-(x-1)(x-2)2N+
(x-r-1)(x-r-2)p-(x-r)(x-r-1)p2Npq
=2(x-1)2N-2(x-r-1)2Nq
=x-1N-x-r-1Nq≥0
以及
h(x+1)-h(x)=xN-x-rNq≤0
得到x≤rp+1及x≥rp， 而rp+1-rp=1, x是正整数，可知当x取不小于rp+1的最小正整数时，h(x)取得最大值。特别的，当rp为整数时，h(x)在rp和rp+1处同时取得最大值，其中r≥1,x≥r,q=1-p。且有
h(x)≤1+rp+1rp2N-rp+1-rrp+1-r-1p+r(r-1)q2Npq
=1+rp+1q-rp+1-r(1-p)-(r-1)q4q
=1+rp+1-rp+1-r-(r-1)4
=1+rp+1-rp-1+r-r+14
=1.25(3.3)
得证。
图3.1给出了负超几何随机变量乘抽样算法NHGMultiplySample的伪代码。
NHGMultiplySample (N,M,cc)
输入： N,M是整数
输出： x是整数
1. a=1.25, p=MN, q=1-p, r=M2
2. x←NB(r,p,cc)
3. ξ←U(0,1)
4. 假如 xN-M+r，转步骤2
5. H←1+x(x-1)2N-(x-r)(x-r-1)2Nq+r(r-1)2Np
6. 假如 ξ≤H/a，返回x
7. 转步骤2
图3.1NHGMultiplySample算法
利用NHGMultiplySample算法生成了100个负超几何随机值，如图3.2所示，其中N=1000,M=0.4,r=200。
498491498468494494488456528485505473495478
505501493497464486511506500486504499494489
494509522507504480492478506495517524503479
478505491498494512497501508505511501480524
498472520501475496464511507524504494489494
510467488477498489520501506521464500424436
507528533541552571436400399498501569577590
421400
图3.2NHGMultiplySample生成的100个随机数
表3.1给出了利用NHGMultiplySample算法分别生成10000个、1000个和100个随机变量时，时间和内存的消耗。算法在华硕星锐4752G上运行，2.9 GHZ Pentium(R) PC，C 编译器，算法使用的是文献［67］中构造的的抽样算法生成负二项随机变量，利用The GNU MP Bignum Library（GMP）库进行任意精度计算，其中N=1000,M=0.5,r=250。
表3.1N=1000,M=0.5,r=250时NHGMultiplySample算法的执行效率
100个随机值
1000个随机值
10000个随机值
时间耗费（s）
0.067934
0.395207
6.126978
内存耗费362 words
3.2负超几何随机变量的一种精确抽样算法
舍选法（AcceptanceRejection Method）是冯·诺依曼为克服直接抽样法的缺点而提出来的随机变量抽样法，它适用于概率密度函数复杂的分布。目前，负超几何概率的分布函数以及分布函数的反函数都未知，要对它的随机变量进行精确抽样，采用直接抽样法显然不行，必须采用间接抽样法。舍选法就是这样一种间接抽样法，它的抽样思想是为了实现从已知概率密度函数f(x)抽样，选取与f(x)取值范围相同的概率密度函数g(x)，生成服从g(x)的随机数序列{ζi},i=1,2,…,n，对{ζi}进行舍选，舍选的原则是在g(x)值大的地方，保留更多的随机数ζi； 在g(x)值小的地方，保留较少的随机数ζi，使得得到的子样中ζi的分布满足密度函数f(x)。一般来说，选择f的近似分布作为g分布。
因此利用舍选法对负超几何分布进行抽样，必须要找到与其相近的连续型分布，本书2.4节中提出的伽马分布就是符合要求的一种分布。
下面先描述舍选法的操作过程，然后详细介绍负超几何分布的一种高效抽样算法［70］，最后分析抽样算法的抽样效率并证明算法的正确性。
3.2.1舍选抽样法
假设f(x)和g(x)均为集合χ上的概率密度函数，且满足f(x)g(x)≤c,c≥0。x∈χ，舍选法的具体步骤如下［71］：
（1） 生成g的样本X；
（2） 生成ζ～U(0,1)，且ζ和X独立；
（3） 如果ζ≤f(X)/c·g(X)，则输出X（表示“接受”）； 否则转步骤（1）（表示“舍弃”）。
使用舍选抽样法时，应注意如下两点：
（1） f(X)和g(X)是互相独立的，因此f(X)/c·g(X)是和步骤（2）中的ζ是相互独立的；
（2） f(X)/c·g(X)的值在0到1之间，即0
用T表示“成功抽取一个个体”所需要执行的步骤（1）和（2）的次数，则T服从几何分布，即p=P(ζ≤f(X)/c·g(X))； P(T=n)=(1-p)n-1p,n≥1，那么执行步骤（1）和（2）的平均次数等于T的数学期望，μ(T）=1/p。计算p的值，则
pr=P(ζ≤f(X)/c·g(X)|X=x)
=f(x)/c·g(x)(3.4)
由于
pr=∫+∞-∞f(x)c·g(x)·g(x)dx
=1c∫+∞-∞f(x)dx
=1c(3.5)
所以μ(T）=1pr=c，即成功抽取一个个体所需要执行步骤（2）的次数为c次。显而易见，c值越小，抽样效率越高。为了使c值尽可能得小，选择f的近似分布作为g分布，且c=maxf(x)/g(x)。
3.2.2c值的计算
本章把伽马分布Ga(r,λ)作为超几何分布NHGD(r,N,M)的连续型近似。设计抽样算法时，要用到抽样效率c的值，从抽样过程可知,NHGD(r,N,M)/Ga(r,λ)的最大值是c，但直接求NHGD(r,N,M)/Ga(r,λ)的最大值不可行，此处借用负二项分布NBD(r,p)作为中间工具，先计算：
Q1←(NHGD(r,N,M)/NBD(r,p))max(3.6)
式（3.6）表示把NHGD(r,N,M)/NBD(r,p)的值上确界赋值给Q1，接着计算：
Q2←(NBD(r,p)/Ga（r,-ln(1-p)）max(3.7)
式（3.7）表示把NBD(r,p)/Ga（r,-ln(1-p)）的上确界赋值给Q2，最后将
c←Q1·Q2(3.8)
式（3.8）表示把Q1·Q2的值赋给c，通过此法得到值c。
定理3.2设N,M∈N,M≤N,r∈{1,2,…,M}，令p=M/N，那么
NHGD(r,N,M)NBD(r,p)
![《基于负超几何分布的十进制分组加密》[33M]百度网盘|亲测有效|pdf下载](/uploads/s0309/feeab2893e0681aa.jpg)