超几何分布与相关变体

一、简介

【本文主题】在有限总体中进行不放回抽样常会遇到超几何分布，本文主要介绍超几何分布及其变体的基本公式与代码实现。

【本文定义】在本文中以下 4 个字母含义一致，粗体表示观测前的随机变量，普通字体表示随机变量的观测结果。

N：样本总量
M：阳性样本总量
n：抽样数量
m：抽到的阳性样本数量

【本文内容】

超几何分布：已知 N 和 M，抽取 n 个样本计算 m 的分布。
变体1：已知 N 和 n，通过测量 m 计算 M 的分布。
变体2：已知 M 和 n，通过测量 m 计算 N 的分布。
变体3：已知 N 和 M，通过测量 m 计算 n 的分布。

【后续案例】后文将会补充它们在实际场景中的应用案例。例如：

产品抽检合格的概率（超几何分布）。
一打刮刮乐中有多少张有奖的奖券（变体1）。
鱼塘中一共有多少鱼（变体2）。
一把曲别针有多少根（变体3）。

二、理论与推导过程

2.1 超几何分布：已知 N 和 M，抽取 n 个样本计算 m 的分布

【定义】盒子中共有 $N$ 个小球，其中有 $M$ 个红球，不放回地随机抽取 $n$ 个小球，抽到红球的个数 $\mathbf{m}$ 服从超几何分布，记为 $\mathbf{m}\sim h(n,N,M)$。

\[ P(\mathbf{m}=x \;|\;\mathbf{n}=n)=\dfrac{\dbinom{M}{x}\dbinom{N-M}{n-x}}{\dbinom{N}{n}},\quad x=0,1,2,\cdots,n. \]

理解：分母是从 $N$ 个小球中抽取 $n$ 个小球的组合数，分子是满足条件的组合数，即以下两种情况的乘积：①从 $M$ 个红球中抽取 $m$ 个红球，②从 $N-M$ 个其它颜色小球中抽取另外 $n-m$ 个小球的组合数。

【性质】交换 $n$ 和 $M$，概率分布不变：

$$ P(=x ;|;=n,=N,=M)=P(=

x;|;=M,=N,=n) $$

【数字特征】

\[ \mathbf{E}(\mathbf{m})=n\dfrac{M}{N} \]

\[ \mathbf{Var}(\mathbf{m})=\dfrac{nM(N-M)(N-n)}{N^2(N-1)} \]

【超几何分布的二项近似】

当 $n\ll N$ 时，即抽取个数远远小于盒中小球总数，每次不放回抽取后，红球占比 $p=M/N$ 几乎不变，这时不放回抽样可近似地看成放回抽样：

\[ \dfrac{\dbinom{M}{m}\dbinom{N-M}{n-m}}{\dbinom{N}{n}}\approx\dbinom{n}{m}p^{m}(1-p)^{n-m},\quad m=0,1,2,\cdots,n. \]

2.2 变体1：已知 N 和 n，通过测量 m 计算 M 的分布

【定义】盒子中共有 $N$ 个小球，其中红球个数 $M$ 待确定，不放回地随机抽取 $n$ 个小球，抽到了 $m$ 个红球，求盒中红球总数 $\mathbf{M}$ 的分布。

由于贝叶斯公式可得：

\[ \begin{align} P(\mathbf{M}=x \;|\; \mathbf{m}=m) &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{M}=x)P(\mathbf{M}=x) }{ \sum_\limits{i=0}^N P(\mathbf{m}=m \;|\; \mathbf{M}=i)P(\mathbf{M}=i) } \\ \\ &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{M}=x) }{ \sum_\limits{i=0}^N P(\mathbf{m}=m \;|\; \mathbf{M}=i) } \\ \\ &=\dfrac{ \dfrac{\dbinom{x}{m}\dbinom{N-x}{n-m}}{\dbinom{N}{n}} }{ \sum_\limits{i=0}^N \dfrac{\dbinom{i}{m}\dbinom{N-i}{n-m}}{\dbinom{N}{n}} } \\ \\ &=\dfrac{ \dbinom{x}{m}\dbinom{N-x}{n-m} }{ \sum_\limits{i=0}^N \dbinom{i}{m}\dbinom{N-i}{n-m} } \\ \\ &={\color{blue} \dfrac{ \dbinom{x}{m}\dbinom{N-x}{n-m} }{ \dbinom{N+1}{n+1} },\quad x=m,m+1,\cdots,N-(n-m). } \end{align} \]

