生成泊松變量++

我實現這個功能，以生成泊松隨機變量生成泊松變量++

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是梅森旋轉算法的隨機數發生器。我發現，隨着我增加lambda，預期的分佈將會出現錯誤，平均值達到750左右。這是由於數值逼近還是我犯了錯誤？

來源

2011-04-14 Bob

IIRC，泊松變量具有指數分佈。因此，這是http://stackoverflow.com/questions/2106503/pseudorandom-number-generator-exponential-distribution的精確副本。但即使我錯了，那裏給出的方法也應該起作用。 – MSalters 2011-04-14 09:09:18

@ MSalters：泊松分佈是離散的 - 它只需要整數值。指數分佈是連續的。所以他們不一樣（雖然他們是相關的）。 – TonyK 2011-04-14 11:16:20

維基百科的權利：「如果給定時間間隔[0，t]內的到達數量遵循泊松分佈，平均值=λt，那麼到達間隔時間的長度將遵循指數分佈，平均值爲1/λ「。這是兩者之間的有效轉換，在結構上類似於我在下面提出的算法。 – MSalters 2011-04-14 11:57:12

從another question I asked earlier，似乎你也可以近似poisson(750)爲poisson(375) + poisson(375)。

來源

2011-04-14 12:06:32 MSalters

exp（-750）是一個非常小的數字，非常接近最小可能的雙倍，所以你的問題是數值。無論如何，你的複雜度在lambda中是線性的，所以對於高lambda來說算法不是很有效率。除非你有很好的理由來自己編寫代碼，否則使用現有的庫實現可能是有意義的，因爲這些數值算法對於你遇到的精確問題而言往往會很敏感。

來源

2011-04-14 02:39:18 bsdfish

我想我會使用正常的近似值，因爲在我的情況下，lambda總是一個很大的數值。 – Bob 2011-04-14 03:28:29

由於在表達式(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前面

來源

2011-04-14 09:21:33 MSalters

在這樣的情況下，你不需要多次調用隨機數發生器。 u需要是累積概率的一個表：

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+1] = (p[i]*lambda)/(i+1)

和i < k使用

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

這確保了最大的概率有最大可能的精確計算p[i]爲i > k。

現在只需使用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

來源

2011-04-14 11:50:47 TonyK

難道你不能存儲'log（p [k]）'而不是？這就是'（k log（λ））/（λ* log（k！））'，並且計算出來並不困難（參見http://en.wikipedia.org/wiki/Factorial#Rate_of_growth for log k！）'） – MSalters 2011-04-14 12:12:41

這是一個倒退的步驟。 log（k！）的精度隨着k的增加而下降，而我們希望最準確的值在k-lambda的平均值附近。另外，根本不需要記錄或者展示。 – TonyK 2011-04-14 12:24:44

如果你走「現有的庫」路線，你的編譯器可能已經支持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

來源

2011-04-14 13:45:45

生成泊松變量++

回答

相關問題