dist_graphs.cpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

//  (C) Copyright John Maddock 2008.
//  Copyright Paul A. Bristow 2008
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifdef _MSC_VER
#  pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
#  pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
#  pragma warning (disable : 4512) // assignment operator could not be generated
#  pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
#  pragma warning (disable : 4127) // conditional expression is constant
#endif

#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error

#include <boost/math/distributions.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/svg_plot/svg_2d_plot.hpp>

#include <list>
#include <map>
#include <string>

template <class Dist>
struct is_discrete_distribution 
   : public boost::mpl::false_{};

template<class T, class P>
struct is_discrete_distribution<boost::math::bernoulli_distribution<T,P> > 
   : public boost::mpl::true_{};
template<class T, class P>
struct is_discrete_distribution<boost::math::binomial_distribution<T,P> > 
   : public boost::mpl::true_{};
template<class T, class P>
struct is_discrete_distribution<boost::math::negative_binomial_distribution<T,P> > 
   : public boost::mpl::true_{};
template<class T, class P>
struct is_discrete_distribution<boost::math::poisson_distribution<T,P> > 
   : public boost::mpl::true_{};


template <class Dist>
struct value_finder
{
   value_finder(Dist const& d, typename Dist::value_type v) 
      : m_dist(d), m_value(v) {}

   inline typename Dist::value_type operator()(const typename Dist::value_type& x)
   {
      return pdf(m_dist, x) - m_value;
   }

private:
   Dist m_dist;
   typename Dist::value_type m_value;
};

