math/distributions/triangular.hpp

0001 //  Copyright John Maddock 2006, 2007.
0002 //  Copyright Paul A. Bristow 2006, 2007.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_STATS_TRIANGULAR_HPP
0008 #define BOOST_STATS_TRIANGULAR_HPP
0009
0010 // http://mathworld.wolfram.com/TriangularDistribution.html
0011 // Note that the 'constructors' defined by Wolfram are difference from those here,
0012 // for example
0013 // N[variance[triangulardistribution{1, +2}, 1.5], 50] computes
0014 // 0.041666666666666666666666666666666666666666666666667
0015 // TriangularDistribution{1, +2}, 1.5 is the analog of triangular_distribution(1, 1.5, 2)
0016
0017 // http://en.wikipedia.org/wiki/Triangular_distribution
0018
0019 #include <boost/math/distributions/fwd.hpp>
0020 #include <boost/math/special_functions/expm1.hpp>
0021 #include <boost/math/distributions/detail/common_error_handling.hpp>
0022 #include <boost/math/distributions/complement.hpp>
0023 #include <boost/math/constants/constants.hpp>
0024
0025 #include <utility>
0026
0027 namespace boost{ namespace math
0028 {
0029   namespace detail
0030   {
0031     template <class RealType, class Policy>
0032     inline bool check_triangular_lower(
0033       const char* function,
0034       RealType lower,
0035       RealType* result, const Policy& pol)
0036     {
0037       if((boost::math::isfinite)(lower))
0038       { // Any finite value is OK.
0039         return true;
0040       }
0041       else
0042       { // Not finite: infinity or NaN.
0043         *result = policies::raise_domain_error<RealType>(
0044           function,
0045           "Lower parameter is %1%, but must be finite!", lower, pol);
0046         return false;
0047       }
0048     } // bool check_triangular_lower(
0049
0050     template <class RealType, class Policy>
0051     inline bool check_triangular_mode(
0052       const char* function,
0053       RealType mode,
0054       RealType* result, const Policy& pol)
0055     {
0056       if((boost::math::isfinite)(mode))
0057       { // any finite value is OK.
0058         return true;
0059       }
0060       else
0061       { // Not finite: infinity or NaN.
0062         *result = policies::raise_domain_error<RealType>(
0063           function,
0064           "Mode parameter is %1%, but must be finite!", mode, pol);
0065         return false;
0066       }
0067     } // bool check_triangular_mode(
0068
0069     template <class RealType, class Policy>
0070     inline bool check_triangular_upper(
0071       const char* function,
0072       RealType upper,
0073       RealType* result, const Policy& pol)
0074     {
0075       if((boost::math::isfinite)(upper))
0076       { // any finite value is OK.
0077         return true;
0078       }
0079       else
0080       { // Not finite: infinity or NaN.
0081         *result = policies::raise_domain_error<RealType>(
0082           function,
0083           "Upper parameter is %1%, but must be finite!", upper, pol);
0084         return false;
0085       }
0086     } // bool check_triangular_upper(
0087
0088     template <class RealType, class Policy>
0089     inline bool check_triangular_x(
0090       const char* function,
0091       RealType const& x,
0092       RealType* result, const Policy& pol)
0093     {
0094       if((boost::math::isfinite)(x))
0095       { // Any finite value is OK
0096         return true;
0097       }
0098       else
0099       { // Not finite: infinity or NaN.
0100         *result = policies::raise_domain_error<RealType>(
0101           function,
0102           "x parameter is %1%, but must be finite!", x, pol);
0103         return false;
0104       }
0105     } // bool check_triangular_x
0106
0107     template <class RealType, class Policy>
0108     inline bool check_triangular(
0109       const char* function,
0110       RealType lower,
0111       RealType mode,
0112       RealType upper,
0113       RealType* result, const Policy& pol)
0114     {
0115       if ((check_triangular_lower(function, lower, result, pol) == false)
0116         || (check_triangular_mode(function, mode, result, pol) == false)
0117         || (check_triangular_upper(function, upper, result, pol) == false))
0118       { // Some parameter not finite.
0119         return false;
0120       }
0121       else if (lower >= upper) // lower == upper NOT useful.
0122       { // lower >= upper.
0123         *result = policies::raise_domain_error<RealType>(
0124           function,
0125           "lower parameter is %1%, but must be less than upper!", lower, pol);
0126         return false;
0127       }
0128       else
0129       { // Check lower <= mode <= upper.
0130         if (mode < lower)
0131         {
0132           *result = policies::raise_domain_error<RealType>(
0133             function,
0134             "mode parameter is %1%, but must be >= than lower!", lower, pol);
0135           return false;
0136         }
0137         if (mode > upper)
0138         {
0139           *result = policies::raise_domain_error<RealType>(
0140             function,
0141             "mode parameter is %1%, but must be <= than upper!", upper, pol);
0142           return false;
0143         }
0144         return true; // All OK.
0145       }
0146     } // bool check_triangular
0147   } // namespace detail
0148
0149   template <class RealType = double, class Policy = policies::policy<> >
0150   class triangular_distribution
0151   {
0152   public:
0153     typedef RealType value_type;
0154     typedef Policy policy_type;
0155
0156     triangular_distribution(RealType l_lower = -1, RealType l_mode = 0, RealType l_upper = 1)
0157       : m_lower(l_lower), m_mode(l_mode), m_upper(l_upper) // Constructor.
0158     { // Evans says 'standard triangular' is lower 0, mode 1/2, upper 1,
0159       // has median sqrt(c/2) for c <=1/2 and 1 - sqrt(1-c)/2 for c >= 1/2
0160       // But this -1, 0, 1 is more useful in most applications to approximate normal distribution,
0161       // where the central value is the most likely and deviations either side equally likely.
0162       RealType result;
0163       detail::check_triangular("boost::math::triangular_distribution<%1%>::triangular_distribution",l_lower, l_mode, l_upper, &result, Policy());
0164     }
0165     // Accessor functions.
0166     RealType lower()const
0167     {
0168       return m_lower;
0169     }
0170     RealType mode()const
0171     {
0172       return m_mode;
0173     }
0174     RealType upper()const
0175     {
0176       return m_upper;
0177     }
0178   private:
0179     // Data members:
0180     RealType m_lower;  // distribution lower aka a
0181     RealType m_mode;  // distribution mode aka c
0182     RealType m_upper;  // distribution upper aka b
0183   }; // class triangular_distribution
0184
0185   typedef triangular_distribution<double> triangular;
0186
0187   #ifdef __cpp_deduction_guides
0188   template <class RealType>
0189   triangular_distribution(RealType)->triangular_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0190   template <class RealType>
0191   triangular_distribution(RealType,RealType)->triangular_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0192   template <class RealType>
0193   triangular_distribution(RealType,RealType,RealType)->triangular_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0194   #endif
0195
0196   template <class RealType, class Policy>
0197   inline const std::pair<RealType, RealType> range(const triangular_distribution<RealType, Policy>& /* dist */)
0198   { // Range of permissible values for random variable x.
0199     using boost::math::tools::max_value;
0200     return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
0201   }
0202
0203   template <class RealType, class Policy>
0204   inline const std::pair<RealType, RealType> support(const triangular_distribution<RealType, Policy>& dist)
0205   { // Range of supported values for random variable x.
0206     // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0207     return std::pair<RealType, RealType>(dist.lower(), dist.upper());
0208   }
0209
0210   template <class RealType, class Policy>
0211   RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
0212   {
0213     static const char* function = "boost::math::pdf(const triangular_distribution<%1%>&, %1%)";
0214     RealType lower = dist.lower();
0215     RealType mode = dist.mode();
0216     RealType upper = dist.upper();
0217     RealType result = 0; // of checks.
0218     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
0219     {
0220       return result;
0221     }
0222     if(false == detail::check_triangular_x(function, x, &result, Policy()))
0223     {
0224       return result;
0225     }
0226     if((x < lower) || (x > upper))
0227     {
0228       return 0;
0229     }
0230     if (x == lower)
0231     { // (mode - lower) == 0 which would lead to divide by zero!
0232       return (mode == lower) ? 2 / (upper - lower) : RealType(0);
0233     }
0234     else if (x == upper)
0235     {
0236       return (mode == upper) ? 2 / (upper - lower) : RealType(0);
0237     }
0238     else if (x <= mode)
0239     {
0240       return 2 * (x - lower) / ((upper - lower) * (mode - lower));
0241     }
0242     else
0243     {  // (x > mode)
0244       return 2 * (upper - x) / ((upper - lower) * (upper - mode));
0245     }
0246   } // RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
0247
0248   template <class RealType, class Policy>
0249   inline RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
0250   {
0251     static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
0252     RealType lower = dist.lower();
0253     RealType mode = dist.mode();
0254     RealType upper = dist.upper();
0255     RealType result = 0; // of checks.
0256     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
0257     {
0258       return result;
0259     }
0260     if(false == detail::check_triangular_x(function, x, &result, Policy()))
0261     {
0262       return result;
0263     }
0264     if((x <= lower))
0265     {
0266       return 0;
0267     }
0268     if (x >= upper)
0269     {
0270       return 1;
0271     }
0272     // else lower < x < upper
0273     if (x <= mode)
0274     {
0275       return ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
0276     }
0277     else
0278     {
0279       return 1 - (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
0280     }
0281   } // RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
0282
0283   template <class RealType, class Policy>
0284   RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& p)
0285   {
0286     BOOST_MATH_STD_USING  // for ADL of std functions (sqrt).
0287     static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
0288     RealType lower = dist.lower();
0289     RealType mode = dist.mode();
0290     RealType upper = dist.upper();
0291     RealType result = 0; // of checks
0292     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
0293     {
0294       return result;
0295     }
0296     if(false == detail::check_probability(function, p, &result, Policy()))
0297     {
0298       return result;
0299     }
0300     if(p == 0)
0301     {
0302       return lower;
0303     }
0304     if(p == 1)
0305     {
0306       return upper;
0307     }
0308     RealType p0 = (mode - lower) / (upper - lower);
0309     RealType q = 1 - p;
0310     if (p < p0)
0311     {
0312       result = sqrt((upper - lower) * (mode - lower) * p) + lower;
0313     }
0314     else if (p == p0)
0315     {
0316       result = mode;
0317     }
0318     else // p > p0
0319     {
0320       result = upper - sqrt((upper - lower) * (upper - mode) * q);
0321     }
0322     return result;
0323
0324   } // RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& q)
0325
0326   template <class RealType, class Policy>
0327   RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
0328   {
0329     static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
0330     RealType lower = c.dist.lower();
0331     RealType mode = c.dist.mode();
0332     RealType upper = c.dist.upper();
0333     RealType x = c.param;
0334     RealType result = 0; // of checks.
0335     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
0336     {
0337       return result;
0338     }
0339     if(false == detail::check_triangular_x(function, x, &result, Policy()))
0340     {
0341       return result;
0342     }
0343     if (x <= lower)
0344     {
0345       return 1;
0346     }
0347     if (x >= upper)
0348     {
0349       return 0;
0350     }
0351     if (x <= mode)
0352     {
0353       return 1 - ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
0354     }
0355     else
0356     {
0357       return (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
0358     }
0359   } // RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
0360
0361   template <class RealType, class Policy>
0362   RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
0363   {
0364     BOOST_MATH_STD_USING  // Aid ADL for sqrt.
0365     static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
0366     RealType l = c.dist.lower();
0367     RealType m = c.dist.mode();
0368     RealType u = c.dist.upper();
0369     RealType q = c.param; // probability 0 to 1.
0370     RealType result = 0; // of checks.
0371     if(false == detail::check_triangular(function, l, m, u, &result, Policy()))
0372     {
0373       return result;
0374     }
0375     if(false == detail::check_probability(function, q, &result, Policy()))
0376     {
0377       return result;
0378     }
0379     if(q == 0)
0380     {
0381       return u;
0382     }
0383     if(q == 1)
0384     {
0385       return l;
0386     }
0387     RealType lower = c.dist.lower();
0388     RealType mode = c.dist.mode();
0389     RealType upper = c.dist.upper();
0390
0391     RealType p = 1 - q;
0392     RealType p0 = (mode - lower) / (upper - lower);
0393     if(p < p0)
0394     {
0395       RealType s = (upper - lower) * (mode - lower);
0396       s *= p;
0397       result = sqrt((upper - lower) * (mode - lower) * p) + lower;
0398     }
0399     else if (p == p0)
0400     {
0401       result = mode;
0402     }
0403     else // p > p0
0404     {
0405       result = upper - sqrt((upper - lower) * (upper - mode) * q);
0406     }
0407     return result;
0408   } // RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
0409
0410   template <class RealType, class Policy>
0411   inline RealType mean(const triangular_distribution<RealType, Policy>& dist)
0412   {
0413     static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
0414     RealType lower = dist.lower();
0415     RealType mode = dist.mode();
0416     RealType upper = dist.upper();
0417     RealType result = 0;  // of checks.
0418     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
0419     {
0420       return result;
0421     }
0422     return (lower + upper + mode) / 3;
0423   } // RealType mean(const triangular_distribution<RealType, Policy>& dist)
0424
0425
0426   template <class RealType, class Policy>
0427   inline RealType variance(const triangular_distribution<RealType, Policy>& dist)
0428   {
0429     static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
0430     RealType lower = dist.lower();
0431     RealType mode = dist.mode();
0432     RealType upper = dist.upper();
0433     RealType result = 0; // of checks.
0434     if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
0435     {
0436       return result;
0437     }
0438     return (lower * lower + upper * upper + mode * mode - lower * upper - lower * mode - upper * mode) / 18;
0439   } // RealType variance(const triangular_distribution<RealType, Policy>& dist)
0440
0441   template <class RealType, class Policy>
0442   inline RealType mode(const triangular_distribution<RealType, Policy>& dist)
0443   {
0444     static const char* function = "boost::math::mode(const triangular_distribution<%1%>&)";
0445     RealType mode = dist.mode();
0446     RealType result = 0; // of checks.
0447     if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
0448     { // This should never happen!
0449       return result;
0450     }
0451     return mode;
0452   } // RealType mode
0453
0454   template <class RealType, class Policy>
0455   inline RealType median(const triangular_distribution<RealType, Policy>& dist)
0456   {
0457     BOOST_MATH_STD_USING // ADL of std functions.
0458     static const char* function = "boost::math::median(const triangular_distribution<%1%>&)";
0459     RealType mode = dist.mode();
0460     RealType result = 0; // of checks.
0461     if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
0462     { // This should never happen!
0463       return result;
0464     }
0465     RealType lower = dist.lower();
0466     RealType upper = dist.upper();
0467     if (mode >= (upper + lower) / 2)
0468     {
0469       return lower + sqrt((upper - lower) * (mode - lower)) / constants::root_two<RealType>();
0470     }
0471     else
0472     {
0473       return upper - sqrt((upper - lower) * (upper - mode)) / constants::root_two<RealType>();
0474     }
0475   } // RealType mode
0476
0477   template <class RealType, class Policy>
0478   inline RealType skewness(const triangular_distribution<RealType, Policy>& dist)
0479   {
0480     BOOST_MATH_STD_USING  // for ADL of std functions
0481     using namespace boost::math::constants; // for root_two
0482     static const char* function = "boost::math::skewness(const triangular_distribution<%1%>&)";
0483
0484     RealType lower = dist.lower();
0485     RealType mode = dist.mode();
0486     RealType upper = dist.upper();
0487     RealType result = 0; // of checks.
0488     if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy()))
0489     {
0490       return result;
0491     }
0492     return root_two<RealType>() * (lower + upper - 2 * mode) * (2 * lower - upper - mode) * (lower - 2 * upper + mode) /
0493       (5 * pow((lower * lower + upper * upper + mode * mode
0494         - lower * upper - lower * mode - upper * mode), RealType(3)/RealType(2)));
0495     // #11768: Skewness formula for triangular distribution is incorrect -  corrected 29 Oct 2015 for release 1.61.
0496   } // RealType skewness(const triangular_distribution<RealType, Policy>& dist)
0497
0498   template <class RealType, class Policy>
0499   inline RealType kurtosis(const triangular_distribution<RealType, Policy>& dist)
0500   { // These checks may be belt and braces as should have been checked on construction?
0501     static const char* function = "boost::math::kurtosis(const triangular_distribution<%1%>&)";
0502     RealType lower = dist.lower();
0503     RealType upper = dist.upper();
0504     RealType mode = dist.mode();
0505     RealType result = 0;  // of checks.
0506     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
0507     {
0508       return result;
0509     }
0510     return static_cast<RealType>(12)/5; //  12/5 = 2.4;
0511   } // RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
0512
0513   template <class RealType, class Policy>
0514   inline RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
0515   { // These checks may be belt and braces as should have been checked on construction?
0516     static const char* function = "boost::math::kurtosis_excess(const triangular_distribution<%1%>&)";
0517     RealType lower = dist.lower();
0518     RealType upper = dist.upper();
0519     RealType mode = dist.mode();
0520     RealType result = 0;  // of checks.
0521     if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
0522     {
0523       return result;
0524     }
0525     return static_cast<RealType>(-3)/5; // - 3/5 = -0.6
0526     // Assuming mathworld really means kurtosis excess?  Wikipedia now corrected to match this.
0527   }
0528
0529   template <class RealType, class Policy>
0530   inline RealType entropy(const triangular_distribution<RealType, Policy>& dist)
0531   {
0532     using std::log;
0533     return constants::half<RealType>() + log((dist.upper() - dist.lower())/2);
0534   }
0535
0536 } // namespace math
0537 } // namespace boost
0538
0539 // This include must be at the end, *after* the accessors
0540 // for this distribution have been defined, in order to
0541 // keep compilers that support two-phase lookup happy.
0542 #include <boost/math/distributions/detail/derived_accessors.hpp>
0543
0544 #endif // BOOST_STATS_TRIANGULAR_HPP
0545
0546
0547