Skip to content

Instantly share code, notes, and snippets.

@sftrabbit
Last active July 24, 2024 10:22
Show Gist options
  • Save sftrabbit/5068941 to your computer and use it in GitHub Desktop.
Save sftrabbit/5068941 to your computer and use it in GitHub Desktop.
A beta distribution random number distribution for C++11. Conforms to the requirements for a random number distribution (§26.5.1.6 [rand.req.dist]) and is therefore compatible with the <random> C++ library header. Implemented in terms of two gamma distributions (std::gamma_distribution): if X and Y are independently distributed Γ(α, θ) and Γ(β, …
#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; }
RealType a() const { return a_gamma.alpha(); }
RealType 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;
}
}
@horizon-z40
Copy link

Can this be implemented on c++98 ?

@nghiadt05
Copy link

Thank you very much for the code.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment