1
\$\begingroup\$

Following from this post on SO, I provide logarithmic_distribution to satisfy the C++ RandomNumberDistribution concept implementing the PDF

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

with a <= x <= b. Parameters are m, k, a, b, while N is a normalisation constant. Here is the header file

// file logarithmic_distribution.h
#include <iostream>
#include <random>
namespace my_random {
  /// random distribution  p(x) = N*(m + k * log(x))
  ///
  /// satisfying the RandomNumberDistribution concept, see
  /// http://en.cppreference.com/w/cpp/concept/RandomNumberDistribution
  struct logarithmic_distribution
  {
    // public types
    using result_type = double;
    struct param_type  // required by RandomNumberDistribution concept
    {
      using distribution_type = logarithmic_distribution;
      param_type(double m, double k, double a, double b);
      param_type(double m, double k, double b)
        : param_type(m,k,1,b) {}
      param_type()
        : param_type(0,1,1,2) {}
      bool operator==(param_type const&p) const
      { return m()==p.m() && k()==p.k() && min()==p.min() && max()==p.max(); }
      double m() const { return _M_m; }
      double k() const { return _M_k; }
      double min() const { return _M_a; }
      double max() const { return _M_b; }
      friend std::ostream&operator<<(std::ostream&, param_type const&);
      friend std::istream&operator>>(std::istream&, param_type&);
    protected:
      double _M_m, _M_k, _M_a, _M_b;
    };
  private:
    param_type const& par() const { return _M_aux; }
    // actual implementation (could this all be hidden?)
    struct auxiliary
      : param_type
    {
      auxiliary(param_type const&);
      template<typename Generator>
      result_type operator()(Generator&g)
      { return sample(_M_unif(g)); }
      double pdf(double) const;
      double cdf(double) const;
      double mean() const;
      double var() const;
    private:
      double _M_mk, _M_alna, _M_cdfb, _M_icdfb;
      std::uniform_real_distribution<double> _M_unif;
      double _pdf(double) const;
      double _cdf(double,double) const;
      result_type sample(double) const;
    } _M_aux;
  public:
    // constructors
    logarithmic_distribution(param_type const&p)
      : _M_aux(p) {}
    logarithmic_distribution(double m, double k, double a, double b)
      : logarithmic_distribution(param_type(m,k,a,b)) {}
    logarithmic_distribution(double m, double k, double b)
      : logarithmic_distribution(param_type(m,k,b)) {}
    logarithmic_distribution()
      : logarithmic_distribution(param_type()) {}
    // getters
    double m() const { return _M_aux.m(); }
    double k() const { return _M_aux.k(); }
    double min() const { return _M_aux.min(); }
    double max() const { return _M_aux.max(); }
    param_type param() const { return _M_aux; }
    // comparison
    bool operator==(logarithmic_distribution const&other) const
    { return par() == other.par(); }
    bool operator!=(logarithmic_distribution const&other) const
    { return !(operator==(other)); }
    // setters
    void param(param_type const&p)
    { _M_aux = auxiliary(p); }
    void reset() const {}
    // I/O of state
    friend std::ostream&operator<<(std::ostream&os, logarithmic_distribution const&d)
    { return os << d.par(); }
    friend std::istream&operator>>(std::istream&is, logarithmic_distribution&d)
    {
      param_type p;
      is >> p;
      d.param(p);
      return is;
    }
    // random number generation
    template<typename Generator>
    result_type operator()(Generator&g)
    { return _M_aux(g); }
    template<typename Generator>
    result_type operator()(Generator&g, param_type const&p)
    { return auxiliary(p)(g); }
    // normalised pdf
    double pdf(double x) const { return _M_aux.pdf(x); }
    // normalised cdf
    double cdf(double x) const { return _M_aux.cdf(x); }
    // mean
    double mean() const { return _M_aux.mean(); }
    // variance
    double variance() const { return _M_aux.var(); }
  };
}

And here the source code:

// file logarithmic_distribution.cc
#include "logarithmic_distribution.h"
#include <exception>
#include <cmath>

namespace my_random {

