在c ++中生成泊松变量

| 我实现了此功能以生成泊松随机变量

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);
}

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

已邀请:

5 个回复

师埠女

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

佩疵瓦

如果使用“现有库”路径，则您的编译器可能已经支持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\';
}

我在上面使用了两种方法：我试图模仿您现有的界面。如果您创建带有平均值的std :: poisson_distribution，则以相同的平均值反复使用该分布会更有效（如main（）中所做的那样）。这是我的示例输出：

751
730
779

梆晨灸碾

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

瞧叮

由于只在表达式(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);

“ 11”只有11位的指数，但有52位的尾数。因此，这极大地增加了数值稳定性。付出的代价是每次迭代都需要log，而不是单笔exp。

夏瓤跋棘

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

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很大，这将变得非常不准确。但是，我们可以在这里使用一个技巧：以value21ѭ开始于（或接近）最大值，并假装p[k]等于1。然后使用递归关系为i > k计算p[i]

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

对于for27ѭ使用

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

这样可以确保最大的概率具有最大的精度。现在只需使用c[i+1] = c[i] + p[i+1]计算c[i]，直至c[i+1]与c[i]相同。然后，可以通过除以该极限值c[i]来归一化数组。或者您可以保留数组原样，并使用随机数number34ѭ。参见：http://en.wikipedia.org/wiki/Inverse_transform_sampling

要回复问题请先登录或注册

在c ++中生成泊松变量

5 个回复

发起人

c++

probability

poisson

knuth

问题状态

在c ++中生成泊松变量

与内容相关的链接

5 个回复

发起人

c++

probability

poisson

knuth

问题状态