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

回答

相关问题