由于计算大阶乘时会遇到计算机上溢问题，一般的超几何分布表给的值也是参数小于100的情形，限制了它的使用。这里参考石家庄军械学院王顺富的论文“超几何分布精确值的计算方法”，编写了计算的MATLAB程序，通过与自带函数的对比，说明程序正确有效。
代码如下：
function pmfk = hypgeopmf(k,n,D,N)
% hypgeopmf 超几何分布的分布律
% Y = hypgeopmf(k,n,D,N) returns the hypergeometric probability
% mass function at k with integer parameters n,D,N.
%%计算初值hypgeopmf(0,n,D,N)
d = max(0,n+D-N);
if d == 0
HdnDN = 1;
for j = 0:n-1
HdnDN = HdnDN*(N-D-j)/(N-j);
end
else%d>0
HdnDN = 1;
for j = d:n-1
HdnDN = HdnDN*(j+1)/(N-j+d);
end
end
pmfk = HdnDN; %初值
% 迭代公式hypgeopmf(k,n,D,N) = hypgeopmf(k-1,n,D,N)*(n-k+1)/(k*(N-D-n+k))
r = min(n,D);
if (k<=r)&&(k>=d)
for i = d+1:k
pmfk = pmfk*(D-i+1)*(n-i+1)/(i*(N-D-n+i));%公式(11)
end
else
pmfk = 0;%error('k is no more than min(n,D)');
end
为了检验程序是否正确，可以与MATLAB自带的函数hygepdf( )进行对比：
>> k = 20;n =80; D = 40; N = 200; pmfk = hypgeopmf(k,n,D,N),y = hygepdf(k,N,D,n)
pmfk =
0.0508
y =
0.0508
>> k = 3;n = 6; D = 4; N = 10; pmfk = hypgeopmf(k,n,D,N),y = hygepdf(k,N,D,n)
pmfk =
0.3810
y =
0.3810
可见，自编程序与MATLAB自带程序得到的结果相同。
又通过type hygepdf可以看到，内置函数用的是斯特林公式处理的，精度也很高。
>> format long
>> k = 20;n =80; D = 40000; N = 200000; pmfk = hypgeopmf(k,n,D,N),y = hygepdf(k,N,D,n)
pmfk =
0.056809627126298
y =
0.056809627126298
转载本文请联系原作者获取授权，同时请注明本文来自王福昌科学网博客。
链接地址：https://wap.sciencenet.cn/blog-292361-1331261.html?mobile=1
收藏
