I provide logarithmic_distribution to satisfy the C++ RandomNumberDistribution concept implementing the PDF
p(x) = Norm * (floor + logfac * log(x) )
with min <= x <= max. Parameters are floor, logfac, min, and max, while Norm is a normalisation constant.
This is the third, improved implementation, see also my previous post (the second version was removed in edit before any answer was received; this version is more similar to the first).
Apart from general fitness of the code, I have some specific questions (marked QUESTION in the code). Note that I'm not asking about the mathematical correctness of the implementation.
header file:
#include <iostream>
#include <random>
#include <tuple>
namespace my_random {
class logarithmic_distribution;
class parameter_for_logarithmic_distribution
: std::tuple<double,double,double,double>
{
using base = std::tuple<double,double,double,double>;
using self = parameter_for_logarithmic_distribution;
friend class logarithmic_distribution;
friend std::ostream&operator<<(std::ostream&, self const&);
friend std::istream&operator>>(std::istream&, self&);
void check_validity() const;
public:
using distribution_type = logarithmic_distribution;
// constructors: must support the same as for distribution_type
// (including auto-generated copy constructtor and assignment)
// QUESTION: I check validity of parameters here (and in stream input)
// Is this the correct/preferred approach? Or should we only
// check later, when the parameters are actually used (i.e. in
// the constructor of logarithmic_distribution::implementation)?
parameter_for_logarithmic_distribution(double floor, double logfac, double min, double max)
: base(floor,logfac,min,max)
{
check_validity();
}
parameter_for_logarithmic_distribution(double floor, double logfac, double max)
: parameter_for_logarithmic_distribution(floor,logfac,1.0,max)
{}
parameter_for_logarithmic_distribution()
: parameter_for_logarithmic_distribution(0,1,1,2)
{}
// comparison (strictly only == is required)
friend bool operator==(self const&p1, self const&p2)
{
return static_cast<base const&>(p1) == static_cast<base const&>(p2);
}
friend bool operator!=(self const&p1, self const&p2)
{
return !(p1==p2);
}
// getter methods
// will be inherited by distribution, thus automatically satisfying the
// corresponding requirement from concept RandomNumberDistribution
double floor() const { return std::get<0>(*this); }
double logfac() const { return std::get<1>(*this); }
double min() const { return std::get<2>(*this); }
double max() const { return std::get<3>(*this); }
}; // class my_random::parameter_for_logarithmic_distribution
// random distribution p(x) = Norm * (floor + logfac * log(x))
// satisfying the RandomNumberDistribution concept, see
// http://en.cppreference.com/w/cpp/concept/RandomNumberDistribution
class logarithmic_distribution
// QUESTION: I derive logarithmic_distribution from its param_type, which
// has the benefit that all the getters are inherited too, thus
// satisfying the corresponding requirement. The same holds for
// comparison and stream output.
// However, is the alternative (having a member 'parameters')
// perhaps preferrable?
: public parameter_for_logarithmic_distribution
{
// further data
double min_log_min;
double cdf_at_max;
std::uniform_real_distribution<double> uniform;
public:
// public types required for RandomNumberDistribution
using result_type = double;
using param_type = parameter_for_logarithmic_distribution;
// constructors and assignment as required for RandomNumberDistribution
// (including auto-generated copy constructor and assignment)
logarithmic_distribution(param_type const&);
logarithmic_distribution(double floor, double logfac, double min, double max)
: logarithmic_distribution(param_type(floor,logfac,min,max))
{}
logarithmic_distribution(double floor, double logfac, double max)
: logarithmic_distribution(param_type(floor,logfac,max))
{}
logarithmic_distribution()
: logarithmic_distribution(param_type())
{}
// setters as required for RandomNumberDistribution
void param(param_type const&p)
{
if(static_cast<param_type const&>(*this) != p)
operator=(logarithmic_distribution(p));
}
void reset() const {}
// getters as required for RandomNumberDistribution (including inherited)
param_type param() const { return *this; }
// state input (differs from parameter input)
friend std::istream&operator>>(std::istream&is, logarithmic_distribution&d)
{
param_type p;
is >> p;
d.param(p);
return is;
}
// random number generation interface (as required)
template<typename Generator>
result_type operator()(Generator &g)
{
return sample(uniform(g));
}
template<typename Generator>
result_type operator()(Generator &g, param_type const &p)
{
return logarithmic_distribution(p)(g);
}
// miscellaneous getters, not required
double pdf(double) const;
double cdf(double) const;
double mean() const;
double variance() const;
private:
double _pdf(double) const;
double _cdf(double, double) const;
double sample(double) const;
};// class my_random::logarithmic_distribution::implementation
} // namespace my_random
source file
#include "logarithmic_distribution.h"
#include <exception>
#include <cmath>
namespace {
// definite integral of PDF(x) from min to x
inline double log_dist_cdf(double x, double logx, double floor, double logfac, double min, double min_log_min)
{
return (floor-logfac) * (x-min) + logfac * (x*logx - min_log_min);
}
// indefinite integral of 2 * x * PDF(x), not normalised
inline double log_dist_for_mean(double x, double logx, double floor, double logfac)
{
return x*x*(floor+logfac*(logx-0.5));
}
// indefinite integral of 3 * x^2 * PDF(x), not normalised
inline double log_dist_for_var(double x, double logx, double floor, double logfac)
{
return x*x*x*(floor+logfac*(logx-0.3333333333333333333));
}
} // namespace {
namespace my_random {
void parameter_for_logarithmic_distribution::check_validity() const
{
if(min() <= 0.)
throw std::runtime_error("logarithmic pdf: min <= 0");
if(max() <= min())
throw std::runtime_error("logarithmic pdf: min >= max");
if((floor()+logfac()*std::log(min())) < 0.0)
throw std::runtime_error("logarithmic pdf <0 at x=min");
}
// I/O of parameters according to RandomNumberDistribution concept
//
// actually, the requirement is for the distribution_type, but we delegate
// to the param_type.
//
// note: must not change flags of I/O streams
std::ostream&operator<<(std::ostream&os, parameter_for_logarithmic_distribution const&p)
{
using ios_base = std::ostream::ios_base;
const auto flags = os.flags();
const auto fill = os.fill();
const auto prec = os.precision();
const auto space = os.widen(' ');
os.flags(ios_base::scientific | ios_base::left);
os.fill(space);
os.precision(std::numeric_limits<double>::max_digits10);
os << p.floor() << space
<< p.logfac() << space
<< p.min() << space
<< p.max();
os.flags(flags);
os.fill(fill);
os.precision(prec);
return os;
}
std::istream&operator>>(std::istream&is, parameter_for_logarithmic_distribution &p)
{
using ios_base = std::istream::ios_base;
const auto flags = is.flags();
is.flags(ios_base::dec | ios_base::skipws);
is >> std::get<0>(p)
>> std::get<1>(p)
>> std::get<2>(p)
>> std::get<3>(p);
is.flags(flags);
p.check_validity();
return is;
}
inline double logarithmic_distribution::_pdf(double logx) const
{
return floor() + logfac() * logx;
}
inline double logarithmic_distribution::_cdf(double x, double logx) const
{
return log_dist_cdf(x,logx,floor(),logfac(),min(),min_log_min);
}
logarithmic_distribution::logarithmic_distribution(param_type const&p)
: param_type (p)
, min_log_min (p.min() * std::log(p.min()))
, cdf_at_max (log_dist_cdf(p.max(),std::log(p.max()),
p.floor(),p.logfac(),p.min(),min_log_min))
, uniform (0.0,cdf_at_max)
{}
double logarithmic_distribution::pdf(double x) const
{
return _pdf(std::log(x))/cdf_at_max;
}
double logarithmic_distribution::cdf(double x) const
{
return _cdf(x,std::log(x))/cdf_at_max;
}
double logarithmic_distribution::mean() const
{
return 0.5*(log_dist_for_mean(max(),std::log(max()),floor(),logfac()) -
log_dist_for_mean(min(),std::log(min()),floor(),logfac())
)/cdf_at_max;
}
double logarithmic_distribution::variance() const
{
auto logmin = std::log(min());
auto logmax = std::log(max());
auto mom1 = 0.5*
(log_dist_for_mean(max(),logmax,floor(),logfac()) -
log_dist_for_mean(min(),logmin,floor(),logfac()))/cdf_at_max;
auto mom2 = 0.33333333333333333*
(log_dist_for_var (max(),logmax,floor(),logfac()) -
log_dist_for_var (min(),logmin,floor(),logfac()))/cdf_at_max;
return mom2-mom1*mom1;
}
// given a random value C for the non-normalised CDF(x), find x
double logarithmic_distribution::sample(const double C) const
{
using std::abs;
if(C <= 0)
return min();
if(C >= cdf_at_max)
return max();
// bracketed Newton-Raphson for f(x) = CDF(x)-C (not normalised)
// see also Numerical Recipies
const unsigned max_iterations = 100;
const double eps = 10*std::numeric_limits<double>::epsilon();
auto xl = min();
auto xh = max();
auto dxo= xh-xl;
auto dx = dxo;
auto x = xl+C*dx/cdf_at_max; // linear interpolation
auto lx = std::log(x);
auto df = _pdf(lx); // df/dx = PDF(x)
auto f = _cdf(x,lx) - C; // f(x) = CDF(x)-C
for(unsigned it=0; it!=max_iterations; ++it) {
auto xo = x;
if((f>0? ((x-xh)*f>df) : ((x-xl)*f<df)) // NR outside of bracket
|| (abs(f+f) > abs(dxo*df))) { // NR convergence is slow
dxo= dx; // bi-section
dx = 0.5*(xh-xl);
x = xl+dx;
} else {
dxo= dx; // Newton-Raphson
dx = f/df;
x -= dx;
}
if((xo<=x && xo>=x) || abs(dx)<=eps*abs(x)) // converged?
return x;
lx = std::log(x);
df = _pdf(lx);
f = _cdf(x,lx) - C;
if(f<0) xl=x; else xh=x; // maintain bracket
}
throw
std::runtime_error("logarithmic_distribution: exceeded 100 iterations");
}
} // namespace my_random
