# Implementing a random distribution satisfying RandomNumberDistribution concept

Improving on this code review, I provide logarithmic_distribution to satisfy the C++ RandomNumberDistribution concept implementing the PDF

p(x) = N*(m + k * log(x))


with min <= x <= max. Parameters are m, k, min, max, while N is a normalization constant (to be determined from the parameters). Here is the header file

// file logarithmic_distribution.h
#ifndef logarithmic_distribution_h
#define logarithmic_distribution_h

#include <iostream>
#include <random>
#include <tuple>

namespace my_random {
/**
random distribution  p(x) = N*(m + k * log(x))   x in [min,max]

satisfying the RandomNumberDistribution concept, see
http://en.cppreference.com/w/cpp/named_req/RandomNumberDistribution
*/
struct logarithmic_distribution
{
using result_type = double;

struct param_type  // required by RandomNumberDistribution concept
{
using distribution_type = logarithmic_distribution;

/// constructor throws if pdf is not well defined, in particular if
/// min < 0, min >= max, or pdf(xmin) < 0.
param_type(double m, double k, double min, double max);

param_type(double m=0.0, double k=1.0, double max=2.0)
: param_type(m,k,1,max) {}

bool operator== (param_type const&p) const noexcept
{
return std::make_tuple(  m(),   k(),   min(),   max()) ==
std::make_tuple(p.m(), p.k(), p.min(), p.max());
}

bool operator!= (param_type const&p) const noexcept
{ return ! operator==(p); }

double m() const noexcept
{ return m_m; }

double k() const noexcept
{ return m_k; }

double min() const noexcept
{ return m_min; }

double max() const noexcept
{ return m_max; }

friend std::ostream&operator<<(std::ostream&, param_type const&);
friend std::istream&operator>>(std::istream&, param_type&);
protected:
double m_m;
double m_k;
double m_min, m_max;
};  // struct my_random::logarithmic_distribution::param_type

private:
struct auxiliary  // the actual implementation
: param_type
{
auxiliary(param_type const&) noexcept;

template<typename Generator>
result_type operator() (Generator &g) const
{ return sample(m_unif(g)); }

double pdf(double x) const noexcept;
double cdf(double x) const noexcept;
double mean() const noexcept;
double var() const noexcept;
private:
double m_mk;
double m_alna;
double m_cdfb, m_icdfb;
std::uniform_real_distribution<double> m_unif;
double aux_pdf(double) const noexcept;
double aux_cdf(double,double) const noexcept;
double aux_mean(double) const noexcept;
double aux_var(double) const noexcept;
result_type sample(double) const;
} m_aux;

public:
logarithmic_distribution(param_type const&p) noexcept
: m_aux(p) {}

logarithmic_distribution(double m, double k, double min, double max)
: logarithmic_distribution(param_type(m,k,min,max)) {}

logarithmic_distribution(double m=0.0, double k=1.0, double max=2.0)
: logarithmic_distribution(param_type(m,k,max)) {}

double m() const noexcept
{ return m_aux.m(); }

double k() const noexcept
{ return m_aux.k(); }

double min() const noexcept
{ return m_aux.min(); }

double max() const noexcept
{ return m_aux.max(); }

param_type const& param() const noexcept
{ return m_aux; }

bool operator==(logarithmic_distribution const&other) const noexcept
{ return param() == other.param(); }

bool operator!=(logarithmic_distribution const&other) const noexcept
{ return param() != other.param(); }

void param(param_type const&p)
{ m_aux = auxiliary(p); }

void reset() const noexcept {}

friend std::ostream&operator<<(std::ostream&os, logarithmic_distribution const&d)
{ return os << d.param(); }

friend std::istream&operator>>(std::istream&is, logarithmic_distribution&d)
{
param_type p;
is >> p;
d.param(p);
return is;
}

template<typename Generator>
result_type operator() (Generator &g) const
{ return m_aux(g); }

template<typename Generator>
result_type operator() (Generator &g, param_type const&p) const
{ return auxiliary(p)(g); }

double pdf(double x) const
{ return m_aux.pdf(x); }

double cdf(double x) const
{ return m_aux.cdf(x); }

double mean() const
{ return m_aux.mean(); }

double variance() const
{ return m_aux.var(); }
};  // struct my_random::logarithmic_distribution
}   // namespace my_random
#endif // logarithmic_distribution_h


And here the source code:

// file logarithmic_distribution.cc

#include "logarithmic_distribution.h"
#include <exception>
#include <string>
#include <limits>
#include <cmath>

