Beta分配的隨機數發生器

Question

我需要一個像betarand(a,b)這樣的函數的c或C++源代碼，它可以產生帶有beta分佈的隨機數。我知道我可以使用boost庫，但是我將把它移植到CUDA架構中，所以我需要代碼。有人能幫助我嗎？
同時我有betapdf（Beta概率密度函數）。但我不知道如何使用它來創建隨機數字:)。

Answer 1

的C++ 11的隨機數庫不提供Beta分佈。但是，Beta分佈可以用兩個伽馬分佈來建模，其中庫確實提供了。按照std::gamma_distribution爲您實施beta_distribution。據我所知，它完全符合隨機數分佈的要求。

#include <iostream> 
#include <sstream> 
#include <string> 
#include <random> 

namespace sftrabbit { 

    template <typename RealType = double> 
    class beta_distribution 
    { 
    public: 
     typedef RealType result_type; 

     class param_type 
     { 
     public: 
      typedef beta_distribution distribution_type; 

      explicit param_type(RealType a = 2.0, RealType b = 2.0) 
      : a_param(a), b_param(b) { } 

      RealType a() const { return a_param; } 
      RealType b() const { return b_param; } 

      bool operator==(const param_type& other) const 
      { 
      return (a_param == other.a_param && 
        b_param == other.b_param); 
      } 

      bool operator!=(const param_type& other) const 
      { 
      return !(*this == other); 
      } 

     private: 
      RealType a_param, b_param; 
     }; 

     explicit beta_distribution(RealType a = 2.0, RealType b = 2.0) 
     : a_gamma(a), b_gamma(b) { } 
     explicit beta_distribution(const param_type& param) 
     : a_gamma(param.a()), b_gamma(param.b()) { } 

     void reset() { } 

     param_type param() const 
     { 
     return param_type(a(), b()); 
     } 

     void param(const param_type& param) 
     { 
     a_gamma = gamma_dist_type(param.a()); 
     b_gamma = gamma_dist_type(param.b()); 
     } 

     template <typename URNG> 
     result_type operator()(URNG& engine) 
     { 
     return generate(engine, a_gamma, b_gamma); 
     } 

     template <typename URNG> 
     result_type operator()(URNG& engine, const param_type& param) 
     { 
     gamma_dist_type a_param_gamma(param.a()), 
         b_param_gamma(param.b()); 
     return generate(engine, a_param_gamma, b_param_gamma); 
     } 

     result_type min() const { return 0.0; } 
     result_type max() const { return 1.0; } 

     result_type a() const { return a_gamma.alpha(); } 
     result_type b() const { return b_gamma.alpha(); } 

     bool operator==(const beta_distribution<result_type>& other) const 
     { 
     return (param() == other.param() && 
       a_gamma == other.a_gamma && 
       b_gamma == other.b_gamma); 
     } 

     bool operator!=(const beta_distribution<result_type>& other) const 
     { 
     return !(*this == other); 
     } 

    private: 
     typedef std::gamma_distribution<result_type> gamma_dist_type; 

     gamma_dist_type a_gamma, b_gamma; 

     template <typename URNG> 
     result_type generate(URNG& engine, 
     gamma_dist_type& x_gamma, 
     gamma_dist_type& y_gamma) 
     { 
     result_type x = x_gamma(engine); 
     return x/(x + y_gamma(engine)); 
     } 
    }; 

    template <typename CharT, typename RealType> 
    std::basic_ostream<CharT>& operator<<(std::basic_ostream<CharT>& os, 
    const beta_distribution<RealType>& beta) 
    { 
    os << "~Beta(" << beta.a() << "," << beta.b() << ")"; 
    return os; 
    } 

    template <typename CharT, typename RealType> 
    std::basic_istream<CharT>& operator>>(std::basic_istream<CharT>& is, 
    beta_distribution<RealType>& beta) 
    { 
    std::string str; 
    RealType a, b; 
    if (std::getline(is, str, '(') && str == "~Beta" && 
     is >> a && is.get() == ',' && is >> b && is.get() == ')') { 
     beta = beta_distribution<RealType>(a, b); 
    } else { 
     is.setstate(std::ios::failbit); 
    } 
    return is; 
    } 

} 

使用它，像這樣：

std::random_device rd; 
std::mt19937 gen(rd()); 
sftrabbit::beta_distribution<> beta(2, 2); 
for (int i = 0; i < 10000; i++) { 
    std::cout << beta(gen) << std::endl; 
} 

Answer 2

查覈在NumPy隨機數發生器實現：NumPy distributions source

他們是用C語言實現，並且工作非常快。

Answer 3

也許您可以使用gsl用於生成隨機數字的代碼與beta分佈。他們使用一種奇怪的方式來生產它們，因爲你必須將一個隨機數發生器傳遞給函數，但肯定你可以得到你所需要的。

這裏的documentation和web page

Answer 4

升壓 「逆完全β」 是模擬貝塔另一個快速（簡單）的方式。

#include <random> 
#include <boost/math/special_functions/beta.hpp> 
template<typename URNG> 
double beta_sample(URNG& engine, double a, double b) 
{ 
    static std::uniform_real_distribution<double> unif(0,1); 
    double p = unif(engine); 
    return boost::math::ibeta_inv(a, b, p); 
    // Use Boost policies if it's not fast enough 
}