  logarithmic_distribution::param_type::
  param_type(double m, double k, double a, double b)
    : _M_m(m), _M_k(k), _M_a(a), _M_b(b)
  {
    if(a <= 0.)
      throw std::runtime_error("logarithmic pdf: min <= 0");
    if(b <= a)
      throw std::runtime_error("logarithmic pdf: min >= max");
    if((m+k*std::log(a)) < 0.)
      throw std::runtime_error("logarithmic pdf <0 at x=min");
  }

  logarithmic_distribution::auxiliary::auxiliary(param_type const&p)
    : param_type(p)
    , _M_mk   (_M_m - _M_k)
    , _M_alna (_M_a * std::log(_M_a) )
    , _M_cdfb (_M_mk*(_M_b-_M_a)+_M_k*(_M_b*std::log(_M_b)-_M_alna))
    , _M_icdfb(1/_M_cdfb)
    , _M_unif (0.0,_M_cdfb) {}

  inline double logarithmic_distribution::auxiliary::_pdf(double logx) const
  { return _M_m + _M_k * logx; }

  inline double logarithmic_distribution::auxiliary::_cdf(double x,
                              double logx) const
  { return _M_mk * (x-_M_a) + _M_k * (x*logx - _M_alna); }

  double logarithmic_distribution::auxiliary::pdf(double x) const
  { return _M_icdfb*_pdf(std::log(x)); }

  double logarithmic_distribution::auxiliary::cdf(double x) const
  { return _M_icdfb*_cdf(x,std::log(x)); }

  double logarithmic_distribution::auxiliary::mean() const
  {
    auto func = [=](double x) { return x*x*(_M_m + _M_k*(std::log(x)-0.5)); };
    return 0.5*(func(_M_b)-func(_M_a))*_M_icdfb;
  }

  double logarithmic_distribution::auxiliary::var() const
  {
    const double third = 0.33333333333333333;
    auto func=[=](double x) { return x*x*x*(_M_m + _M_k*(std::log(x)-third)); };
    auto mom1=mean();
    auto mom2=third*(func(_M_b)-func(_M_a))*_M_icdfb;
    return mom2-mom1*mom1;
  }

  double logarithmic_distribution::auxiliary::sample(const double C) const
  {
    using std::abs;
    const unsigned max_iterations = 100;
    const double eps = 10*std::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
    auto xl = min();
    auto xh = max();
    auto dxo= xh-xl;
    auto dx = dxo;
    auto x  = xl+C*dx/_M_cdfb;                 // linear interpolation
    auto lx = std::log(x);
    auto df = _pdf(lx);
    auto f  = _cdf(x,lx) - C;
    for(unsigned it=0; it!=max_iterations; ++it) {
      auto xo = x;
      if((f>0? ((x-xh)*f>df) : ((x-xl)*f<df))  // NR goes outside of bracket
     || (abs(f+f) > abs(dxo*df))) {        // NR convergence is slow
    dxo= dx;                               // bisection
    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))                  // test for convergence
    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");
  }
  //
  std::ostream&operator<<(std::ostream&os,
              logarithmic_distribution::param_type 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._M_m << space
       << p._M_k << space
       << p._M_a << space
       << p._M_b;
    os.flags(flags);
    os.fill(fill);
    os.precision(prec);
    return os;
  }
  //
  std::istream&operator>>(std::istream&is,
              logarithmic_distribution::param_type &p)
  {
    using ios_base = std::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_a >> p._M_b;
    is.flags(flags);
    return is;
  }
}

I would welcome constructive critique. In particular is this the correct way of satisfying the RandomNumberDistribution concept? Is is possible to make the implementation more efficient and/or hide more details in the source file? Is a pimpl implementation preferrable?

Note, I'm not asking whether the code is mathematically correct. I have tested that myself already.

\$\endgroup\$
  • \$\begingroup\$ There is already a logarithmic distribution for random numbers in the standard: cplusplus.com/reference/random/lognormal_distribution \$\endgroup\$ – Martin York Apr 27 '17 at 14:53
  • \$\begingroup\$ I doubt anybody here has enough maths skills to check that the random numbers are logarithmically distributed. Can you provide a test harness so we can run it a 10 million times and plot the distribution. That way we should be able to see to an appropriate graph. \$\endgroup\$ – Martin York Apr 27 '17 at 14:54
  • \$\begingroup\$ @LokiAstari std::lognormal_distribution implements the log-normal distribution (i.e. a pdf where log(x) follows a normal distribution) while I'm implementing a distribution with the pdf given. \$\endgroup\$ – Walter Apr 27 '17 at 14:58
  • \$\begingroup\$ @LokiAstari I doubt your doubts. Also, I'm not asking for mathematical validation. \$\endgroup\$ – Walter Apr 27 '17 at 14:59