namespace my_random {
using namespace std;

logarithmic_distribution::param_type::
param_type(double m, double k, double min, double max)
: m_m(m), m_k(k), m_min(min), m_max(max)
{
if(min <= 0.)
throw runtime_error("logarithmic pdf: min <= 0");
if(max <= min)
throw runtime_error("logarithmic pdf: min >= max");
if((m+k*log(min)) < 0.)
throw runtime_error("logarithmic pdf <0 at x=min");
}

logarithmic_distribution::auxiliary::auxiliary(param_type const&p) noexcept
: param_type(p)
, m_mk   (m_m - m_k)
, m_alna (m_min * log(m_min) )
, m_cdfb (m_mk*(m_max-m_min)+m_k*(m_max*log(m_max)-m_alna))
, m_icdfb(1/m_cdfb)
, m_unif (0.0,m_cdfb) {}

inline
double logarithmic_distribution::auxiliary::aux_pdf(double logx) const noexcept
{ return m_m + m_k * logx; }

double logarithmic_distribution::auxiliary::pdf(double x) const noexcept
{ return m_icdfb * aux_pdf(log(x)); }

inline
double logarithmic_distribution::auxiliary::aux_cdf(double x, double logx) const noexcept
{ return m_mk * (x-m_min) + m_k * (x*logx - m_alna); }

double logarithmic_distribution::auxiliary::cdf(double x) const noexcept
{ return m_icdfb * aux_cdf(x,log(x)); }

inline
double logarithmic_distribution::auxiliary::aux_mean(double x) const noexcept
{ return x*x*(m_m + m_k*(log(x)-0.5)); }

double logarithmic_distribution::auxiliary::mean() const noexcept
{ return 0.5*( aux_mean(m_max) - aux_mean(m_min) )*m_icdfb; }

inline
double logarithmic_distribution::auxiliary::aux_var(double x) const noexcept
{ return x*x*x*(m_m + m_k*(log(x)-0.33333333333333333)); }

double logarithmic_distribution::auxiliary::var() const noexcept
{
double mom1=mean();
double mom2=0.33333333333333333*(aux_var(m_max)-aux_var(m_min))*m_icdfb;
return mom2-mom1*mom1;
}

double logarithmic_distribution::auxiliary::sample(const double C) const
{
const unsigned max_iterations = 100;
const double eps = 10*numeric_limits<double>::epsilon();
if(C <= 0) return min();
if(C >= m_cdfb) return max();

// find x where cdf(x)=C using bracketed Newton-Raphson (NR)
double xl = min();
double xh = max();
double dxo= xh-xl;
double dx = dxo;
double x  = xl+C*dx/m_cdfb;                     // linear interpolation
double lx = log(x);
double df = aux_pdf(lx);
double f  = aux_cdf(x,lx) - C;
for(unsigned it=0; it!=max_iterations; ++it) {
double xo = x;
if((f>0? ((x-xh)*f>df) : ((x-xl)*f<df)) ||  // if NR goes outside of bracket
(abs(f+f) > abs(dxo*df))) {              // or NR convergence is slow
dxo= dx;                                //   bisection
dx = 0.5*(xh-xl);
x  = xl+dx;
} else {                                    // otherwise
dxo= dx;                                //   Newton-Raphson
dx = f/df;
x -= dx;
}
if((xo<=x && xo>=x) || abs(dx)<=eps*abs(x)) // test for convergence
return x;
lx = log(x);
df = aux_pdf(lx);
f  = aux_cdf(x,lx) - C;
if(f<0) xl=x; else xh=x;                    // maintain bracket
}
throw runtime_error("logarithmic_distribution: exceeded " +
to_string(max_iterations) + " iterations");
}

ostream&operator<<(ostream&os, logarithmic_distribution::param_type const&p)
{
using ios_base = 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(numeric_limits<double>::max_digits10);
os << p.m_m << space
<< p.m_k << space
<< p.m_min << space
<< p.m_max;
os.flags(flags);
os.fill(fill);
os.precision(prec);
return os;
}

istream&operator>>(istream&is, logarithmic_distribution::param_type &p)
{
using ios_base = istream::ios_base;
const auto flags = is.flags();
is.flags(ios_base::dec | ios_base::skipws);
is >> p.m_m >> p.m_k >> p.m_min >> p.m_max;
is.flags(flags);
return is;
}
} // namespace my_random


One particular question I have regards the design of logarithmic_distribution, which follows the following pattern

struct logarithmic_distribution
{
struct param_type { ... };
private:
struct auxiliary : param_type { /* full implementation */ }
m_aux;
public:
/* functionality implemented using m_aux */
};


An alternative would be

namespace details {
struct logarithmic_distribution_parameters { ... };
}
struct logarithmic_distribution
: details::logarithmic_distribution_parameters
{
using param_type = details::logarithmic_distribution_parameters;
/* full implementation similar to auxiliary above */
};


This avoids the re-direction (which however is fully inlined), but is somewhat arkward to implement (e.g.\ construction is best done using a setup() function). I wonder whether there are any principle design considerations that favour this alternative, or whether there is another yet better way.

• I don't advise importing all the names from std into your own namespace, even if only in the implementation file. That's a recipe for confusion and errors - in the worst case, a function added to later standards could be an unambiguously better match at a call site, as explained on Stack Overflow.

• I think that instead of std::runtime_error, it would be better to use the more specific std::range_error when out-of-range parameters are passed.

• Perhaps your constant 0.33333333333333333 would be better expressed as

constexpr auto one_third = 1.0 / 3.0;

• Does it make sense to use std::uniform_real_distribution<double> internally, or would it be better to use std::uniform_real_distribution<result_type>, in preparation for making it generic?

I had to fix a bug: because std::​uniform_real_distribution::​operator() is unfortunately non-const, we have to make my_random::​logarithmic_distribution::​operator() and my_random::​logarithmic_distribution::​auxiliary::​operator() non-const, too.

With that fixed, I made a test program, which generates plausible-looking results:

int main()
{
static constexpr std::size_t bins = 20;
std::array<std::size_t,bins> histogram{};

my_random::logarithmic_distribution dist{};

std::random_device rd{};
std::minstd_rand gen{rd()};

for (auto i = 0;  i < 1000000;  ++i) {
auto const val = dist(gen);
// store it into the right histogram bin
auto bin = (val - dist.min()) / (dist.max() - dist.min()) * histogram.size();
++histogram[bin];
}

auto i = 0;
for (auto bin: histogram) {
std::cout << "bin " << i++ << ": " << bin << '\n';
}
}


(partial review - I hope to add to this later)