// Copyright John Maddock 2006, 2007. | |
// Copyright Paul A. Bristow 2006, 2007. | |
// 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) | |
#ifndef BOOST_STATS_TRIANGULAR_HPP | |
#define BOOST_STATS_TRIANGULAR_HPP | |
// http://mathworld.wolfram.com/TriangularDistribution.html | |
// http://en.wikipedia.org/wiki/Triangular_distribution | |
#include <boost/math/distributions/fwd.hpp> | |
#include <boost/math/special_functions/expm1.hpp> | |
#include <boost/math/distributions/detail/common_error_handling.hpp> | |
#include <boost/math/distributions/complement.hpp> | |
#include <boost/math/constants/constants.hpp> | |
#include <utility> | |
namespace boost{ namespace math | |
{ | |
namespace detail | |
{ | |
template <class RealType, class Policy> | |
inline bool check_triangular_lower( | |
const char* function, | |
RealType lower, | |
RealType* result, const Policy& pol) | |
{ | |
if((boost::math::isfinite)(lower)) | |
{ // Any finite value is OK. | |
return true; | |
} | |
else | |
{ // Not finite: infinity or NaN. | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"Lower parameter is %1%, but must be finite!", lower, pol); | |
return false; | |
} | |
} // bool check_triangular_lower( | |
template <class RealType, class Policy> | |
inline bool check_triangular_mode( | |
const char* function, | |
RealType mode, | |
RealType* result, const Policy& pol) | |
{ | |
if((boost::math::isfinite)(mode)) | |
{ // any finite value is OK. | |
return true; | |
} | |
else | |
{ // Not finite: infinity or NaN. | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"Mode parameter is %1%, but must be finite!", mode, pol); | |
return false; | |
} | |
} // bool check_triangular_mode( | |
template <class RealType, class Policy> | |
inline bool check_triangular_upper( | |
const char* function, | |
RealType upper, | |
RealType* result, const Policy& pol) | |
{ | |
if((boost::math::isfinite)(upper)) | |
{ // any finite value is OK. | |
return true; | |
} | |
else | |
{ // Not finite: infinity or NaN. | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"Upper parameter is %1%, but must be finite!", upper, pol); | |
return false; | |
} | |
} // bool check_triangular_upper( | |
template <class RealType, class Policy> | |
inline bool check_triangular_x( | |
const char* function, | |
RealType const& x, | |
RealType* result, const Policy& pol) | |
{ | |
if((boost::math::isfinite)(x)) | |
{ // Any finite value is OK | |
return true; | |
} | |
else | |
{ // Not finite: infinity or NaN. | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"x parameter is %1%, but must be finite!", x, pol); | |
return false; | |
} | |
} // bool check_triangular_x | |
template <class RealType, class Policy> | |
inline bool check_triangular( | |
const char* function, | |
RealType lower, | |
RealType mode, | |
RealType upper, | |
RealType* result, const Policy& pol) | |
{ | |
if ((check_triangular_lower(function, lower, result, pol) == false) | |
|| (check_triangular_mode(function, mode, result, pol) == false) | |
|| (check_triangular_upper(function, upper, result, pol) == false)) | |
{ // Some parameter not finite. | |
return false; | |
} | |
else if (lower >= upper) // lower == upper NOT useful. | |
{ // lower >= upper. | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"lower parameter is %1%, but must be less than upper!", lower, pol); | |
return false; | |
} | |
else | |
{ // Check lower <= mode <= upper. | |
if (mode < lower) | |
{ | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"mode parameter is %1%, but must be >= than lower!", lower, pol); | |
return false; | |
} | |
if (mode > upper) | |
{ | |
*result = policies::raise_domain_error<RealType>( | |
function, | |
"mode parameter is %1%, but must be <= than upper!", upper, pol); | |
return false; | |
} | |
return true; // All OK. | |
} | |
} // bool check_triangular | |
} // namespace detail | |
template <class RealType = double, class Policy = policies::policy<> > | |
class triangular_distribution | |
{ | |
public: | |
typedef RealType value_type; | |
typedef Policy policy_type; | |
triangular_distribution(RealType lower = -1, RealType mode = 0, RealType upper = 1) | |
: m_lower(lower), m_mode(mode), m_upper(upper) // Constructor. | |
{ // Evans says 'standard triangular' is lower 0, mode 1/2, upper 1, | |
// has median sqrt(c/2) for c <=1/2 and 1 - sqrt(1-c)/2 for c >= 1/2 | |
// But this -1, 0, 1 is more useful in most applications to approximate normal distribution, | |
// where the central value is the most likely and deviations either side equally likely. | |
RealType result; | |
detail::check_triangular("boost::math::triangular_distribution<%1%>::triangular_distribution",lower, mode, upper, &result, Policy()); | |
} | |
// Accessor functions. | |
RealType lower()const | |
{ | |
return m_lower; | |
} | |
RealType mode()const | |
{ | |
return m_mode; | |
} | |
RealType upper()const | |
{ | |
return m_upper; | |
} | |
private: | |
// Data members: | |
RealType m_lower; // distribution lower aka a | |
RealType m_mode; // distribution mode aka c | |
RealType m_upper; // distribution upper aka b | |
}; // class triangular_distribution | |
typedef triangular_distribution<double> triangular; | |
template <class RealType, class Policy> | |
inline const std::pair<RealType, RealType> range(const triangular_distribution<RealType, Policy>& /* dist */) | |
{ // Range of permissible values for random variable x. | |
using boost::math::tools::max_value; | |
return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); | |
} | |
template <class RealType, class Policy> | |
inline const std::pair<RealType, RealType> support(const triangular_distribution<RealType, Policy>& dist) | |
{ // Range of supported values for random variable x. | |
// This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
return std::pair<RealType, RealType>(dist.lower(), dist.upper()); | |
} | |
template <class RealType, class Policy> | |
RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) | |
{ | |
static const char* function = "boost::math::pdf(const triangular_distribution<%1%>&, %1%)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
if(false == detail::check_triangular_x(function, x, &result, Policy())) | |
{ | |
return result; | |
} | |
if((x < lower) || (x > upper)) | |
{ | |
return 0; | |
} | |
if (x == lower) | |
{ // (mode - lower) == 0 which would lead to divide by zero! | |
return (mode == lower) ? 2 / (upper - lower) : RealType(0); | |
} | |
else if (x == upper) | |
{ | |
return (mode == upper) ? 2 / (upper - lower) : RealType(0); | |
} | |
else if (x <= mode) | |
{ | |
return 2 * (x - lower) / ((upper - lower) * (mode - lower)); | |
} | |
else | |
{ // (x > mode) | |
return 2 * (upper - x) / ((upper - lower) * (upper - mode)); | |
} | |
} // RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) | |
template <class RealType, class Policy> | |
inline RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) | |
{ | |
static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
if(false == detail::check_triangular_x(function, x, &result, Policy())) | |
{ | |
return result; | |
} | |
if((x <= lower)) | |
{ | |
return 0; | |
} | |
if (x >= upper) | |
{ | |
return 1; | |
} | |
// else lower < x < upper | |
if (x <= mode) | |
{ | |
return ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower)); | |
} | |
else | |
{ | |
return 1 - (upper - x) * (upper - x) / ((upper - lower) * (upper - mode)); | |
} | |
} // RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x) | |
template <class RealType, class Policy> | |
RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& p) | |
{ | |
BOOST_MATH_STD_USING // for ADL of std functions (sqrt). | |
static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks | |
if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
if(false == detail::check_probability(function, p, &result, Policy())) | |
{ | |
return result; | |
} | |
if(p == 0) | |
{ | |
return lower; | |
} | |
if(p == 1) | |
{ | |
return upper; | |
} | |
RealType p0 = (mode - lower) / (upper - lower); | |
RealType q = 1 - p; | |
if (p < p0) | |
{ | |
result = sqrt((upper - lower) * (mode - lower) * p) + lower; | |
} | |
else if (p == p0) | |
{ | |
result = mode; | |
} | |
else // p > p0 | |
{ | |
result = upper - sqrt((upper - lower) * (upper - mode) * q); | |
} | |
return result; | |
} // RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& q) | |
template <class RealType, class Policy> | |
RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) | |
{ | |
static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)"; | |
RealType lower = c.dist.lower(); | |
RealType mode = c.dist.mode(); | |
RealType upper = c.dist.upper(); | |
RealType x = c.param; | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
if(false == detail::check_triangular_x(function, x, &result, Policy())) | |
{ | |
return result; | |
} | |
if (x <= lower) | |
{ | |
return 1; | |
} | |
if (x >= upper) | |
{ | |
return 0; | |
} | |
if (x <= mode) | |
{ | |
return 1 - ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower)); | |
} | |
else | |
{ | |
return (upper - x) * (upper - x) / ((upper - lower) * (upper - mode)); | |
} | |
} // RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) | |
template <class RealType, class Policy> | |
RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) | |
{ | |
BOOST_MATH_STD_USING // Aid ADL for sqrt. | |
static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)"; | |
RealType l = c.dist.lower(); | |
RealType m = c.dist.mode(); | |
RealType u = c.dist.upper(); | |
RealType q = c.param; // probability 0 to 1. | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, l, m, u, &result, Policy())) | |
{ | |
return result; | |
} | |
if(false == detail::check_probability(function, q, &result, Policy())) | |
{ | |
return result; | |
} | |
if(q == 0) | |
{ | |
return u; | |
} | |
if(q == 1) | |
{ | |
return l; | |
} | |
RealType lower = c.dist.lower(); | |
RealType mode = c.dist.mode(); | |
RealType upper = c.dist.upper(); | |
RealType p = 1 - q; | |
RealType p0 = (mode - lower) / (upper - lower); | |
if(p < p0) | |
{ | |
RealType s = (upper - lower) * (mode - lower); | |
s *= p; | |
result = sqrt((upper - lower) * (mode - lower) * p) + lower; | |
} | |
else if (p == p0) | |
{ | |
result = mode; | |
} | |
else // p > p0 | |
{ | |
result = upper - sqrt((upper - lower) * (upper - mode) * q); | |
} | |
return result; | |
} // RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c) | |
template <class RealType, class Policy> | |
inline RealType mean(const triangular_distribution<RealType, Policy>& dist) | |
{ | |
static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
return (lower + upper + mode) / 3; | |
} // RealType mean(const triangular_distribution<RealType, Policy>& dist) | |
template <class RealType, class Policy> | |
inline RealType variance(const triangular_distribution<RealType, Policy>& dist) | |
{ | |
static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
return (lower * lower + upper * upper + mode * mode - lower * upper - lower * mode - upper * mode) / 18; | |
} // RealType variance(const triangular_distribution<RealType, Policy>& dist) | |
template <class RealType, class Policy> | |
inline RealType mode(const triangular_distribution<RealType, Policy>& dist) | |
{ | |
static const char* function = "boost::math::mode(const triangular_distribution<%1%>&)"; | |
RealType mode = dist.mode(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular_mode(function, mode, &result, Policy())) | |
{ // This should never happen! | |
return result; | |
} | |
return mode; | |
} // RealType mode | |
template <class RealType, class Policy> | |
inline RealType median(const triangular_distribution<RealType, Policy>& dist) | |
{ | |
BOOST_MATH_STD_USING // ADL of std functions. | |
static const char* function = "boost::math::median(const triangular_distribution<%1%>&)"; | |
RealType mode = dist.mode(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular_mode(function, mode, &result, Policy())) | |
{ // This should never happen! | |
return result; | |
} | |
RealType lower = dist.lower(); | |
RealType upper = dist.upper(); | |
if (mode < (upper - lower) / 2) | |
{ | |
return lower + sqrt((upper - lower) * (mode - lower)) / constants::root_two<RealType>(); | |
} | |
else | |
{ | |
return upper - sqrt((upper - lower) * (upper - mode)) / constants::root_two<RealType>(); | |
} | |
} // RealType mode | |
template <class RealType, class Policy> | |
inline RealType skewness(const triangular_distribution<RealType, Policy>& dist) | |
{ | |
BOOST_MATH_STD_USING // for ADL of std functions | |
using namespace boost::math::constants; // for root_two | |
static const char* function = "boost::math::skewness(const triangular_distribution<%1%>&)"; | |
RealType lower = dist.lower(); | |
RealType mode = dist.mode(); | |
RealType upper = dist.upper(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
return root_two<RealType>() * (lower + upper - 2 * mode) * (2 * lower - upper - mode) * (lower - 2 * upper + mode) / | |
(5 * pow((lower * lower + upper + upper + mode * mode - lower * upper - lower * mode - upper * mode), RealType(3)/RealType(2))); | |
} // RealType skewness(const triangular_distribution<RealType, Policy>& dist) | |
template <class RealType, class Policy> | |
inline RealType kurtosis(const triangular_distribution<RealType, Policy>& dist) | |
{ // These checks may be belt and braces as should have been checked on construction? | |
static const char* function = "boost::math::kurtosis(const triangular_distribution<%1%>&)"; | |
RealType lower = dist.lower(); | |
RealType upper = dist.upper(); | |
RealType mode = dist.mode(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
return static_cast<RealType>(12)/5; // 12/5 = 2.4; | |
} // RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist) | |
template <class RealType, class Policy> | |
inline RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist) | |
{ // These checks may be belt and braces as should have been checked on construction? | |
static const char* function = "boost::math::kurtosis_excess(const triangular_distribution<%1%>&)"; | |
RealType lower = dist.lower(); | |
RealType upper = dist.upper(); | |
RealType mode = dist.mode(); | |
RealType result; // of checks. | |
if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) | |
{ | |
return result; | |
} | |
return static_cast<RealType>(-3)/5; // - 3/5 = -0.6 | |
// Assuming mathworld really means kurtosis excess? Wikipedia now corrected to match this. | |
} | |
} // namespace math | |
} // namespace boost | |
// This include must be at the end, *after* the accessors | |
// for this distribution have been defined, in order to | |
// keep compilers that support two-phase lookup happy. | |
#include <boost/math/distributions/detail/derived_accessors.hpp> | |
#endif // BOOST_STATS_TRIANGULAR_HPP | |