注1：第二个等式，假设事先对 $\mathbf{M}$ 没有任何了解，即 $\mathbf{M}$ 的每种取值都等可能出现。

注2：最后一个等式，参见《统计趣题-前k个次序统计量的平均值的期望》第2.1节。

2.3 变体2：已知 M 和 n，通过测量 m 计算 N 的分布

【定义】盒子中小球总数 $N$ 待确定，其中有 $M$ 红球，不放回地随机抽取 $n$ 个小球，抽到了 $m$ 个红球，求盒中小球总数 $\mathbf{N}$ 的分布。

由于贝叶斯公式可得：

\[ \begin{align} P(\mathbf{N}=x \;|\; \mathbf{m}=m) &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{N}=x)P(\mathbf{N}=x) }{ \sum_\limits{i=N_0}^\infty P(\mathbf{m}=m \;|\; \mathbf{N}=i)P(\mathbf{N}=i) } \\ \\ &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{N}=x) }{ \sum_\limits{i=N_0}^\infty P(\mathbf{m}=m \;|\; \mathbf{N}=i) } \\ \\ &=\dfrac{ \dfrac{\dbinom{M}{m}\dbinom{x-M}{n-m}}{\dbinom{x}{n}} }{ \sum_\limits{i=N_0}^\infty \dfrac{\dbinom{M}{m}\dbinom{i-M}{n-m}}{\dbinom{i}{n}} } \\ \\ &={\color{blue} \dfrac{\dbinom{M}{m}\dbinom{x-M}{n-m}}{\dbinom{x}{n}}\dfrac{m(m-1)}{nM},\quad x=N_0,N_0+1,\cdots,\infty. } \end{align} \]

其中 $N_0=\mathbf{max}(M,n,m)$。

注1：第二个等式，假设事先对 $\mathbf{N}$ 没有任何了解，即 $\mathbf{N}$ 的每种取值都等可能出现。

注2：最后一个等式，用到了以下公式（证明候补）：

\[ \sum_\limits{i=\max(M,n,m)}^\infty \dfrac{\dbinom{M}{m}\dbinom{i-M}{n-m}}{\dbinom{i}{n}}=\dfrac{nM}{m(m-1)}. \]

2.4 变体3：已知 N 和 M，通过测量 m 计算 n 的分布

【定义】盒子中共有 $N$ 个小球，其中有 $M$ 个红球，随机抽取一批样本，其中有 $m$ 个红球，求抽取样本量 $\mathbf{n}$ 的分布。

由于贝叶斯公式可得：

\[ \begin{align} P(\mathbf{n}=x \;|\; \mathbf{m}=m) &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{n}=x)P(\mathbf{n}=x) }{ \sum_\limits{i=m}^N P(\mathbf{m}=m \;|\; \mathbf{n}=i)P(\mathbf{n}=i) } \\ \\ &=\dfrac{ P(\mathbf{m}=m \;|\; \mathbf{n}=x) }{ \sum_\limits{i=m}^N P(\mathbf{m}=m \;|\; \mathbf{n}=i) } \\ \\ &=\dfrac{ \dfrac{\dbinom{M}{m}\dbinom{N-M}{x-m}}{\dbinom{N}{x}} }{ \sum_\limits{i=m}^N \dfrac{\dbinom{M}{m}\dbinom{N-M}{i-m}}{\dbinom{N}{i}} } \\ \\ &={\color{blue} \dfrac{\dbinom{M}{m}\dbinom{N-M}{x-m}}{\dbinom{N}{x}}\dfrac{M+1}{N+1},\quad x=m,m+1,\cdots,N. } \end{align} \]

注1：第二个等式，假设事先对 $\mathbf{n}$ 没有任何了解，即 $\mathbf{n}$ 的每种取值都等可能出现。

注2：最后一个等式，用到了以下公式（证明候补）：

