1
\$\begingroup\$

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.

\$\endgroup\$
3
\$\begingroup\$
  • 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)

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.