template <class Dist>
class distribution_plotter
{
public:
   distribution_plotter() : m_pdf(true), m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0) {}
   distribution_plotter(bool pdf) : m_pdf(pdf), m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0) {}

   void add(const Dist& d, const std::string& name)
   {
      //
      // Add to our list for later:
      //
      m_distributions.push_back(std::make_pair(name, d));
      //
      // Get the extent:
      //
      double a, b;
      std::tr1::tie(a, b) = support(d);
      //
      // PDF maximimum is at the mode:
      //
      double mod;
      try
      { 
         mod = mode(d); 
      } 
      catch(const std::domain_error& )
      {
         mod = a;
      }
      if((mod <= a) && !is_discrete_distribution<Dist>::value)
      {
         if((a != 0) && (fabs(a) > 1e-2))
            mod = a * (1 + 1e-2);
         else
            mod = 1e-2;
      }
      double peek_y = pdf(d, mod);
      double min_y = peek_y / 20;
      //
      // If the extent is "infinite" then find out how large it 
      // has to be for the PDF to decay to min_y:
      //
      if(a <= -(std::numeric_limits<double>::max)())
      {
         boost::uintmax_t max_iter = 500;
         double guess = mod;
         if((pdf(d, 0) > min_y) || (guess == 0))
            guess = -1e-3;
         a = boost::math::tools::bracket_and_solve_root(
            value_finder<Dist>(d, min_y),
            guess,
            8.0,
            true,
            boost::math::tools::eps_tolerance<double>(10),
            max_iter).first;
      }
      if(b >= (std::numeric_limits<double>::max)())
      {
         boost::uintmax_t max_iter = 500;
         double guess = mod;
         if(a <= 0)
            if((pdf(d, 0) > min_y) || (guess == 0))
               guess = 1e-3;
         b = boost::math::tools::bracket_and_solve_root(
            value_finder<Dist>(d, min_y),
            guess,
            8.0,
            false,
            boost::math::tools::eps_tolerance<double>(10),
            max_iter).first;
      }
      //
      // Recalculate peek_y and location of mod so that
      // it's not too close to one end of the graph: 
      // otherwise we may be shooting off to infinity.
      //
      if(!is_discrete_distribution<Dist>::value)
      {
         if(mod <= a + (b-a)/50)
         {
            mod = a + (b-a)/50;
         }
         if(mod >= b - (b-a)/50)
         {
            mod = b - (b-a)/50;
         }
         peek_y = pdf(d, mod);
      }
      //
      // Now set our limits:
      //
      if(peek_y > m_max_y)
         m_max_y = peek_y;
      if(m_max_x == m_min_x)
      {
         m_max_x = b;
         m_min_x = a;
      }
      else
      {
         if(a < m_min_x) 
            m_min_x = a;
         if(b > m_max_x)
            m_max_x = b;
      }
   }

   void plot(const std::string& title, const std::string& file)
   {
      using namespace boost::svg;

      static const svg_color colors[5] = 
      {
         darkblue,
         darkred,
         darkgreen,
         darkorange,
         chartreuse
      };

      if(m_pdf == false)
      {
         m_min_y = 0;
         m_max_y = 1;
      }

      svg_2d_plot plot;
      plot.image_size(750, 400);
      plot.coord_precision(4); // Avoids any visible steps.
      plot.title_font_size(20);
      plot.legend_title_font_size(15);
      plot.title(title);
      if((m_distributions.size() == 1) && (m_distributions.begin()->first == ""))
         plot.legend_on(false);
      else
         plot.legend_on(true);
      plot.title_on(true);
      //plot.x_major_labels_on(true).y_major_labels_on(true);
      //double x_delta = (m_max_x - m_min_x) / 10;
      double y_delta = (m_max_y - m_min_y) / 10;
      if(is_discrete_distribution<Dist>::value)
         plot.x_range(m_min_x - 0.5, m_max_x + 0.5)
             .y_range(m_min_y, m_max_y + y_delta);
      else
         plot.x_range(m_min_x, m_max_x)
             .y_range(m_min_y, m_max_y + y_delta);
      plot.x_label_on(true).x_label("Random Variable");
      plot.y_label_on(true).y_label("Probability");
      plot.plot_border_color(lightslategray)
          .background_border_color(lightslategray)
          .legend_border_color(lightslategray)
          .legend_background_color(white);
      //
      // Work out axis tick intervals:
      //
      double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
      double interval = std::pow(10.0, (int)l);
      if(((m_max_x - m_min_x) / interval) > 10)
         interval *= 5;
      if(is_discrete_distribution<Dist>::value)
      {
         interval = interval > 1 ? std::floor(interval) : 1;
         plot.x_num_minor_ticks(0);
      }
      plot.x_major_interval(interval);
      l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
      interval = std::pow(10.0, (int)l);
      if(((m_max_y - m_min_y) / interval) > 10)
         interval *= 5;
      plot.y_major_interval(interval);

      int color_index = 0;

      if(!is_discrete_distribution<Dist>::value)
      {
         //
         // Continuous distribution:
         //
         for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
            i != m_distributions.end(); ++i)
         {
            double x = m_min_x;
            double interval = (m_max_x - m_min_x) / 200;
            std::map<double, double> data;
            while(x <= m_max_x)
            {
               data[x] = m_pdf ? pdf(i->second, x) : cdf(i->second, x);
               x += interval;
            }
            plot.plot(data, i->first)
               .line_on(true)
               .line_color(colors[color_index])
               .line_width(1.)
               .shape(none);
               //.bezier_on(true) // Bezier can't cope with badly behaved like uniform & triangular.
            ++color_index;
            color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
         }
      }
      else
      {
         //
         // Discrete distribution:
         //
         double x_width = 0.75 / m_distributions.size();
         double x_off = -0.5 * 0.75;
         for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
            i != m_distributions.end(); ++i)
         {
            double x = ceil(m_min_x);
            double interval = 1;
            std::map<double, double> data;
            while(x <= m_max_x)
            {
               double p;
               try{
                  p = m_pdf ? pdf(i->second, x) : cdf(i->second, x);
               }
               catch(const std::domain_error&)
               {
                  p = 0;
               }
               data[x + x_off] = 0;
               data[x + x_off + 0.00001] = p;
               data[x + x_off + x_width] = p;
               data[x + x_off + x_width + 0.00001] = 0;
               x += interval;
            }
            x_off += x_width;
            svg_2d_plot_series& s = plot.plot(data, i->first);
            s.line_on(true)
               .line_color(colors[color_index])
               .line_width(1.)