具有Beta分布的随机数生成器

Question

具有Beta分布的随机数生成器

9

我需要一个类似于 betarand(a,b) 的 C 或 C++ 源代码的函数，该函数生成符合贝塔分布的随机数。我知道可以使用 boost 库，但我要将其移植到 CUDA 架构中，所以需要该源代码。有人能帮我吗？
同时，我有 betapdf（贝塔概率密度函数）。但我不知道如何使用它来创建随机数 :)

- s4eed

嗯，我想根据定义，如果数字以某种方式产生，可以被描述为分布的一部分，那么它们就不是随机的。:) 但无论如何..这个beta函数对你有用吗？此beta函数 - Mike

@Mike，你指的是哪个定义？你是在暗示线性分布比其他任何分布都“更随机”吗？ - Kos

根据http://en.wikipedia.org/wiki/Beta_distribution#Generating_beta-distributed_random_variates，您可以从2个伽马分布生成Beta分布。由于C++11提供了伽马分布，因此可以使用它们来创建Beta分布。 - Dave S

@sftrabbit~> 对我来说？！ :D 我迫不及待地等待着！ - s4eed

@saeed，提醒一下您——我已经更新了我的答案中的类，使其完全符合随机数分布要求。 - Joseph Mansfield

显示剩余4条评论

4个回答

2

也许您可以使用 gsl 用于生成 Beta 分布随机数的代码。它们使用了一种有点奇怪的方式来生成随机数，因为您需要将一个随机数生成器传递给函数，但是您肯定可以得到所需的内容。

这里是文档和网页。

- Noxbru

1

Boost的“反不完全Beta函数”是另一种快速（且简单）模拟Beta分布的方法。

#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
}

- Tanguy

0

查看在NumPy中的随机数生成器实现：NumPy分布源代码它们是用C语言实现的，并且运行速度非常快。

- meyumer

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Joseph Mansfield · Accepted Answer

C++11的随机数库没有提供beta分布。不过，可以用两个gamma分布来模拟beta分布，而这个库确实提供了gamma分布。我已经为您使用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;
}