\[ \sum_\limits{i=m}^N\dfrac{\dbinom{M}{m}\dbinom{N-M}{i-m}}{\dbinom{N}{i}}=\dfrac{N+1}{M+1},\quad 0\leqslant m\leqslant M\leqslant N. \]

三、计算代码

3.1 最短置信区间（精确求解）

使用条件：分布律是单峰的。

定义函数：计算离散随机变量最短置信区间

(*计算离散随机变量最短置信区间*)
(*使用条件：分布律是单峰的*)

(*参数 pdfN：分布律-分子，以随机变量为参数*)
(*参数 pdfD：分布律-分母，固定值*)
(*参数 xE：随机变量的点估计值*)
(*参数 xMin：随机变量的最小值*)
(*参数 xMax：随机变量的最大值*)
(*参数 cl：置信水平，默认值0.95*)
Clear[fShortestCI];
fShortestCI[pdfN_,pdfD_Integer,xE_Integer,xMin_Integer,xMax_Integer,cl_:0.95]:=Block[{
ca=1-cl,
xL=xMin(*临时置信区间左侧位置*),
xR=xMax(*临时置信区间右侧位置*),
pL=0(*临时左侧单个x的分布律值*),
pR=0(*临时右侧单个x的分布律值*),
lrChange=""(*记录循环中增加的左侧或右侧*),
sumPDFL=0(*左尾分布律求和存储变量*),
sumPDFR=0(*右尾分布律求和存储变量*),
sumPDFLimit=Floor[pdfD(1-cl)](*双尾分布律求和上限*)
},
(*循环从两侧收缩：找到超过置信水平的临界点*)
While[sumPDFL+sumPDFR<sumPDFLimit,
If[xL>=xE&&xR<=xE,Break[]];
If[xL<xE&&xR<=xE,sumPDFL+=pdfN[xL];xL++;lrChange="L";Continue[];];
If[xL>=xE&&xR>xE,sumPDFR+=pdfN[xR];xR--;lrChange="R";Continue[];];

pL=pdfN[xL];
pR=pdfN[xR];
If[pL<=pR(*哪侧增加概率小，收缩哪侧*),
sumPDFL+=pL;xL++;lrChange="L";
,
sumPDFR+=pR;xR--;lrChange="R";
];
];
(*修正循环中可能多计算了一次的参数*)
If[lrChange=="L",xL--;sumPDFL-=pdfN[xL];
(*如果末次收缩的是左侧，且右侧单个概率大于左侧单个概率，且如果换成收缩右侧也不会超过sumPDFLimit，则换成收缩右侧*)
(*If[xL-1>=xMin&&xR-1>xE&&pdfN[xL-1]<pdfN[xR]&&pdfN[xR]-pdfN[xL-1]+sumPDFL+sumPDFR<sumPDFLimit,
xL--;xR--;sumPDFL-=pdfN[xL];sumPDFR+=pdfN[xR+1]];*)
];
If[lrChange=="R",xR++;sumPDFR-=pdfN[xR];
(*如果末次收缩的是右侧，且左侧单个概率大于右侧单个概率，且如果换成收缩左侧也不会超过sumPDFLimit，则换成收缩左侧*)
(*If[xR+1<=xMax&&xL+1<xE&&pdfN[xL]>pdfN[xR+1]&&pdfN[xL]-pdfN[xR+1]+sumPDFL+sumPDFR<sumPDFLimit,
xL++;xR++;sumPDFL+=pdfN[xL-1];sumPDFR-=pdfN[xR]];*)
];
(*返回结果*)
Return[<|
"xL"->xL(*置信区间下限*),
"xE"->xE(*点估计值*),
"xR"->xR(*置信区间上限*),
"xCI"->{xL,xR}(*置信区间*),
"sumPDFL"->sumPDFL(*置信区间外，左侧分布律之和，未除分母*),
"sumPDFR"->sumPDFR(*置信区间外，右侧分布律之和，未除分母*),
"sumPDF"->sumPDFL+sumPDFR(*置信区间外，两侧分布律之和，未除分母*),
"actualCL"->1-(sumPDFL+sumPDFR)/pdfD(*实际置信水平*),
"isShrinkEP"->If[xL!=xR&&xL-1>=xMin&&pdfN[xL-1]==pdfN[xR],{True,"最后一次收缩时两侧概率相同，优先收缩了左侧，以下两置信区间的实际置信水平相同：",{ {xL-1,xR-1},{xL,xR}}},{False,"只存在一种最短置信区间：",{ {xL,xR},{xL,xR}}}](*最后一次收缩，是否为两边等概率收缩*)
|>];
];

