来自 Beta 分布的随机数，C++

Posted 2023-04-15

【中文标题】来自 Beta 分布的随机数，C++

【英文标题】：Random numbers from Beta distribution, C++

【发布时间】：2012-05-08 15:22:14

【问题描述】：

我用 C++ 编写了一个模拟程序，它根据特定的概率分布生成 (1,000,000)^2 个数字，然后对它们进行处理。到目前为止，我使用了指数、正态、伽玛、均匀和泊松分布。这是其中之一的代码：

#include <boost/random.hpp>

...main...

    srand(time(NULL)) ;
    seed = rand();
    boost::random::mt19937 igen(seed) ;
    boost::random::variate_generator<boost::random::mt19937, boost::random::normal_distribution<> >
    norm_dist(igen, boost::random::normal_distribution<>(mu,sigma)) ;

现在我需要为 Beta 发行版运行它。到目前为止，我所做的所有分发都需要 10-15 个小时。 Beta 发行版不在 boost/random 包中，所以我不得不使用 boost/math/distributions 包。我找到了this page on ***，它提出了一个解决方案。这是（复制粘贴）：

#include <boost/math/distributions.hpp> 
using namespace boost::math;  
double alpha, beta, randFromUnif;  
//parameters and the random value on (0,1) you drew  
beta_distribution<> dist(alpha, beta); 
double randFromDist = quantile(dist, randFromUnif); 

我复制了它并且它有效。我的模拟的运行时间估计是线性且可准确预测的。他们说这将运行 25 天。我看到两种可能性： 1. 提出的方法不如我以前用于其他发行版的方法 2. Beta 分布更难从中生成随机数

请记住，我对 C++ 编码的了解不足，所以我提出的问题可能很愚蠢。我迫不及待地要等一个月才能完成这个模拟，所以我能做些什么来改进它吗？也许使用我使用的初始方法并修改它以使用 boost/math/distributions 包？我什至不知道这是否可能。

另一条可能有用的信息是，对于我需要生成的所有 (1,000,000)^2 数字，参数都是相同的。我这么说是因为 Beta 发行版确实有一个讨厌的 PDF，也许参数固定的知识可以以某种方式用于简化流程？只是随机猜测。

【参考方案1】：

beta 分布与 gamma 分布有关。令 X 是从 Gamma(α,1) 中抽取的随机数，Y 是从 Gamma(β,1) 中抽取的随机数，其中 gamma 分布的第一个参数是形状参数。则 Z=X/(X+Y) 具有分布 Beta(α,β)。通过这种转换，它只需要两倍于您的 gamma 分布测试的时间。

注意：以上假设了伽马分布的最常见表示，Gamma(shape,scale)。请注意，伽马分布随机生成器的不同实现会因参数的含义和顺序而异。

【讨论】：

谢谢，我不知道那个属性。我还使用了伽玛分布只是指数之和的事实。您的建议将运行时间从 600 小时更改为 40 小时。谢谢！

使用标准符号，这不应该是 Gamma(α,1) 和 Gamma(β,1) 而不是 Gamma(1,α) 和 Gamma(1,β)。也就是说，α和β应该是形状参数，而不是尺度参数。

@BenRI - 有多种方法可以参数化伽马分布。使用更广泛的一个是gamma(shape_factor, scale_factor)，使用这种用法，你是对的。该比例因子也可以是比率或平均值。其他人将形状指定为第二个参数（但这显然没有被广泛使用）。我将更新我的答案以反映更常见的符号并表明 shape 参数是关键。

【参考方案2】：

如果您想要一个非常类似于 Beta 但具有非常简单的封闭形式逆 CDF 的分布，则值得考虑使用 Kumaraswamy 分布：

http://en.wikipedia.org/wiki/Kumaraswamy_distribution

当需要快速大量随机样本时，它被用作 Beta 分布的替代方案。

【参考方案3】：

尝试使用优化进行编译。使用标志 -O3 通常会加快速度。请参阅此post on optimisation flags 或this overview 了解更多详细信息。

【讨论】：

你应该先在评论中询问 OP 他是否真的在使用优化标志。