使用R中Stats包dmultinom计算概率返回0的问题求助

问题描述 #

使用R的stats包中dmultinom()函数计算真实生物数据的多项分布概率时，始终返回0，但随机生成的小数据集可正常运行。相关数据如下：

观测值向量：

observed_values <- c(23, 14, 81, 0, 0, 0, 154, 0, 0, 54, 0, 0, 0, 50, 174, 7, 1719, 0, 0, 3, 0, 0, 738, 28, 1653, 17, 77, 2, 44, 242, 59, 0, 0, 23, 108, 14, 268, 22, 0, 22, 25, 0, 0, 268, 213, 0, 377, 10, 0, 200, 91, 0, 0, 0, 1231, 0, 122, 0, 0, 462, 322, 0, 0, 0, 0, 24, 233, 0, 31, 382, 484, 82, 31, 81, 48, 0, 0, 0, 3, 0, 0, 93, 0, 0, 0, 4954, 111, 0, 0, 473)

对应概率向量：

probs <- c(0.0013982689, 0.0059828432, 0.0059828432, 0.0199428108, 0.0059828432, 0.0001747836,
           0.0518513081, 0.0008739181, 0.0358970594, 0.0179485297, 0.0119656865, 0.0059828432,
           0.0219370919, 0.0059828432, 0.0006991345, 0.0019942811, 0.0139599676, 0.0239313730,
           0.0039885622, 0.0001747836, 0.0219370919, 0.0159542486, 0.0119656865, 0.0039885622,
           0.0038452396, 0.0039885622, 0.0003495672, 0.0290595243, 0.0008739181, 0.0047191577,
           0.0013982689, 0.0059828432, 0.0059828432, 0.0199428108, 0.0059828432, 0.0001747836,
           0.0518513081, 0.0008739181, 0.0358970594, 0.0179485297, 0.0119656865, 0.0059828432,
           0.0219370919, 0.0059828432, 0.0006991345, 0.0019942811, 0.0139599676, 0.0239313730,
           0.0039885622, 0.0001747836, 0.0219370919, 0.0159542486, 0.0119656865, 0.0039885622,
           0.0038452396, 0.0039885622, 0.0003495672, 0.0290595243, 0.0008739181, 0.0047191577,
           0.0013982689, 0.0059828432, 0.0059828432, 0.0199428108, 0.0059828432, 0.0001747836,
           0.0518513081, 0.0008739181, 0.0358970594, 0.0179485297, 0.0119656865, 0.0059828432,
           0.0219370919, 0.0059828432, 0.0006991345, 0.0019942811, 0.0139599676, 0.0239313730,
           0.0039885622, 0.0001747836, 0.0219370919, 0.0159542486, 0.0119656865, 0.0039885622,
           0.0038452396, 0.0039885622, 0.0003495672, 0.0290595243, 0.0008739181, 0.0047191577)

原因分析 #

核心问题是数值下溢：

• 多项分布概率公式包含大整数阶乘（观测值总和为15947，15947!是一个极其庞大的数）和高次幂运算，这些计算结果远超出R中双精度浮点数的最小可表示范围（约1e-308）。
• 当概率值小于这个下限，R会直接返回0，导致计算失效。小数据集的阶乘和幂运算结果在浮点数范围内，因此能正常返回非零值。

解决方案 #

1. 使用对数概率计算 #

调用dmultinom()时设置log = TRUE，直接返回对数形式的概率，避免数值下溢：

log_prob <- dmultinom(x = observed_values, prob = probs, log = TRUE)

如果需要实际概率值，可通过exp(log_prob)转换，但需注意：若对数概率过小（如小于-700），exp()仍会返回0，但对数概率本身可用于似然比检验、模型拟合等后续分析，无需转换为原始概率。

2. 手动计算对数似然 #

多项分布的对数概率公式为：
$$\log P = \log(n!) - \sum_{i=1}^k \log(x_i!) + \sum_{i=1}^k x_i \log(p_i)$$
其中$n = \sum x_i$，$x_i$为观测值，$p_i$为对应概率。手动实现该公式：

n <- sum(observed_values)
log_fact_n <- lfactorial(n)
sum_log_fact_xi <- sum(lfactorial(observed_values))
sum_xi_log_pi <- sum(observed_values * log(probs))
log_prob_manual <- log_fact_n - sum_log_fact_xi + sum_xi_log_pi

该方法与dmultinom(..., log=TRUE)结果一致，且更直观。

替代包推荐 #

• VGAM：提供dmultinomial()函数，支持对数概率计算，同时针对离散分布提供更多扩展功能。
• statmod：包含针对大样本和数值稳定性优化的统计函数，可用于替代基础包中的分布计算。
• extraDistr：提供多种扩展分布函数，其中的多项分布实现也考虑了数值稳定性问题。

内容的提问来源于stack exchange，提问作者Hasinala Ramangason