ceph/src/boost/libs/math/include/boost/math/distributions/triangular.hpp

   1 //  Copyright John Maddock 2006, 2007.
   2 //  Copyright Paul A. Bristow 2006, 2007.
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_STATS_TRIANGULAR_HPP
   8 #define BOOST_STATS_TRIANGULAR_HPP
   9
  10 // http://mathworld.wolfram.com/TriangularDistribution.html
  11 // Note that the 'constructors' defined by Wolfram are difference from those here,
  12 // for example
  13 // N[variance[triangulardistribution{1, +2}, 1.5], 50] computes
  14 // 0.041666666666666666666666666666666666666666666666667
  15 // TriangularDistribution{1, +2}, 1.5 is the analog of triangular_distribution(1, 1.5, 2)
  16
  17 // http://en.wikipedia.org/wiki/Triangular_distribution
  18
  19 #include <boost/math/distributions/fwd.hpp>
  20 #include <boost/math/special_functions/expm1.hpp>
  21 #include <boost/math/distributions/detail/common_error_handling.hpp>
  22 #include <boost/math/distributions/complement.hpp>
  23 #include <boost/math/constants/constants.hpp>
  24
  25 #include <utility>
  26
  27 namespace boost{ namespace math
  28 {
  29   namespace detail
  30   {
  31     template <class RealType, class Policy>
  32     inline bool check_triangular_lower(
  33       const char* function,
  34       RealType lower,
  35       RealType* result, const Policy& pol)
  36     {
  37       if((boost::math::isfinite)(lower))
  38       { // Any finite value is OK.
  39         return true;
  40       }
  41       else
  42       { // Not finite: infinity or NaN.
  43         *result = policies::raise_domain_error<RealType>(
  44           function,
  45           "Lower parameter is %1%, but must be finite!", lower, pol);
  46         return false;
  47       }
  48     } // bool check_triangular_lower(
  49
  50     template <class RealType, class Policy>
  51     inline bool check_triangular_mode(
  52       const char* function,
  53       RealType mode,
  54       RealType* result, const Policy& pol)
  55     {
  56       if((boost::math::isfinite)(mode))
  57       { // any finite value is OK.
  58         return true;
  59       }
  60       else
  61       { // Not finite: infinity or NaN.
  62         *result = policies::raise_domain_error<RealType>(
  63           function,
  64           "Mode parameter is %1%, but must be finite!", mode, pol);
  65         return false;
  66       }
  67     } // bool check_triangular_mode(
  68
  69     template <class RealType, class Policy>
  70     inline bool check_triangular_upper(
  71       const char* function,
  72       RealType upper,
  73       RealType* result, const Policy& pol)
  74     {
  75       if((boost::math::isfinite)(upper))
  76       { // any finite value is OK.
  77         return true;
  78       }
  79       else
  80       { // Not finite: infinity or NaN.
  81         *result = policies::raise_domain_error<RealType>(
  82           function,
  83           "Upper parameter is %1%, but must be finite!", upper, pol);
  84         return false;
  85       }
  86     } // bool check_triangular_upper(
  87
  88     template <class RealType, class Policy>
  89     inline bool check_triangular_x(
  90       const char* function,
  91       RealType const& x,
  92       RealType* result, const Policy& pol)
  93     {
  94       if((boost::math::isfinite)(x))
  95       { // Any finite value is OK
  96         return true;
  97       }
  98       else
  99       { // Not finite: infinity or NaN.
 100         *result = policies::raise_domain_error<RealType>(
 101           function,
 102           "x parameter is %1%, but must be finite!", x, pol);
 103         return false;
 104       }
 105     } // bool check_triangular_x
 106
 107     template <class RealType, class Policy>
 108     inline bool check_triangular(
 109       const char* function,
 110       RealType lower,
 111       RealType mode,
 112       RealType upper,
 113       RealType* result, const Policy& pol)
 114     {
 115       if ((check_triangular_lower(function, lower, result, pol) == false)
 116         || (check_triangular_mode(function, mode, result, pol) == false)
 117         || (check_triangular_upper(function, upper, result, pol) == false))
 118       { // Some parameter not finite.
 119         return false;
 120       }
 121       else if (lower >= upper) // lower == upper NOT useful.
 122       { // lower >= upper.
 123         *result = policies::raise_domain_error<RealType>(
 124           function,
 125           "lower parameter is %1%, but must be less than upper!", lower, pol);
 126         return false;
 127       }
 128       else
 129       { // Check lower <= mode <= upper.
 130         if (mode < lower)
 131         {
 132           *result = policies::raise_domain_error<RealType>(
 133             function,
 134             "mode parameter is %1%, but must be >= than lower!", lower, pol);
 135           return false;
 136         }
 137         if (mode > upper)
 138         {
 139           *result = policies::raise_domain_error<RealType>(
 140             function,
 141             "mode parameter is %1%, but must be <= than upper!", upper, pol);
 142           return false;
 143         }
 144         return true; // All OK.
 145       }
 146     } // bool check_triangular
 147   } // namespace detail
 148
 149   template <class RealType = double, class Policy = policies::policy<> >
 150   class triangular_distribution
 151   {
 152   public:
 153     typedef RealType value_type;
 154     typedef Policy policy_type;
 155
 156     triangular_distribution(RealType l_lower = -1, RealType l_mode = 0, RealType l_upper = 1)
 157       : m_lower(l_lower), m_mode(l_mode), m_upper(l_upper) // Constructor.
 158     { // Evans says 'standard triangular' is lower 0, mode 1/2, upper 1,
 159       // has median sqrt(c/2) for c <=1/2 and 1 - sqrt(1-c)/2 for c >= 1/2
 160       // But this -1, 0, 1 is more useful in most applications to approximate normal distribution,
 161       // where the central value is the most likely and deviations either side equally likely.
 162       RealType result;
 163       detail::check_triangular("boost::math::triangular_distribution<%1%>::triangular_distribution",l_lower, l_mode, l_upper, &result, Policy());
 164     }
 165     // Accessor functions.
 166     RealType lower()const
 167     {
 168       return m_lower;
 169     }
 170     RealType mode()const
 171     {
 172       return m_mode;
 173     }
 174     RealType upper()const
 175     {
 176       return m_upper;
 177     }
 178   private:
 179     // Data members:
 180     RealType m_lower;  // distribution lower aka a
 181     RealType m_mode;  // distribution mode aka c
 182     RealType m_upper;  // distribution upper aka b
 183   }; // class triangular_distribution
 184
 185   typedef triangular_distribution<double> triangular;
 186
 187   template <class RealType, class Policy>
 188   inline const std::pair<RealType, RealType> range(const triangular_distribution<RealType, Policy>& /* dist */)
 189   { // Range of permissible values for random variable x.
 190     using boost::math::tools::max_value;
 191     return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
 192   }
 193
 194   template <class RealType, class Policy>
 195   inline const std::pair<RealType, RealType> support(const triangular_distribution<RealType, Policy>& dist)
 196   { // Range of supported values for random variable x.
 197     // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 198     return std::pair<RealType, RealType>(dist.lower(), dist.upper());
 199   }
 200
 201   template <class RealType, class Policy>
 202   RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
 203   {
 204     static const char* function = "boost::math::pdf(const triangular_distribution<%1%>&, %1%)";
 205     RealType lower = dist.lower();
 206     RealType mode = dist.mode();
 207     RealType upper = dist.upper();
 208     RealType result = 0; // of checks.
 209     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
 210     {
 211       return result;
 212     }
 213     if(false == detail::check_triangular_x(function, x, &result, Policy()))
 214     {
 215       return result;
 216     }
 217     if((x < lower) || (x > upper))
 218     {
 219       return 0;
 220     }
 221     if (x == lower)
 222     { // (mode - lower) == 0 which would lead to divide by zero!
 223       return (mode == lower) ? 2 / (upper - lower) : RealType(0);
 224     }
 225     else if (x == upper)
 226     {
 227       return (mode == upper) ? 2 / (upper - lower) : RealType(0);
 228     }
 229     else if (x <= mode)
 230     {
 231       return 2 * (x - lower) / ((upper - lower) * (mode - lower));
 232     }
 233     else
 234     {  // (x > mode)
 235       return 2 * (upper - x) / ((upper - lower) * (upper - mode));
 236     }
 237   } // RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
 238
 239   template <class RealType, class Policy>
 240   inline RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
 241   {
 242     static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
 243     RealType lower = dist.lower();
 244     RealType mode = dist.mode();
 245     RealType upper = dist.upper();
 246     RealType result = 0; // of checks.
 247     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
 248     {
 249       return result;
 250     }
 251     if(false == detail::check_triangular_x(function, x, &result, Policy()))
 252     {
 253       return result;
 254     }
 255     if((x <= lower))
 256     {
 257       return 0;
 258     }
 259     if (x >= upper)
 260     {
 261       return 1;
 262     }
 263     // else lower < x < upper
 264     if (x <= mode)
 265     {
 266       return ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
 267     }
 268     else
 269     {
 270       return 1 - (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
 271     }
 272   } // RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
 273
 274   template <class RealType, class Policy>
 275   RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& p)
 276   {
 277     BOOST_MATH_STD_USING  // for ADL of std functions (sqrt).
 278     static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
 279     RealType lower = dist.lower();
 280     RealType mode = dist.mode();
 281     RealType upper = dist.upper();
 282     RealType result = 0; // of checks
 283     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
 284     {
 285       return result;
 286     }
 287     if(false == detail::check_probability(function, p, &result, Policy()))
 288     {
 289       return result;
 290     }
 291     if(p == 0)
 292     {
 293       return lower;
 294     }
 295     if(p == 1)
 296     {
 297       return upper;
 298     }
 299     RealType p0 = (mode - lower) / (upper - lower);
 300     RealType q = 1 - p;
 301     if (p < p0)
 302     {
 303       result = sqrt((upper - lower) * (mode - lower) * p) + lower;
 304     }
 305     else if (p == p0)
 306     {
 307       result = mode;
 308     }
 309     else // p > p0
 310     {
 311       result = upper - sqrt((upper - lower) * (upper - mode) * q);
 312     }
 313     return result;
 314
 315   } // RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& q)
 316
 317   template <class RealType, class Policy>
 318   RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
 319   {
 320     static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
 321     RealType lower = c.dist.lower();
 322     RealType mode = c.dist.mode();
 323     RealType upper = c.dist.upper();
 324     RealType x = c.param;
 325     RealType result = 0; // of checks.
 326     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
 327     {
 328       return result;
 329     }
 330     if(false == detail::check_triangular_x(function, x, &result, Policy()))
 331     {
 332       return result;
 333     }
 334     if (x <= lower)
 335     {
 336       return 1;
 337     }
 338     if (x >= upper)
 339     {
 340       return 0;
 341     }
 342     if (x <= mode)
 343     {
 344       return 1 - ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
 345     }
 346     else
 347     {
 348       return (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
 349     }
 350   } // RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
 351
 352   template <class RealType, class Policy>
 353   RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
 354   {
 355     BOOST_MATH_STD_USING  // Aid ADL for sqrt.
 356     static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
 357     RealType l = c.dist.lower();
 358     RealType m = c.dist.mode();
 359     RealType u = c.dist.upper();
 360     RealType q = c.param; // probability 0 to 1.
 361     RealType result = 0; // of checks.
 362     if(false == detail::check_triangular(function, l, m, u, &result, Policy()))
 363     {
 364       return result;
 365     }
 366     if(false == detail::check_probability(function, q, &result, Policy()))
 367     {
 368       return result;
 369     }
 370     if(q == 0)
 371     {
 372       return u;
 373     }
 374     if(q == 1)
 375     {
 376       return l;
 377     }
 378     RealType lower = c.dist.lower();
 379     RealType mode = c.dist.mode();
 380     RealType upper = c.dist.upper();
 381
 382     RealType p = 1 - q;
 383     RealType p0 = (mode - lower) / (upper - lower);
 384     if(p < p0)
 385     {
 386       RealType s = (upper - lower) * (mode - lower);
 387       s *= p;
 388       result = sqrt((upper - lower) * (mode - lower) * p) + lower;
 389     }
 390     else if (p == p0)
 391     {
 392       result = mode;
 393     }
 394     else // p > p0
 395     {
 396       result = upper - sqrt((upper - lower) * (upper - mode) * q);
 397     }
 398     return result;
 399   } // RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
 400
 401   template <class RealType, class Policy>
 402   inline RealType mean(const triangular_distribution<RealType, Policy>& dist)
 403   {
 404     static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
 405     RealType lower = dist.lower();
 406     RealType mode = dist.mode();
 407     RealType upper = dist.upper();
 408     RealType result = 0;  // of checks.
 409     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
 410     {
 411       return result;
 412     }
 413     return (lower + upper + mode) / 3;
 414   } // RealType mean(const triangular_distribution<RealType, Policy>& dist)
 415
 416
 417   template <class RealType, class Policy>
 418   inline RealType variance(const triangular_distribution<RealType, Policy>& dist)
 419   {
 420     static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
 421     RealType lower = dist.lower();
 422     RealType mode = dist.mode();
 423     RealType upper = dist.upper();
 424     RealType result = 0; // of checks.
 425     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
 426     {
 427       return result;
 428     }
 429     return (lower * lower + upper * upper + mode * mode - lower * upper - lower * mode - upper * mode) / 18;
 430   } // RealType variance(const triangular_distribution<RealType, Policy>& dist)
 431
 432   template <class RealType, class Policy>
 433   inline RealType mode(const triangular_distribution<RealType, Policy>& dist)
 434   {
 435     static const char* function = "boost::math::mode(const triangular_distribution<%1%>&)";
 436     RealType mode = dist.mode();
 437     RealType result = 0; // of checks.
 438     if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
 439     { // This should never happen!
 440       return result;
 441     }
 442     return mode;
 443   } // RealType mode
 444
 445   template <class RealType, class Policy>
 446   inline RealType median(const triangular_distribution<RealType, Policy>& dist)
 447   {
 448     BOOST_MATH_STD_USING // ADL of std functions.
 449     static const char* function = "boost::math::median(const triangular_distribution<%1%>&)";
 450     RealType mode = dist.mode();
 451     RealType result = 0; // of checks.
 452     if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
 453     { // This should never happen!
 454       return result;
 455     }
 456     RealType lower = dist.lower();
 457     RealType upper = dist.upper();
 458     if (mode >= (upper + lower) / 2)
 459     {
 460       return lower + sqrt((upper - lower) * (mode - lower)) / constants::root_two<RealType>();
 461     }
 462     else
 463     {
 464       return upper - sqrt((upper - lower) * (upper - mode)) / constants::root_two<RealType>();
 465     }
 466   } // RealType mode
 467
 468   template <class RealType, class Policy>
 469   inline RealType skewness(const triangular_distribution<RealType, Policy>& dist)
 470   {
 471     BOOST_MATH_STD_USING  // for ADL of std functions
 472     using namespace boost::math::constants; // for root_two
 473     static const char* function = "boost::math::skewness(const triangular_distribution<%1%>&)";
 474
 475     RealType lower = dist.lower();
 476     RealType mode = dist.mode();
 477     RealType upper = dist.upper();
 478     RealType result = 0; // of checks.
 479     if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy()))
 480     {
 481       return result;
 482     }
 483     return root_two<RealType>() * (lower + upper - 2 * mode) * (2 * lower - upper - mode) * (lower - 2 * upper + mode) /
 484       (5 * pow((lower * lower + upper * upper + mode * mode
 485         - lower * upper - lower * mode - upper * mode), RealType(3)/RealType(2)));
 486     // #11768: Skewness formula for triangular distribution is incorrect -  corrected 29 Oct 2015 for release 1.61.
 487   } // RealType skewness(const triangular_distribution<RealType, Policy>& dist)
 488
 489   template <class RealType, class Policy>
 490   inline RealType kurtosis(const triangular_distribution<RealType, Policy>& dist)
 491   { // These checks may be belt and braces as should have been checked on construction?
 492     static const char* function = "boost::math::kurtosis(const triangular_distribution<%1%>&)";
 493     RealType lower = dist.lower();
 494     RealType upper = dist.upper();
 495     RealType mode = dist.mode();
 496     RealType result = 0;  // of checks.
 497     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
 498     {
 499       return result;
 500     }
 501     return static_cast<RealType>(12)/5; //  12/5 = 2.4;
 502   } // RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
 503
 504   template <class RealType, class Policy>
 505   inline RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
 506   { // These checks may be belt and braces as should have been checked on construction?
 507     static const char* function = "boost::math::kurtosis_excess(const triangular_distribution<%1%>&)";
 508     RealType lower = dist.lower();
 509     RealType upper = dist.upper();
 510     RealType mode = dist.mode();
 511     RealType result = 0;  // of checks.
 512     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
 513     {
 514       return result;
 515     }
 516     return static_cast<RealType>(-3)/5; // - 3/5 = -0.6
 517     // Assuming mathworld really means kurtosis excess?  Wikipedia now corrected to match this.
 518   }
 519
 520 } // namespace math
 521 } // namespace boost
 522
 523 // This include must be at the end, *after* the accessors
 524 // for this distribution have been defined, in order to
 525 // keep compilers that support two-phase lookup happy.
 526 #include <boost/math/distributions/detail/derived_accessors.hpp>
 527
 528 #endif // BOOST_STATS_TRIANGULAR_HPP
 529
 530
 531