函数使用示例：

nN=20;(*已知：样品总量*)
nn=12;(*已知：抽样数量*)
nm=5;(*已知：抽到的次品量*)

cCL=0.95;(*置信水平 Confidence Level*)

pdfN[x_]:=Binomial[x,nm]Binomial[nN-x,nn-nm];
pdfD=Binomial[nN+1,nn+1];
xE=Round[nN nm/nn];
xMin=nm;
xMax=nN-nn+nm;
cl=cCL;

ci=fShortestCI[pdfN,pdfD,xE,xMin,xMax,cl]

ci["sumPDFL"]/pdfD//N
ci["sumPDFR"]/pdfD//N
ci["actualCL"]//N

3.2 等尾置信区间（精确求解）

定义函数：计算离散随机变量等尾置信区间

(*计算离散随机变量等尾置信区间*)

(*参数 pdfN：分布律-分子，以随机变量为参数*)
(*参数 pdfD：分布律-分母，固定值*)
(*参数 xE：随机变量的点估计值*)
(*参数 xMin：随机变量的最小值*)
(*参数 xMax：随机变量的最大值*)
(*参数 cl：置信水平，默认值0.95*)
Clear[fEqualCI];
fEqualCI[pdfN_,pdfD_Integer,xE_Integer,xMin_Integer,xMax_Integer,cl_:0.95]:=Block[{
ca=1-cl,
xL=xMin(*临时置信区间左侧位置*),
xR=xMax(*临时置信区间右侧位置*),
pL=0(*临时左侧单个x的分布律值*),
pR=0(*临时右侧单个x的分布律值*),
lrChange=""(*记录循环中增加的左侧或右侧*),
sumPDFL=0(*左尾分布律求和存储变量*),
sumPDFR=0(*右尾分布律求和存储变量*),
sumPDFLimit=Floor[pdfD(1-cl)](*双尾分布律求和上限*)
},
(*循环从两侧收缩：找到超过置信水平的临界点*)
While[True,
(*当收缩到均值时，该侧不再收缩*)
If[xL>=xE&&xR<=xE,Break[]];
If[xL<xE&&xR<=xE,
pL=pdfN[xL];
If[sumPDFL+sumPDFR+pL<sumPDFLimit,sumPDFL+=pL;xL++;lrChange="L";Continue[];]];
If[xL>=xE&&xR>xE,
pR=pdfN[xR];
If[sumPDFL+sumPDFR+pR<sumPDFLimit,sumPDFR+=pR;xR--;lrChange="R";Continue[];]];

pL=pdfN[xL];
pR=pdfN[xR];
(*当两侧都可以收缩时，优先收缩概率小且可以收缩的那一侧，相等时优选收缩左侧*)
If[sumPDFL+sumPDFR+pL<sumPDFLimit&&sumPDFL+sumPDFR+pR<sumPDFLimit,
If[sumPDFL<=sumPDFR,
sumPDFL+=pL;xL++;lrChange="L";Continue[];
,
sumPDFR+=pR;xR--;lrChange="R";Continue[];
]
];
(*当仅有一侧可以收缩时，收缩能收缩的那一侧*)
If[sumPDFL+sumPDFR+pL<sumPDFLimit&&sumPDFL+sumPDFR+pR>=sumPDFLimit,sumPDFL+=pL;xL++;lrChange="L";Continue[];];
If[sumPDFL+sumPDFR+pL>=sumPDFLimit&&sumPDFL+sumPDFR+pR<sumPDFLimit,sumPDFR+=pR;xR--;lrChange="R";Continue[];];
(*不能再收缩时，跳出循环*)
If[sumPDFL+sumPDFR+pL>=sumPDFLimit&&sumPDFL+sumPDFR+pR>=sumPDFLimit,Break[]];
];

(*返回结果*)
Return[<|
"xL"->xL(*置信区间下限*),
"xE"->xE(*点估计值*),
"xR"->xR(*置信区间上限*),
"xCI"->{xL,xR}(*置信区间*),
"sumPDFL"->sumPDFL(*置信区间外，左侧分布律之和，未除分母*),
"sumPDFR"->sumPDFR(*置信区间外，右侧分布律之和，未除分母*),
"sumPDF"->sumPDFL+sumPDFR(*置信区间外，两侧分布律之和，未除分母*),
"actualCL"->1-(sumPDFL+sumPDFR)/pdfD(*实际置信水平*),
"isShrinkEP"->If[xL!=xR&&xL-1>=xMin&&pdfN[xL-1]==pdfN[xR],{True,"最后一次收缩时两侧概率相同，优先收缩了左侧，以下两置信区间的实际置信水平相同：",{ {xL-1,xR-1},{xL,xR}}},{False,"只存在一种最短置信区间：",{ {xL,xR},{xL,xR}}}](*最后一次收缩，是否为两边等概率收缩*)
|>];
];

