C++ gamma_distribution returns infinity

我最近寫了一些代碼，需要從伽馬分佈中抽取大量的代碼。我使用標準的gamma_distribution方法實現了這個功能，但偶爾會發現（一旦在藍色月亮中），這會返回值爲1.＃INF的原因不明。C++ gamma_distribution returns infinity

下面是一個最小的例子，我發現展示這個問題。在我的測試中，我發現問題通常發生在大約十億次迭代左右。

#include "stdafx.h" 
#include <iostream> 
#include <random> 

std::random_device rd; 
std::default_random_engine generator(rd()); 

int _tmain(int argc, _TCHAR* argv[]) 
{ 
    // create gamma distribution and random variable x 
    std::gamma_distribution<double> rgamma(2.0,1.0); 
    double x; 

    // loop through a large number of iterations 
    for (unsigned long int i=0; i<int(4e9); i++) { 

     // print update to console every million iterations 
     if ((i+1)%int(1e6)==0) 
      std::cout << "iteration: " << (i+1)/1e6 << " million\n"; 

     // draw new value of x from gamma distribution 
     x = rgamma(generator); 

     // if x==infinity then break 
     if ((1.0/x)==0) { 
      std::cout << "Error at iteration " << (i+1) << ": x=" << x << "\n"; 
      std::cin.get(); 
      exit(1); 
     } 
    } 
    // print message if reach end of loop 
    std::cout << "end\n"; 
    std::cin.get(); 
    return 0; 
}

我很想知道別人是否可以複製這個問題。我不知道這是相關的，但上面的程序在Visual Studio 2010中被編寫爲Win32應用程序，並在具有8核英特爾處理器的Windows計算機上運行。

目前，我已經通過捕獲無限值並使其成爲大數來修補此問題。但如果任何人有任何洞察，爲什麼/如何發生這將不勝感激！

來源

2016-01-24 bobverity

據我所見，沒有禁止無限的規則。我會看看概率是否匹配。看看'rgamma.max（）'返回的是什麼。根據cplusplus「numeric_limits :: max（）或numeric_limits :: infinity（）」。這是實現定義的，因此其他人無法重現它 – WorldSEnder

我剛剛運行你的程序（在64位GNU/Linux上，用g ++編譯）進行了21.47億次迭代，而沒有停止在你的'exit（1）'處。然後它打印出'end'並正常退出（0）。這是隨機的，也許它可能發生了......但我沒有見證它。 – e0k

使用您的參數，伽馬函數變爲$ xe^{ - x} $，並將導致積分到$ e^{ - x} *（ - x-1）$，必須在'numeric_limits :: max（） '以產生你的回報價值更大的可能性。請注意，這是非常系統相關的，所以你必須在你的系統上執行它 – WorldSEnder

我所做的隨機種子確定型有：

//std::random_device rd; 
std::default_random_engine generator(-246744094);

，並在832million迭代可以一致地重現此錯誤。我用不同的編譯器（VS2013 32/64位和Intel C++ 2016 64位）以及不同的機器（我的桌面i7和羣集節點Xeon E5）獲得了相同的結果。我的rgamma.max（）也和Bob一致。

然後我用MinGW（g ++ main.cpp -o main.exe -std = C++ 0x -O3）編譯它 - 它運行時沒有生成任何infinites。

SO - 我推測：無論使用MS還是英特爾編譯器，它都是通過Visual Studio鏈接的庫中的某種錯誤。不確定是否可以報告/修復 - 或者是否可以更快地找到一些可靠的gammas源代碼。

來源

2016-01-24 20:53:40

C++ gamma_distribution returns infinity

回答

相關問題