Rust-Random项目中的超几何分布采样问题分析与修复

问题背景

在Rust生态系统中，rust-random/rand是一个广泛使用的随机数生成库。最近，开发团队发现其超几何分布(Hypergeometric distribution)采样实现存在严重问题。具体表现为，当使用参数Hypergeometric::new(100,50,49)时，生成的样本明显不符合预期的概率分布。

问题现象

通过卡方检验(Chi-squared test)可以清楚地观察到问题：对于100万次采样，得到的频率分布与理论分布存在显著差异。实际采样结果中，某些区间的样本数量为0，而理论上这些区间应该有非零概率。这种偏差导致统计检验的p值为0，表明采样结果极不可能来自预期的超几何分布。

问题根源分析

深入调查发现，问题主要存在于拒绝-接受采样(Rejection-Acceptance)方法的实现中。具体有两个关键问题：

1. 算法实现错误：在计算接受概率时，分母部分的计算顺序错误。原始代码中使用了错误的公式，导致接受概率计算不准确。这与R语言多年前发现的bug类似。

2. 阶乘对数近似精度不足：实现中使用了斯特林公式(Stirling's approximation)来计算阶乘的对数值，这对于小数值的近似精度不足。例如，对于ln(7!)，近似结果为6.62，而实际值应为8.53。

解决方案

针对上述问题，开发团队提出了以下修复方案：

1. 修正接受概率计算公式：将原来的计算方式改为正确的形式，确保分子分母的顺序正确。具体修改是将分母部分的计算调整为(n1 - i + 1) * (k - i + 1)。

2. 关于阶乘近似的考量：虽然斯特林公式的近似误差确实存在，但初步测试表明，在修正了主要算法错误后，这种近似对整体结果的影响相对较小。不过，对于需要高精度计算的场景，可以考虑使用更精确的阶乘计算方法。

影响范围

这个问题不仅影响(100,50,49)这一组参数，还会影响其他使用拒绝-接受采样方法的参数组合，例如(65,30,28)、(48,25,20)和(40,20,19)等。相比之下，使用逆变换采样(Inverse-transform sampling)方法的参数组合则表现正常。

技术启示

这个案例给我们几个重要的技术启示：

1. 算法移植需谨慎：当从其他语言或论文中移植算法时，需要特别注意实现细节，即使是微小的差异也可能导致完全错误的结果。

2. 统计测试的重要性：像卡方检验这样的统计测试是验证随机数生成器正确性的有效工具，应当在开发过程中充分使用。

3. 数值稳定性考量：在概率计算中，特别是涉及大数阶乘时，对数变换和近似方法的选择需要慎重考虑精度与性能的平衡。

总结

rust-random/rand库中的超几何分布采样问题展示了随机数生成算法实现中的潜在陷阱。通过深入分析问题根源并实施针对性修复，不仅解决了当前问题，也为未来类似问题的预防和处理提供了宝贵经验。对于使用者而言，升级到修复后的版本将确保超几何分布采样的正确性。