4
\$\begingroup\$

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 distriubtion, 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 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
\$\endgroup\$

Your Answer

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

Browse other questions tagged or ask your own question.