防灾数学分享 http://blog.sciencenet.cn/u/fzmath 防灾科技学院数学教研室

博文

利用MATLAB生成超几何概率分布表

已有 1615 次阅读 2022-3-27 13:02 |系统分类:教学心得

由于计算大阶乘时会遇到计算机上溢问题,一般的超几何分布表给的值也是参数小于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://blog.sciencenet.cn/blog-292361-1331261.html

上一篇:利用MATLAB生成概率分布表xlsx文件
收藏 IP: 221.192.180.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-9-25 07:27

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部