2
\$\begingroup\$

Overview

I am sure the maths is fine.
But overall very sloppy programming and as a result a maintenance nightmare for anybody to take over.

Validity

Looking at: http://en.cppreference.com/w/cpp/concept/RandomNumberDistribution seems like your functionality conforms.

Review

Your code is a bit compressed for my liking.

This is OK.

  double m() const { return _M_m; }
  double k() const { return _M_k; }
  double min() const { return _M_a; }
  double max() const { return _M_b; }

But when you try and do it on two lines it gets a bit crowded:

  bool operator==(param_type const&p) const
  { return m()==p.m() && k()==p.k() && min()==p.min() && max()==p.max(); }

Especially since you have very little (none) vertical white space. This makes the code very hard to read.

Please one variable per line.

  double _M_m, _M_k, _M_a, _M_b;

Also, each of those variables are reserved for the implementation (technically making this an illegal program). You should prefer not to use _ as the first character in an identifier the rules are complex and most people get it wrong. stick to m_<Name>.

Also this variables names are not very descriptive. Most of the time I would complain about it. But I have a feeling they are part of a very common maths expression and the names are standard (in which case it is fine but if it is not then you should be more expressive with the names).

You have an equality but not the inequality in param_type

  bool operator==(param_type const&p) const
  { return m()==p.m() && k()==p.k() && min()==p.min() && max()==p.max(); }

You can make this easier to write using tuples:

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

Note: std::tuple also implements operator< so it is useful to implement the strictly less than operator.

Did a teacher punish you for wasting space in another life?

result_type operator()(Generator&g)
                               ^^^^

Useless comments are worse than no comments.

// normalised pdf
double pdf(double x) const { return _M_aux.pdf(x); }
// normalised cdf
double cdf(double x) const { return _M_aux.cdf(x); }
// mean
double mean() const { return _M_aux.mean(); }
// variance
double variance() const { return _M_aux.var(); }

All the above are useless and give you no more information than reading the code. The trouble with useless comments is that they fall out of sync with the code and become a maintenance nightmare.

This does not make the code easier to read:

auto func=[=](double x) { return x*x*x*(_M_m + _M_k*(std::log(x)-third)); };

You are abusing lambda's here. Write a member function and pass parameters and call it like you would a normal function.

Tidy up your indetation this is impossible to follow:

    for(unsigned it=0; it!=max_iterations; ++it) {
      auto xo = x;
      if((f>0? ((x-xh)*f>df) : ((x-xl)*f<df))  // NR goes outside of bracket
     || (abs(f+f) > abs(dxo*df))) {        // NR convergence is slow
    dxo= dx;                               // bisection
    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))                  // test for convergence
    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");
  }
\$\endgroup\$
  • \$\begingroup\$ Thanks for (some or) your comments. The identation issues arose from pasting here. The _M_k was a slip and arose by taking guidance from stdlib implementations. lambda-use is not abuse: if that function is not used anywhere else, why define it? The condensed form is not my usual and arose to avoid too long a post. \$\endgroup\$ – Walter Apr 27 '17 at 17:03
  • \$\begingroup\$ Is it customary to improve/change the code after an answer has been given? \$\endgroup\$ – Walter Apr 27 '17 at 17:04
  • \$\begingroup\$ It is abuse of the lambda. The whole point in writting code in a higher level language is to make it readable by humans. That is the whole concept of self documenting code and having code reviews so that other people can read your code. The compiler and thus the machine will always be able to read it. You write in a style that does not make it easy for humans to read. \$\endgroup\$ – Martin York Apr 27 '17 at 19:05
  • \$\begingroup\$ Half the function I write are only used in one place. The point of writing them as a separate function is well discussed in many articles and can be summarized by the term "Self Documenting Code". \$\endgroup\$ – Martin York Apr 27 '17 at 19:07
  • \$\begingroup\$ As it stands I would never allow this to be checked into a larger code base. \$\endgroup\$ – Martin York Apr 27 '17 at 19:08

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.