17-SAS 与负二项分布

作者

Simon Zhou

发布于

2025年6月9日

修改于

2025年6月9日

代码 
%load_ext saspy.sas_magic

1 负二项分布中的参数估计

负二项分布( negative binomial distribution )是一种离散型分布,常用于描述生物的群聚性:在毒理学的显性致死试验或致癌试验中也都有应用。

二项分布中的 n 是固定的,当 n 不固定,并用 x+k 来替换 n 后,所得到的在 x+k 次试验中得到此种结果恰为 k 次的概率,这时的概率函数就是负二项分布,所以在是负二项分布中的一个重要的参数。

计算参数k的常用方法有动差法、频数法、零频数法、最大似然法等。这里介绍相对较为简单的动差法。

1.1 动差法示例

在研究某种毒物的致死作用时,对 60 只小白鼠进行了显性致死试验,得到数据资料见表。若该样本计数服从负二项分布,试估计其参数 μ 和 k。

胚胎死亡数 0 1 2 3 4 5 6 合计
观察雌鼠数 30 14 8 4 2 0 2 60

1.1.1 动差法原理（Method of Moments）

动差法是一种参数估计方法，通过令样本矩（如样本均值、样本方差）等于理论分布的对应矩，解方程得到参数估计。

对于负二项分布：

  • 参数：

    • \(\mu = \text{E}(X)\)
    • \(\text{Var}(X) = \mu + \mu^2/k\)

由样本数据可得样本均值 \(\bar{x}\) 和样本方差 \(s^2\)，则可解出：

  • \(\hat{k} = \mu^2 / (s^2 - \mu)\)

这个估计在 \(s^2 > \mu\) 时成立（即数据存在过度离散）。

1.1.2 其他参数估计方法

方法 思路 优点 缺点
频数法（Frequency Method） 利用频数分布构造估计量，如通过最大频数位置反推参数 直观，适用于整数型数据 精度差，参数间依赖强
零频数法（Zero Frequency Method） 利用观测中“零”的频率推断参数 简便，只需零频数 精度有限，需大量样本支持
最大似然估计（MLE） 构造似然函数，以数值优化求得参数估计 统计效率高，常用于建模 需要迭代计算，依赖软件
代码 
%%SAS
/*程序17-1*/
data nbinomial;
    input x f @@;
datalines;
0 30 1 14 2 8 3 4 4 2 5 0 6 2
;
run;
proc univariate;
    var x;
    freq f;
output out= mv2 mean = mu var = v;
run;
data k;
    set mv2;
    k = mu**2/(v-mu);
proc print;
    var mu k;
run;
SAS 输出

SAS 系统

UNIVARIATE 过程

变量: x

 

频数: f

矩
数目 60 权重总和 60
均值 1.03333333 观测总和 62
标准差 1.43759058 方差 2.06666667
偏度 1.78111198 峰度 3.3122114
未校平方和 186 校正平方和 121.933333
变异系数 139.121669 标准误差均值 0.18559215
基本统计测度
位置 变异性
均值 1.033333 标准差 1.43759
中位数 0.500000 方差 2.06667
众数 0.000000 极差 6.00000
    四分位间距 2.00000
位置检验: Mu0=0
检验 统计量 p 值
Student t t 5.567764 Pr > |t| <.0001
符号 M 15 Pr >= |M| <.0001
符号秩 S 232.5 Pr >= |S| <.0001
分位数（定义 5）
水平 分位数
100% 最大值 6.0
99% 6.0
95% 4.0
90% 3.0
75% Q3 2.0
50% 中位数 0.5
25% Q1 0.0
10% 0.0
5% 0.0
1% 0.0
0% 最小值 0.0
极值观测
最小值 最大值
值 频数 观测 值 频数 观测
0 30 1 1 14 2
1 14 2 2 8 3
2 8 3 3 4 4
3 4 4 4 2 5
4 2 5 6 2 7

SAS 系统

观测 mu k
1 1.03333 1.03333

1.1.3 程序说明

数据集中的x和f分别表示胚胎死亡数和雌鼠数,首先通过 univariate 过程计算均数和方差,并将该两项指标输出到 mv2 数据集中,再用数据集k调用mv2的内容,用专用公式计算k的值。

1.1.4 结果说明

univariate 过程的输出结果不再叙述,最后输出的两个参数分别为u=1.033 33,k=1.033 33.

1.2 零频数法

理论上，

\[ P(X=0) = \left( \frac{k}{k+\mu} \right)^k \]

设观察零频率 \(f_0 = 30/60 = 0.5\)，解此方程估计 k。

代码 
%%SAS
data zerofreq;
    f0 = 0.5;
    mu = 1.0;  /* 可以先用动差法得出的均值 */
    do k = 0.01 to 20 by 0.01;
        p0 = (k/(k+mu))**k;
        diff = abs(p0 - f0);
        output;
    end;
run;

proc sort data=zerofreq;
    by diff;
run;

proc print data=zerofreq(obs=1);
    var mu k p0;
run;
SAS 输出

SAS 系统

观测 mu k p0
1 1 1 0.5