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乱数ジェネレーターです。ラムダを増やすと、期待される分布が間違ってしまい、平均が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';
}

上記の2つの方法で使用しました。

  1. 私はあなたの既存のインターフェースを模倣しようとしました。

  2. 平均を使用してstd::poisson_distributionを作成する場合、同じ平均に対してその分布を繰り返し使用する方が効率的です（main（）で行われるように）。

これが私のためのサンプル出力です：

751
730
779
于 2011-04-14T13:45:45.537 に答える
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 に答える
2

exp(-750) は非常に小さい数値であり、可能な限り最小の double に非常に近いため、問題は数値です。いずれにせよ、複雑さはラムダで線形になるため、ラムダが高い場合、アルゴリズムはあまり効率的ではありません。これを自分でコーディングする大きな理由がない限り、既存のライブラリ実装を使用することはおそらく理にかなっています。これらの数値アルゴリズムは、遭遇している精度の問題に対して正確に扱いにくい傾向があるためです。

于 2011-04-14T02:39:18.993 に答える
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最初の整数を取ります。これは二分探索で 見つけることができます。Xc[X] > rX

このテーブルを生成するには、個々の確率が必要です

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

が大きい場合lambda、あなたが見つけたように、これは非常に不正確になります。ただし、ここでトリックを使用できます。最大値（またはその近く）から開始し、で、と等しいk = floor[lambda]瞬間のふりをします。次に、漸化式 を使用して計算しますp[k]1p[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 に答える