|
由于计算大阶乘时会遇到计算机上溢问题,一般的超几何分布表给的值也是参数小于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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-12-3 15:57
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社