由于计算大阶乘时会遇到计算机上溢问题，一般的超几何分布表给的值也是参数小于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