函数使用示例：

nN=100;(*已知：样本总量*)
nn=20;(*已知：抽样数量*)
nm=7;(*已知：抽到的阳性样本数量*)

cCL=0.95;(*置信水平 Confidence Level*)

pdfN[x_]:=Binomial[x,nm]Binomial[nN-x,nn-nm];
pdfD=Binomial[nN+1,nn+1];
xE=Round[nN nm/nn];
xMin=nm;
xMax=nN-nn+nm;
cl=cCL;
ci=fEqualCI[pdfN,pdfD,xE,xMin,xMax,cl]
ci["sumPDFL"]/pdfD//N
ci["sumPDFR"]/pdfD//N
ci["actualCL"]//N

3.3 等尾置信区间（中心极限定理，正态分布近似求解）

用正态分布近似计算总体阳性占比的置信区间。

(*用正态分布近似计算总体阳性占比的置信区间*)
nN=100;(*已知：样本总量*)
nn=20;(*已知：抽样数量*)
nm=7;(*已知：抽到的阳性样本数量*)

cCL=0.95;(*置信水平 Confidence Level*)

(*中间计算*)
z=InverseCDF[NormalDistribution[0,1],1-(1-cCL)/2];
a=nn+z^2;
b=-(2 nm+z^2);
c=nn (nm/nn)^2;

(*总体阳性占比的置信区间*)
r={(-b-Sqrt[b^2-4 a c])/(2 a),N[nm/nn],(-b+Sqrt[b^2-4 a c])/(2 a)};

(*总阳性样本量的置信区间*)
nM=r nN;
(*输出*)
Grid[{ {"指标","左置信区间","均值","右置信区间"},
Flatten[{"占比",PercentForm[#,{8,2},NumberPadding->{" ","0"}]&/@r}],
Flatten[{"规模",NumberForm[#,{8,2},NumberPadding->{" ","0"}]&/@nM}]},
Alignment->Center,ItemSize->8,Frame->All]

3.4 等尾置信区间（贝叶斯定理，Beta分布近似求解）

用Beta分布近似计算总体阳性占比的置信区间。

(*用Beta分布近似计算总体阳性占比的置信区间*)
nN=100;(*已知：样本总量*)
nn=20;(*已知：抽样数量*)
nm=7;(*已知：抽到的阳性样本数量*)

(*总体阳性占比的置信区间*)
r={InverseCDF[BetaDistribution[nm+1,nn-nm+1],(1-cCL)/2],N[nm/nn],InverseCDF[BetaDistribution[nm+1,nn-nm+1],1-(1-cCL)/2]};

(*总阳性样本量的置信区间*)
nM=r nN;

(*输出*)
Grid[{ {"指标","左置信区间","均值","右置信区间"},
Flatten[{"占比",PercentForm[#,{8,2},NumberPadding->{" ","0"}]&/@r}],
Flatten[{"规模",NumberForm[#,{8,2},NumberPadding->{" ","0"}]&/@nM}]},
Alignment->Center,ItemSize->8,Frame->All]