在 c++ 中生成泊松变量

Question

我实现了这个函数来生成泊松随机变量，

typedef long unsigned int luint;
luint poisson(luint lambda) {
    double L = exp(-double(lambda));
    luint k = 0;
    double p = 1;
    do {
        k++;
        p *= mrand.rand();
    } while( p > L);
    return (k-1);
}

其中 mrand 是 MersenneTwister 随机数生成器。我发现，当我增加 lambda 时，预期分布将会出错，均值在 750 左右饱和。这是由于数值近似还是我犯了任何错误？

Answer 1

如果您选择“现有库”路线，您的编译器可能已经支持 C++11 std::random 包。以下是您如何使用它：

#include <random>
#include <ctime>
#include <iostream>

std::mt19937 mrand(std::time(0));  // seed however you want

typedef long unsigned int luint;

luint poisson(luint lambda)
{
    std::poisson_distribution<luint> d(lambda);
    return d(mrand);
}

int main()
{
    std::cout << poisson(750) << '\n';
    std::poisson_distribution<luint> d(750);
    std::cout << d(mrand) << '\n';
    std::cout << d(mrand) << '\n';
}

我在上面使用了两种方式：

1. 我试图模仿您现有的界面。

2. 如果您创建一个带有均值的 std::poisson_distribution，那么反复使用该分布来获得相同的均值（如 main() 中所做的那样）会更有效。

这是我的输出示例：

751
730
779

Answer 2

exp(-750) 是一个非常小的数字，非常接近最小可能的双精度数，因此您的问题是数值问题。无论如何，您的复杂性在 lambda 中都是线性的，因此该算法对于高 lambda 来说不是很有效。除非您有充分的理由自己编写此代码，否则使用现有的库实现可能是有意义的，因为这些数值算法往往对于您遇到的精度问题非常敏感。

Answer 3

由于您仅在表达式 (p>L) 中使用 L，因此您实际上是在测试 (log(p) > -lambda)。这不是一个很有帮助的转变。当然，您不再需要 exp(-750)，但您只会溢出 p 。

现在，p 就是 Π(mrand.rand())，log(p) 就是 log(Π(mrand.rand())) 就是 Σ(log(mrand.rand())。这为您提供了必要的转换：

double logp = 0;
do {
    k++;
    logp += log(mrand.rand());
} while( logp > -lambda);

double 只有 11 位指数，但有 52 位尾数，因此，这是数值稳定性的巨大提高，代价是您需要一个 log。 在每次迭代中，而不是前面的单个 exp 。

Answer 4

来自 另一个问题我之前问过，似乎您也可以将 poisson(750) 近似为 poisson(375) + poisson(375)。

Answer 5

在此类情况下，您无需多次调用随机数生成器。您所需要的只是一个累积概率表：

double c[k] = // the probability that X <= k (k = 0,...)

然后生成一个随机数 0 <= r < 1，并取第一个整数X，使得c[X]> r。您可以通过二分搜索找到这个X。

为了生成这个表，我们需要单独的概率。

p[k] = lambda^k / (k! e^lambda) // // the probability that X = k

如果 lambda 很大，这会变得非常不准确，正如您所发现的。但我们可以在这里使用一个技巧：从（或接近）最大值开始，使用 k = Floor[lambda]，并暂时假装 p[k]等于1。然后计算 p[i] for i > k 使用递推关系

p[i+1] = (p[i]*lambda) / (i+1)

并且对于 i k using

p[i-1] = (p[i]*i)/lambda

这可确保最大的概率具有最大可能的精度。

现在只需使用 c[i+1] = c[i] + p[i+1] 计算 c[i]，直到 c [i+1] 与 c[i] 相同。然后你可以通过除以这个限制值c[i]来标准化数组；或者您可以保留数组不变，并使用随机数 0 <= r < c[i]。

请参阅：http://en.wikipedia.org/wiki/Inverse_transform_sampling