5

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

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 左右饱和。是由于数值近似还是我犯了任何错误？

4

5 回答 5

3

如果您采用“现有库”路线，您的编译器可能已经支持 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
于 2011-04-14T13:45:45.537 回答
2

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

于 2011-04-14T02:39:18.993 回答
2

由于您仅L在表达式中使用(p>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。

于 2011-04-14T09:21:33.897 回答
1

从我之前提出的另一个问题来看，您似乎也可以近似poisson(750)为poisson(375) + poisson(375)。

于 2011-04-14T12:06:32.850 回答
0

在这种情况下，您不需要多次调用随机数生成器。您只需要一张累积概率表： 

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]使用i > k递归关系 

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

并用于i < k使用

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

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

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

见：http ://en.wikipedia.org/wiki/Inverse_transform_sampling

于 2011-04-14T11:50:47.437 回答