ceph/src/boost/libs/math/doc/graphs/dist_graphs.cpp

   1 /*! \file dist_graphs.cpp
   2     \brief Produces Scalable Vector Graphic (.svg) files for all distributions.
   3     \details These files can be viewed using most browsers,
   4     though MS Internet Explorer requires a plugin from Adobe.
   5     These file can be converted to .png using Inkscape
   6     (see www.inkscape.org) Export Bit option which by default produces
   7     a Portable Network Graphic file with that same filename but .png suffix instead of .svg.
   8     Using Python, generate.sh does this conversion automatically for all .svg files in a folder.
   9
  10     \author John Maddock and Paul A. Bristow
  11   */
  12 //  Copyright John Maddock 2008.
  13 //  Copyright Paul A. Bristow 2008, 2009, 2012
  14 //  Use, modification and distribution are subject to the
  15 //  Boost Software License, Version 1.0. (See accompanying file
  16 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  17
  18 #ifdef _MSC_VER
  19 #  pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
  20 #  pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
  21 #  pragma warning (disable : 4512) // assignment operator could not be generated
  22 #  pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
  23 #  pragma warning (disable : 4127) // conditional expression is constant
  24 #endif
  25
  26 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  27
  28 #include <boost/math/distributions.hpp>
  29 #include <boost/math/tools/roots.hpp>
  30 #include <boost/svg_plot/svg_2d_plot.hpp>
  31
  32 #include <list>
  33 #include <map>
  34 #include <string>
  35
  36 template <class Dist>
  37 struct is_discrete_distribution
  38    : public boost::mpl::false_{};
  39
  40 template<class T, class P>
  41 struct is_discrete_distribution<boost::math::bernoulli_distribution<T,P> >
  42    : public boost::mpl::true_{};
  43 template<class T, class P>
  44 struct is_discrete_distribution<boost::math::binomial_distribution<T,P> >
  45    : public boost::mpl::true_{};
  46 template<class T, class P>
  47 struct is_discrete_distribution<boost::math::negative_binomial_distribution<T,P> >
  48    : public boost::mpl::true_{};
  49 template<class T, class P>
  50 struct is_discrete_distribution<boost::math::poisson_distribution<T,P> >
  51    : public boost::mpl::true_{};
  52 template<class T, class P>
  53 struct is_discrete_distribution<boost::math::hypergeometric_distribution<T,P> >
  54    : public boost::mpl::true_{};
  55
  56
  57 template <class Dist>
  58 struct value_finder
  59 {
  60    value_finder(Dist const& d, typename Dist::value_type v)
  61       : m_dist(d), m_value(v) {}
  62
  63    inline typename Dist::value_type operator()(const typename Dist::value_type& x)
  64    {
  65       return pdf(m_dist, x) - m_value;
  66    }
  67
  68 private:
  69    Dist m_dist;
  70    typename Dist::value_type m_value;
  71 };
  72
  73 template <class Dist>
  74 class distribution_plotter
  75 {
  76 public:
  77    distribution_plotter() : m_pdf(true), m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0) {}
  78    distribution_plotter(bool pdf) : m_pdf(pdf), m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0) {}
  79
  80    void add(const Dist& d, const std::string& name)
  81    {
  82       // Add name of distribution to our list for later:
  83       m_distributions.push_back(std::make_pair(name, d));
  84       //
  85       // Get the extent of the distribution from the support:
  86       double a, b;
  87       std::tr1::tie(a, b) = support(d);
  88       //
  89       // PDF maximimum is at the mode (probably):
  90       double mod;
  91       try
  92       {
  93          mod = mode(d);
  94       }
  95       catch(const std::domain_error& )
  96       { // but if not use the lower limit of support.
  97          mod = a;
  98       }
  99       if((mod <= a) && !is_discrete_distribution<Dist>::value)
 100       { // Continuous distribution at or below lower limit of support.
 101         double margin = 1e-2; // Margin of 1% (say) to get lowest off the 'end stop'.
 102          if((a != 0) && (fabs(a) > margin))
 103          {
 104             mod = a * (1 + ((a > 0) ? margin : -margin));
 105          }
 106          else
 107          { // Case of mod near zero?
 108             mod = margin;
 109          }
 110       }
 111       double peek_y = pdf(d, mod);
 112       double min_y = peek_y / 20;
 113       //
 114       // If the extent is "infinite" then find out how large it
 115       // has to be for the PDF to decay to min_y:
 116       //
 117       if(a <= -(std::numeric_limits<double>::max)())
 118       {
 119          boost::uintmax_t max_iter = 500;
 120          double guess = mod;
 121          if((pdf(d, 0) > min_y) || (guess == 0))
 122             guess = -1e-3;
 123          a = boost::math::tools::bracket_and_solve_root(
 124             value_finder<Dist>(d, min_y),
 125             guess,
 126             8.0,
 127             true,
 128             boost::math::tools::eps_tolerance<double>(10),
 129             max_iter).first;
 130       }
 131       if(b >= (std::numeric_limits<double>::max)())
 132       {
 133          boost::uintmax_t max_iter = 500;
 134          double guess = mod;
 135          if(a <= 0)
 136             if((pdf(d, 0) > min_y) || (guess == 0))
 137                guess = 1e-3;
 138          b = boost::math::tools::bracket_and_solve_root(
 139             value_finder<Dist>(d, min_y),
 140             guess,
 141             8.0,
 142             false,
 143             boost::math::tools::eps_tolerance<double>(10),
 144             max_iter).first;
 145       }
 146       //
 147       // Recalculate peek_y and location of mod so that
 148       // it's not too close to one end of the graph:
 149       // otherwise we may be shooting off to infinity.
 150       //
 151       if(!is_discrete_distribution<Dist>::value)
 152       {
 153          if(mod <= a + (b-a)/50)
 154          {
 155             mod = a + (b-a)/50;
 156          }
 157          if(mod >= b - (b-a)/50)
 158          {
 159             mod = b - (b-a)/50;
 160          }
 161          peek_y = pdf(d, mod);
 162       }
 163       //
 164       // Now set our limits:
 165       //
 166       if(peek_y > m_max_y)
 167          m_max_y = peek_y;
 168       if(m_max_x == m_min_x)
 169       {
 170          m_max_x = b;
 171          m_min_x = a;
 172       }
 173       else
 174       {
 175          if(a < m_min_x)
 176             m_min_x = a;
 177          if(b > m_max_x)
 178             m_max_x = b;
 179       }
 180    }
 181
 182    void plot(const std::string& title, const std::string& file)
 183    {
 184       using namespace boost::svg;
 185
 186       static const svg_color colors[5] =
 187       {
 188          darkblue,
 189          darkred,
 190          darkgreen,
 191          darkorange,
 192          chartreuse
 193       };
 194
 195       if(m_pdf == false)
 196       {
 197          m_min_y = 0;
 198          m_max_y = 1;
 199       }
 200
 201       svg_2d_plot plot;
 202       plot.image_x_size(750);
 203       plot.image_y_size(400);
 204       plot.coord_precision(4); // Avoids any visible steps.
 205       plot.title_font_size(20);
 206       plot.legend_title_font_size(15);
 207       plot.title(title);
 208       if((m_distributions.size() == 1) && (m_distributions.begin()->first == ""))
 209          plot.legend_on(false);
 210       else
 211          plot.legend_on(true);
 212       plot.title_on(true);
 213       //plot.x_major_labels_on(true).y_major_labels_on(true);
 214       //double x_delta = (m_max_x - m_min_x) / 10;
 215       double y_delta = (m_max_y - m_min_y) / 10;
 216       if(is_discrete_distribution<Dist>::value)
 217          plot.x_range(m_min_x - 0.5, m_max_x + 0.5)
 218              .y_range(m_min_y, m_max_y + y_delta);
 219       else
 220          plot.x_range(m_min_x, m_max_x)
 221              .y_range(m_min_y, m_max_y + y_delta);
 222       plot.x_label_on(true).x_label("Random Variable");
 223       plot.y_label_on(true).y_label("Probability");
 224       plot.plot_border_color(lightslategray)
 225           .background_border_color(lightslategray)
 226           .legend_border_color(lightslategray)
 227           .legend_background_color(white);
 228       //
 229       // Work out axis tick intervals:
 230       //
 231       double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
 232       double interval = std::pow(10.0, (int)l);
 233       if(((m_max_x - m_min_x) / interval) > 10)
 234          interval *= 5;
 235       if(is_discrete_distribution<Dist>::value)
 236       {
 237          interval = interval > 1 ? std::floor(interval) : 1;
 238          plot.x_num_minor_ticks(0);
 239       }
 240       plot.x_major_interval(interval);
 241       l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
 242       interval = std::pow(10.0, (int)l);
 243       if(((m_max_y - m_min_y) / interval) > 10)
 244          interval *= 5;
 245       plot.y_major_interval(interval);
 246
 247       int color_index = 0;
 248
 249       if(!is_discrete_distribution<Dist>::value)
 250       {
 251          //
 252          // Continuous distribution:
 253          //
 254          for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
 255             i != m_distributions.end(); ++i)
 256          {
 257             double x = m_min_x;
 258             double interval = (m_max_x - m_min_x) / 200;
 259             std::map<double, double> data;
 260             while(x <= m_max_x)
 261             {
 262                data[x] = m_pdf ? pdf(i->second, x) : cdf(i->second, x);
 263                x += interval;
 264             }
 265             plot.plot(data, i->first)
 266                .line_on(true)
 267                .line_color(colors[color_index])
 268                .line_width(1.)
 269                .shape(none);
 270
 271                //.bezier_on(true) // Bezier can't cope with badly behaved like uniform & triangular.
 272             ++color_index;
 273             color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
 274          }
 275       }
 276       else
 277       {
 278          //
 279          // Discrete distribution:
 280          //
 281          double x_width = 0.75 / m_distributions.size();
 282          double x_off = -0.5 * 0.75;
 283          for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
 284             i != m_distributions.end(); ++i)
 285          {
 286             double x = ceil(m_min_x);
 287             double interval = 1;
 288             std::map<double, double> data;
 289             while(x <= m_max_x)
 290             {
 291                double p;
 292                try{
 293                   p = m_pdf ? pdf(i->second, x) : cdf(i->second, x);
 294                }
 295                catch(const std::domain_error&)
 296                {
 297                   p = 0;
 298                }
 299                data[x + x_off] = 0;
 300                data[x + x_off + 0.00001] = p;
 301                data[x + x_off + x_width] = p;
 302                data[x + x_off + x_width + 0.00001] = 0;
 303                x += interval;
 304             }
 305             x_off += x_width;
 306             svg_2d_plot_series& s = plot.plot(data, i->first);
 307             s.line_on(true)
 308                .line_color(colors[color_index])
 309                .line_width(1.)
 310                .shape(none)
 311                .area_fill(colors[color_index]);
 312             ++color_index;
 313             color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
 314          }
 315       }
 316       plot.write(file);
 317    }
 318
 319 private:
 320    bool m_pdf;
 321    std::list<std::pair<std::string, Dist> > m_distributions;
 322    double m_min_x, m_max_x, m_min_y, m_max_y;
 323 };
 324
 325 int main()
 326 {
 327   try
 328   {
 329
 330    distribution_plotter<boost::math::gamma_distribution<> >
 331       gamma_plotter;
 332    gamma_plotter.add(boost::math::gamma_distribution<>(0.75), "shape = 0.75");
 333    gamma_plotter.add(boost::math::gamma_distribution<>(1), "shape = 1");
 334    gamma_plotter.add(boost::math::gamma_distribution<>(3), "shape = 3");
 335    gamma_plotter.plot("Gamma Distribution PDF With Scale = 1", "gamma1_pdf.svg");
 336
 337    distribution_plotter<boost::math::gamma_distribution<> >
 338       gamma_plotter2;
 339    gamma_plotter2.add(boost::math::gamma_distribution<>(2, 0.5), "scale = 0.5");
 340    gamma_plotter2.add(boost::math::gamma_distribution<>(2, 1), "scale = 1");
 341    gamma_plotter2.add(boost::math::gamma_distribution<>(2, 2), "scale = 2");
 342    gamma_plotter2.plot("Gamma Distribution PDF With Shape = 2", "gamma2_pdf.svg");
 343
 344    distribution_plotter<boost::math::normal>
 345       normal_plotter;
 346    normal_plotter.add(boost::math::normal(0, 1), "&#x3BC; = 0, &#x3C3; = 1");
 347    normal_plotter.add(boost::math::normal(0, 0.5), "&#x3BC; = 0, &#x3C3; = 0.5");
 348    normal_plotter.add(boost::math::normal(0, 2), "&#x3BC; = 0, &#x3C3; = 2");
 349    normal_plotter.add(boost::math::normal(-1, 1), "&#x3BC; = -1, &#x3C3; = 1");
 350    normal_plotter.add(boost::math::normal(1, 1), "&#x3BC; = 1, &#x3C3; = 1");
 351    normal_plotter.plot("Normal Distribution PDF", "normal_pdf.svg");
 352
 353    distribution_plotter<boost::math::laplace>
 354       laplace_plotter;
 355    laplace_plotter.add(boost::math::laplace(0, 1), "&#x3BC; = 0, &#x3C3; = 1");
 356    laplace_plotter.add(boost::math::laplace(0, 0.5), "&#x3BC; = 0, &#x3C3; = 0.5");
 357    laplace_plotter.add(boost::math::laplace(0, 2), "&#x3BC; = 0, &#x3C3; = 2");
 358    laplace_plotter.add(boost::math::laplace(-1, 1), "&#x3BC; = -1, &#x3C3; = 1");
 359    laplace_plotter.add(boost::math::laplace(1, 1), "&#x3BC; = 1, &#x3C3; = 1");
 360    laplace_plotter.plot("Laplace Distribution PDF", "laplace_pdf.svg");
 361
 362    distribution_plotter<boost::math::non_central_chi_squared>
 363       nc_cs_plotter;
 364    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 0), "v=20, &#x3BB;=0");
 365    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 1), "v=20, &#x3BB;=1");
 366    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 5), "v=20, &#x3BB;=5");
 367    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 10), "v=20, &#x3BB;=10");
 368    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 20), "v=20, &#x3BB;=20");
 369    nc_cs_plotter.add(boost::math::non_central_chi_squared(20, 100), "v=20, &#x3BB;=100");
 370    nc_cs_plotter.plot("Non Central Chi Squared PDF", "nccs_pdf.svg");
 371
 372    distribution_plotter<boost::math::non_central_beta>
 373       nc_beta_plotter;
 374    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 0), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=0");
 375    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 1), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=1");
 376    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 5), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=5");
 377    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 10), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=10");
 378    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 40), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=40");
 379    nc_beta_plotter.add(boost::math::non_central_beta(10, 15, 100), "&#x3B1;=10, &#x3B2;=15, &#x3B4;=100");
 380    nc_beta_plotter.plot("Non Central Beta PDF", "nc_beta_pdf.svg");
 381
 382    distribution_plotter<boost::math::non_central_f>
 383       nc_f_plotter;
 384    nc_f_plotter.add(boost::math::non_central_f(10, 20, 0), "v1=10, v2=20, &#x3BB;=0");
 385    nc_f_plotter.add(boost::math::non_central_f(10, 20, 1), "v1=10, v2=20, &#x3BB;=1");
 386    nc_f_plotter.add(boost::math::non_central_f(10, 20, 5), "v1=10, v2=20, &#x3BB;=5");
 387    nc_f_plotter.add(boost::math::non_central_f(10, 20, 10), "v1=10, v2=20, &#x3BB;=10");
 388    nc_f_plotter.add(boost::math::non_central_f(10, 20, 40), "v1=10, v2=20, &#x3BB;=40");
 389    nc_f_plotter.add(boost::math::non_central_f(10, 20, 100), "v1=10, v2=20, &#x3BB;=100");
 390    nc_f_plotter.plot("Non Central F PDF", "nc_f_pdf.svg");
 391
 392    distribution_plotter<boost::math::non_central_t>
 393       nc_t_plotter;
 394    nc_t_plotter.add(boost::math::non_central_t(10, -10), "v=10, &#x3B4;=-10");
 395    nc_t_plotter.add(boost::math::non_central_t(10, -5), "v=10, &#x3B4;=-5");
 396    nc_t_plotter.add(boost::math::non_central_t(10, 0), "v=10, &#x3B4;=0");
 397    nc_t_plotter.add(boost::math::non_central_t(10, 5), "v=10, &#x3B4;=5");
 398    nc_t_plotter.add(boost::math::non_central_t(10, 10), "v=10, &#x3B4;=10");
 399    nc_t_plotter.add(boost::math::non_central_t(std::numeric_limits<double>::infinity(), 15), "v=inf, &#x3B4;=15");
 400    nc_t_plotter.plot("Non Central T PDF", "nc_t_pdf.svg");
 401
 402    distribution_plotter<boost::math::non_central_t>
 403      nc_t_CDF_plotter(false);
 404    nc_t_CDF_plotter.add(boost::math::non_central_t(10, -10), "v=10, &#x3B4;=-10");
 405    nc_t_CDF_plotter.add(boost::math::non_central_t(10, -5), "v=10, &#x3B4;=-5");
 406    nc_t_CDF_plotter.add(boost::math::non_central_t(10, 0), "v=10, &#x3B4;=0");
 407    nc_t_CDF_plotter.add(boost::math::non_central_t(10, 5), "v=10, &#x3B4;=5");
 408    nc_t_CDF_plotter.add(boost::math::non_central_t(10, 10), "v=10, &#x3B4;=10");
 409    nc_t_CDF_plotter.add(boost::math::non_central_t(std::numeric_limits<double>::infinity(), 15), "v=inf, &#x3B4;=15");
 410    nc_t_CDF_plotter.plot("Non Central T CDF", "nc_t_cdf.svg");
 411
 412    distribution_plotter<boost::math::beta_distribution<> >
 413       beta_plotter;
 414    beta_plotter.add(boost::math::beta_distribution<>(0.5, 0.5), "alpha=0.5, beta=0.5");
 415    beta_plotter.add(boost::math::beta_distribution<>(5, 1), "alpha=5, beta=1");
 416    beta_plotter.add(boost::math::beta_distribution<>(1, 3), "alpha=1, beta=3");
 417    beta_plotter.add(boost::math::beta_distribution<>(2, 2), "alpha=2, beta=2");
 418    beta_plotter.add(boost::math::beta_distribution<>(2, 5), "alpha=2, beta=5");
 419    beta_plotter.plot("Beta Distribution PDF", "beta_pdf.svg");
 420
 421    distribution_plotter<boost::math::cauchy_distribution<> >
 422       cauchy_plotter;
 423    cauchy_plotter.add(boost::math::cauchy_distribution<>(-5, 1), "location = -5");
 424    cauchy_plotter.add(boost::math::cauchy_distribution<>(0, 1), "location = 0");
 425    cauchy_plotter.add(boost::math::cauchy_distribution<>(5, 1), "location = 5");
 426    cauchy_plotter.plot("Cauchy Distribution PDF (scale = 1)", "cauchy_pdf1.svg");
 427
 428    distribution_plotter<boost::math::cauchy_distribution<> >
 429       cauchy_plotter2;
 430    cauchy_plotter2.add(boost::math::cauchy_distribution<>(0, 0.5), "scale = 0.5");
 431    cauchy_plotter2.add(boost::math::cauchy_distribution<>(0, 1), "scale = 1");
 432    cauchy_plotter2.add(boost::math::cauchy_distribution<>(0, 2), "scale = 2");
 433    cauchy_plotter2.plot("Cauchy Distribution PDF (location = 0)", "cauchy_pdf2.svg");
 434
 435    distribution_plotter<boost::math::chi_squared_distribution<> >
 436       chi_squared_plotter;
 437    //chi_squared_plotter.add(boost::math::chi_squared_distribution<>(1), "v=1");
 438    chi_squared_plotter.add(boost::math::chi_squared_distribution<>(2), "v=2");
 439    chi_squared_plotter.add(boost::math::chi_squared_distribution<>(5), "v=5");
 440    chi_squared_plotter.add(boost::math::chi_squared_distribution<>(10), "v=10");
 441    chi_squared_plotter.plot("Chi Squared Distribution PDF", "chi_squared_pdf.svg");
 442
 443    distribution_plotter<boost::math::exponential_distribution<> >
 444       exponential_plotter;
 445    exponential_plotter.add(boost::math::exponential_distribution<>(0.5), "&#x3BB;=0.5");
 446    exponential_plotter.add(boost::math::exponential_distribution<>(1), "&#x3BB;=1");
 447    exponential_plotter.add(boost::math::exponential_distribution<>(2), "&#x3BB;=2");
 448    exponential_plotter.plot("Exponential Distribution PDF", "exponential_pdf.svg");
 449
 450    distribution_plotter<boost::math::extreme_value_distribution<> >
 451       extreme_value_plotter;
 452    extreme_value_plotter.add(boost::math::extreme_value_distribution<>(-5), "location=-5");
 453    extreme_value_plotter.add(boost::math::extreme_value_distribution<>(0), "location=0");
 454    extreme_value_plotter.add(boost::math::extreme_value_distribution<>(5), "location=5");
 455    extreme_value_plotter.plot("Extreme Value Distribution PDF (shape=1)", "extreme_value_pdf1.svg");
 456
 457    distribution_plotter<boost::math::extreme_value_distribution<> >
 458       extreme_value_plotter2;
 459    extreme_value_plotter2.add(boost::math::extreme_value_distribution<>(0, 0.5), "shape=0.5");
 460    extreme_value_plotter2.add(boost::math::extreme_value_distribution<>(0, 1), "shape=1");
 461    extreme_value_plotter2.add(boost::math::extreme_value_distribution<>(0, 2), "shape=2");
 462    extreme_value_plotter2.plot("Extreme Value Distribution PDF (location=0)", "extreme_value_pdf2.svg");
 463
 464    distribution_plotter<boost::math::fisher_f_distribution<> >
 465       fisher_f_plotter;
 466    fisher_f_plotter.add(boost::math::fisher_f_distribution<>(4, 4), "n=4, m=4");
 467    fisher_f_plotter.add(boost::math::fisher_f_distribution<>(10, 4), "n=10, m=4");
 468    fisher_f_plotter.add(boost::math::fisher_f_distribution<>(10, 10), "n=10, m=10");
 469    fisher_f_plotter.add(boost::math::fisher_f_distribution<>(4, 10), "n=4, m=10");
 470    fisher_f_plotter.plot("F Distribution PDF", "fisher_f_pdf.svg");
 471
 472    distribution_plotter<boost::math::lognormal_distribution<> >
 473       lognormal_plotter;
 474    lognormal_plotter.add(boost::math::lognormal_distribution<>(-1), "location=-1");
 475    lognormal_plotter.add(boost::math::lognormal_distribution<>(0), "location=0");
 476    lognormal_plotter.add(boost::math::lognormal_distribution<>(1), "location=1");
 477    lognormal_plotter.plot("Lognormal Distribution PDF (scale=1)", "lognormal_pdf1.svg");
 478
 479    distribution_plotter<boost::math::lognormal_distribution<> >
 480       lognormal_plotter2;
 481    lognormal_plotter2.add(boost::math::lognormal_distribution<>(0, 0.5), "scale=0.5");
 482    lognormal_plotter2.add(boost::math::lognormal_distribution<>(0, 1), "scale=1");
 483    lognormal_plotter2.add(boost::math::lognormal_distribution<>(0, 2), "scale=2");
 484    lognormal_plotter2.plot("Lognormal Distribution PDF (location=0)", "lognormal_pdf2.svg");
 485
 486    distribution_plotter<boost::math::pareto_distribution<> >
 487       pareto_plotter; // Rely on 2nd parameter shape = 1 default.
 488    pareto_plotter.add(boost::math::pareto_distribution<>(1), "scale=1");
 489    pareto_plotter.add(boost::math::pareto_distribution<>(2), "scale=2");
 490    pareto_plotter.add(boost::math::pareto_distribution<>(3), "scale=3");
 491    pareto_plotter.plot("Pareto Distribution PDF (shape=1)", "pareto_pdf1.svg");
 492
 493    distribution_plotter<boost::math::pareto_distribution<> >
 494       pareto_plotter2;
 495    pareto_plotter2.add(boost::math::pareto_distribution<>(1, 0.5), "shape=0.5");
 496    pareto_plotter2.add(boost::math::pareto_distribution<>(1, 1), "shape=1");
 497    pareto_plotter2.add(boost::math::pareto_distribution<>(1, 2), "shape=2");
 498    pareto_plotter2.plot("Pareto Distribution PDF (scale=1)", "pareto_pdf2.svg");
 499
 500    distribution_plotter<boost::math::rayleigh_distribution<> >
 501       rayleigh_plotter;
 502    rayleigh_plotter.add(boost::math::rayleigh_distribution<>(0.5), "&#x3C3;=0.5");
 503    rayleigh_plotter.add(boost::math::rayleigh_distribution<>(1), "&#x3C3;=1");
 504    rayleigh_plotter.add(boost::math::rayleigh_distribution<>(2), "&#x3C3;=2");
 505    rayleigh_plotter.add(boost::math::rayleigh_distribution<>(4), "&#x3C3;=4");
 506    rayleigh_plotter.add(boost::math::rayleigh_distribution<>(10), "&#x3C3;=10");
 507    rayleigh_plotter.plot("Rayleigh Distribution PDF", "rayleigh_pdf.svg");
 508
 509    distribution_plotter<boost::math::rayleigh_distribution<> >
 510       rayleigh_cdf_plotter(false);
 511    rayleigh_cdf_plotter.add(boost::math::rayleigh_distribution<>(0.5), "&#x3C3;=0.5");
 512    rayleigh_cdf_plotter.add(boost::math::rayleigh_distribution<>(1), "&#x3C3;=1");
 513    rayleigh_cdf_plotter.add(boost::math::rayleigh_distribution<>(2), "&#x3C3;=2");
 514    rayleigh_cdf_plotter.add(boost::math::rayleigh_distribution<>(4), "&#x3C3;=4");
 515    rayleigh_cdf_plotter.add(boost::math::rayleigh_distribution<>(10), "&#x3C3;=10");
 516    rayleigh_cdf_plotter.plot("Rayleigh Distribution CDF", "rayleigh_cdf.svg");
 517
 518    distribution_plotter<boost::math::skew_normal_distribution<> >
 519       skew_normal_plotter;
 520    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(0,1,0), "{0,1,0}");
 521    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(0,1,1), "{0,1,1}");
 522    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(0,1,4), "{0,1,4}");
 523    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(0,1,20), "{0,1,20}");
 524    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(0,1,-2), "{0,1,-2}");
 525    skew_normal_plotter.add(boost::math::skew_normal_distribution<>(-2,0.5,-1), "{-2,0.5,-1}");
 526    skew_normal_plotter.plot("Skew Normal Distribution PDF", "skew_normal_pdf.svg");
 527
 528    distribution_plotter<boost::math::skew_normal_distribution<> >
 529       skew_normal_cdf_plotter(false);
 530    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(0,1,0), "{0,1,0}");
 531    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(0,1,1), "{0,1,1}");
 532    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(0,1,4), "{0,1,4}");
 533    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(0,1,20), "{0,1,20}");
 534    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(0,1,-2), "{0,1,-2}");
 535    skew_normal_cdf_plotter.add(boost::math::skew_normal_distribution<>(-2,0.5,-1), "{-2,0.5,-1}");
 536    skew_normal_cdf_plotter.plot("Skew Normal Distribution CDF", "skew_normal_cdf.svg");
 537
 538    distribution_plotter<boost::math::triangular_distribution<> >
 539       triangular_plotter;
 540    triangular_plotter.add(boost::math::triangular_distribution<>(-1,0,1), "{-1,0,1}");
 541    triangular_plotter.add(boost::math::triangular_distribution<>(0,1,1), "{0,1,1}");
 542    triangular_plotter.add(boost::math::triangular_distribution<>(0,1,3), "{0,1,3}");
 543    triangular_plotter.add(boost::math::triangular_distribution<>(0,0.5,1), "{0,0.5,1}");
 544    triangular_plotter.add(boost::math::triangular_distribution<>(-2,0,3), "{-2,0,3}");
 545    triangular_plotter.plot("Triangular Distribution PDF", "triangular_pdf.svg");
 546
 547    distribution_plotter<boost::math::triangular_distribution<> >
 548       triangular_cdf_plotter(false);
 549    triangular_cdf_plotter.add(boost::math::triangular_distribution<>(-1,0,1), "{-1,0,1}");
 550    triangular_cdf_plotter.add(boost::math::triangular_distribution<>(0,1,1), "{0,1,1}");
 551    triangular_cdf_plotter.add(boost::math::triangular_distribution<>(0,1,3), "{0,1,3}");
 552    triangular_cdf_plotter.add(boost::math::triangular_distribution<>(0,0.5,1), "{0,0.5,1}");
 553    triangular_cdf_plotter.add(boost::math::triangular_distribution<>(-2,0,3), "{-2,0,3}");
 554    triangular_cdf_plotter.plot("Triangular Distribution CDF", "triangular_cdf.svg");
 555
 556    distribution_plotter<boost::math::students_t_distribution<> >
 557       students_t_plotter;
 558    students_t_plotter.add(boost::math::students_t_distribution<>(1), "v=1");
 559    students_t_plotter.add(boost::math::students_t_distribution<>(5), "v=5");
 560    students_t_plotter.add(boost::math::students_t_distribution<>(30), "v=30");
 561    students_t_plotter.plot("Students T Distribution PDF", "students_t_pdf.svg");
 562
 563    distribution_plotter<boost::math::weibull_distribution<> >
 564       weibull_plotter;
 565    weibull_plotter.add(boost::math::weibull_distribution<>(0.75), "shape=0.75");
 566    weibull_plotter.add(boost::math::weibull_distribution<>(1), "shape=1");
 567    weibull_plotter.add(boost::math::weibull_distribution<>(5), "shape=5");
 568    weibull_plotter.add(boost::math::weibull_distribution<>(10), "shape=10");
 569    weibull_plotter.plot("Weibull Distribution PDF (scale=1)", "weibull_pdf1.svg");
 570
 571    distribution_plotter<boost::math::weibull_distribution<> >
 572       weibull_plotter2;
 573    weibull_plotter2.add(boost::math::weibull_distribution<>(3, 0.5), "scale=0.5");
 574    weibull_plotter2.add(boost::math::weibull_distribution<>(3, 1), "scale=1");
 575    weibull_plotter2.add(boost::math::weibull_distribution<>(3, 2), "scale=2");
 576    weibull_plotter2.plot("Weibull Distribution PDF (shape=3)", "weibull_pdf2.svg");
 577
 578    distribution_plotter<boost::math::uniform_distribution<> >
 579       uniform_plotter;
 580    uniform_plotter.add(boost::math::uniform_distribution<>(0, 1), "{0,1}");
 581    uniform_plotter.add(boost::math::uniform_distribution<>(0, 3), "{0,3}");
 582    uniform_plotter.add(boost::math::uniform_distribution<>(-2, 3), "{-2,3}");
 583    uniform_plotter.add(boost::math::uniform_distribution<>(-1, 1), "{-1,1}");
 584    uniform_plotter.plot("Uniform Distribution PDF", "uniform_pdf.svg");
 585
 586    distribution_plotter<boost::math::uniform_distribution<> >
 587       uniform_cdf_plotter(false);
 588    uniform_cdf_plotter.add(boost::math::uniform_distribution<>(0, 1), "{0,1}");
 589    uniform_cdf_plotter.add(boost::math::uniform_distribution<>(0, 3), "{0,3}");
 590    uniform_cdf_plotter.add(boost::math::uniform_distribution<>(-2, 3), "{-2,3}");
 591    uniform_cdf_plotter.add(boost::math::uniform_distribution<>(-1, 1), "{-1,1}");
 592    uniform_cdf_plotter.plot("Uniform Distribution CDF", "uniform_cdf.svg");
 593
 594    distribution_plotter<boost::math::bernoulli_distribution<> >
 595       bernoulli_plotter;
 596    bernoulli_plotter.add(boost::math::bernoulli_distribution<>(0.25), "p=0.25");
 597    bernoulli_plotter.add(boost::math::bernoulli_distribution<>(0.5), "p=0.5");
 598    bernoulli_plotter.add(boost::math::bernoulli_distribution<>(0.75), "p=0.75");
 599    bernoulli_plotter.plot("Bernoulli Distribution PDF", "bernoulli_pdf.svg");
 600
 601    distribution_plotter<boost::math::bernoulli_distribution<> >
 602       bernoulli_cdf_plotter(false);
 603    bernoulli_cdf_plotter.add(boost::math::bernoulli_distribution<>(0.25), "p=0.25");
 604    bernoulli_cdf_plotter.add(boost::math::bernoulli_distribution<>(0.5), "p=0.5");
 605    bernoulli_cdf_plotter.add(boost::math::bernoulli_distribution<>(0.75), "p=0.75");
 606    bernoulli_cdf_plotter.plot("Bernoulli Distribution CDF", "bernoulli_cdf.svg");
 607
 608    distribution_plotter<boost::math::binomial_distribution<> >
 609       binomial_plotter;
 610    binomial_plotter.add(boost::math::binomial_distribution<>(5, 0.5), "n=5 p=0.5");
 611    binomial_plotter.add(boost::math::binomial_distribution<>(20, 0.5), "n=20 p=0.5");
 612    binomial_plotter.add(boost::math::binomial_distribution<>(50, 0.5), "n=50 p=0.5");
 613    binomial_plotter.plot("Binomial Distribution PDF", "binomial_pdf_1.svg");
 614
 615    distribution_plotter<boost::math::binomial_distribution<> >
 616       binomial_plotter2;
 617    binomial_plotter2.add(boost::math::binomial_distribution<>(20, 0.1), "n=20 p=0.1");
 618    binomial_plotter2.add(boost::math::binomial_distribution<>(20, 0.5), "n=20 p=0.5");
 619    binomial_plotter2.add(boost::math::binomial_distribution<>(20, 0.9), "n=20 p=0.9");
 620    binomial_plotter2.plot("Binomial Distribution PDF", "binomial_pdf_2.svg");
 621
 622    distribution_plotter<boost::math::negative_binomial_distribution<> >
 623       negative_binomial_plotter;
 624    negative_binomial_plotter.add(boost::math::negative_binomial_distribution<>(20, 0.25), "n=20 p=0.25");
 625    negative_binomial_plotter.add(boost::math::negative_binomial_distribution<>(20, 0.5), "n=20 p=0.5");
 626    negative_binomial_plotter.add(boost::math::negative_binomial_distribution<>(20, 0.75), "n=20 p=0.75");
 627    negative_binomial_plotter.plot("Negative Binomial Distribution PDF", "negative_binomial_pdf_1.svg");
 628
 629    distribution_plotter<boost::math::negative_binomial_distribution<> >
 630       negative_binomial_plotter2;
 631    negative_binomial_plotter2.add(boost::math::negative_binomial_distribution<>(10, 0.5), "n=10 p=0.5");
 632    negative_binomial_plotter2.add(boost::math::negative_binomial_distribution<>(20, 0.5), "n=20 p=0.5");
 633    negative_binomial_plotter2.add(boost::math::negative_binomial_distribution<>(70, 0.5), "n=70 p=0.5");
 634    negative_binomial_plotter2.plot("Negative Binomial Distribution PDF", "negative_binomial_pdf_2.svg");
 635
 636    distribution_plotter<boost::math::poisson_distribution<> >
 637       poisson_plotter;
 638    poisson_plotter.add(boost::math::poisson_distribution<>(5), "&#x3BB;=5");
 639    poisson_plotter.add(boost::math::poisson_distribution<>(10), "&#x3BB;=10");
 640    poisson_plotter.add(boost::math::poisson_distribution<>(20), "&#x3BB;=20");
 641    poisson_plotter.plot("Poisson Distribution PDF", "poisson_pdf_1.svg");
 642
 643    distribution_plotter<boost::math::hypergeometric_distribution<> >
 644       hypergeometric_plotter;
 645    hypergeometric_plotter.add(boost::math::hypergeometric_distribution<>(30, 50, 500), "N=500, r=50, n=30");
 646    hypergeometric_plotter.add(boost::math::hypergeometric_distribution<>(30, 100, 500), "N=500, r=100, n=30");
 647    hypergeometric_plotter.add(boost::math::hypergeometric_distribution<>(30, 250, 500), "N=500, r=250, n=30");
 648    hypergeometric_plotter.add(boost::math::hypergeometric_distribution<>(30, 400, 500), "N=500, r=400, n=30");
 649    hypergeometric_plotter.add(boost::math::hypergeometric_distribution<>(30, 450, 500), "N=500, r=450, n=30");
 650    hypergeometric_plotter.plot("Hypergeometric Distribution PDF", "hypergeometric_pdf_1.svg");
 651
 652    distribution_plotter<boost::math::hypergeometric_distribution<> >
 653       hypergeometric_plotter2;
 654    hypergeometric_plotter2.add(boost::math::hypergeometric_distribution<>(50, 50, 500), "N=500, r=50, n=50");
 655    hypergeometric_plotter2.add(boost::math::hypergeometric_distribution<>(100, 50, 500), "N=500, r=50, n=100");
 656    hypergeometric_plotter2.add(boost::math::hypergeometric_distribution<>(250, 50, 500), "N=500, r=50, n=250");
 657    hypergeometric_plotter2.add(boost::math::hypergeometric_distribution<>(400, 50, 500), "N=500, r=50, n=400");
 658    hypergeometric_plotter2.add(boost::math::hypergeometric_distribution<>(450, 50, 500), "N=500, r=50, n=450");
 659    hypergeometric_plotter2.plot("Hypergeometric Distribution PDF", "hypergeometric_pdf_2.svg");
 660
 661   }
 662   catch (std::exception ex)
 663   {
 664     std::cout << ex.what() << std::endl;
 665   }
 666
 667
 668
 669    /* these graphs for hyperexponential distribution not used.
 670
 671    distribution_plotter<boost::math::hyperexponential_distribution<> >
 672       hyperexponential_plotter;
 673    {
 674        const double probs1_1[] = {1.0};
 675        const double rates1_1[] = {1.0};
 676        hyperexponential_plotter.add(boost::math::hyperexponential_distribution<>(probs1_1,rates1_1), "&#x3B1=(1.0), &#x3BB=(1.0)");
 677        const double probs2_1[] = {0.1,0.9};
 678        const double rates2_1[] = {0.5,1.5};
 679        hyperexponential_plotter.add(boost::math::hyperexponential_distribution<>(probs2_1,rates2_1), "&#x3B1=(0.1,0.9), &#x3BB=(0.5,1.5)");
 680        const double probs2_2[] = {0.9,0.1};
 681        const double rates2_2[] = {0.5,1.5};
 682        hyperexponential_plotter.add(boost::math::hyperexponential_distribution<>(probs2_2,rates2_2), "&#x3B1=(0.9,0.1), &#x3BB=(0.5,1.5)");
 683        const double probs3_1[] = {0.2,0.3,0.5};
 684        const double rates3_1[] = {0.5,1.0,1.5};
 685        hyperexponential_plotter.add(boost::math::hyperexponential_distribution<>(probs3_1,rates3_1), "&#x3B1=(0.2,0.3,0.5), &#x3BB=(0.5,1.0,1.5)");
 686        const double probs3_2[] = {0.5,0.3,0.2};
 687        const double rates3_2[] = {0.5,1.0,1.5};
 688        hyperexponential_plotter.add(boost::math::hyperexponential_distribution<>(probs3_1,rates3_1), "&#x3B1=(0.5,0.3,0.2), &#x3BB=(0.5,1.0,1.5)");
 689    }
 690    hyperexponential_plotter.plot("Hyperexponential Distribution PDF", "hyperexponential_pdf.svg");
 691
 692    distribution_plotter<boost::math::hyperexponential_distribution<> >
 693       hyperexponential_plotter2;
 694    {
 695        const double rates[] = {0.5,1.5};
 696        const double probs1[] = {0.1,0.9};
 697        hyperexponential_plotter2.add(boost::math::hyperexponential_distribution<>(probs1,rates), "&#x3B1=(0.1,0.9), &#x3BB=(0.5,1.5)");
 698        const double probs2[] = {0.6,0.4};
 699        hyperexponential_plotter2.add(boost::math::hyperexponential_distribution<>(probs2,rates), "&#x3B1=(0.6,0.4), &#x3BB=(0.5,1.5)");
 700        const double probs3[] = {0.9,0.1};
 701        hyperexponential_plotter2.add(boost::math::hyperexponential_distribution<>(probs3,rates), "&#x3B1=(0.9,0.1), &#x3BB=(0.5,1.5)");
 702    }
 703    hyperexponential_plotter2.plot("Hyperexponential Distribution PDF (Different Probabilities, Same Rates)", "hyperexponential_pdf_samerate.svg");
 704
 705    distribution_plotter<boost::math::hyperexponential_distribution<> >
 706       hyperexponential_plotter3;
 707    {
 708        const double probs1[] = {1.0};
 709        const double rates1[] = {2.0};
 710        hyperexponential_plotter3.add(boost::math::hyperexponential_distribution<>(probs1,rates1), "&#x3B1=(1.0), &#x3BB=(2.0)");
 711        const double probs2[] = {0.5,0.5};
 712        const double rates2[] = {0.3,1.5};
 713        hyperexponential_plotter3.add(boost::math::hyperexponential_distribution<>(probs2,rates2), "&#x3B1=(0.5,0.5), &#x3BB=(0.3,1.5)");
 714        const double probs3[] = {1.0/3.0,1.0/3.0,1.0/3.0};
 715        const double rates3[] = {0.2,1.5,3.0};
 716        hyperexponential_plotter3.add(boost::math::hyperexponential_distribution<>(probs2,rates2), "&#x3B1=(1.0/3.0,1.0/3.0,1.0/3.0), &#x3BB=(0.2,1.5,3.0)");
 717    }
 718    hyperexponential_plotter3.plot("Hyperexponential Distribution PDF (Different Number of Phases, Same Mean)", "hyperexponential_pdf_samemean.svg");
 719    */
 720
 721
 722 } // int main()