diff options
Diffstat (limited to 'boost/math/distributions/skew_normal.hpp')
-rw-r--r-- | boost/math/distributions/skew_normal.hpp | 716 |
1 files changed, 716 insertions, 0 deletions
diff --git a/boost/math/distributions/skew_normal.hpp b/boost/math/distributions/skew_normal.hpp new file mode 100644 index 0000000000..58526defdb --- /dev/null +++ b/boost/math/distributions/skew_normal.hpp @@ -0,0 +1,716 @@ +// (C) Benjamin Sobotta 2012 + +// 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_SKEW_NORMAL_HPP +#define BOOST_STATS_SKEW_NORMAL_HPP + +// http://en.wikipedia.org/wiki/Skew_normal_distribution +// http://azzalini.stat.unipd.it/SN/ +// Also: +// Azzalini, A. (1985). "A class of distributions which includes the normal ones". +// Scand. J. Statist. 12: 171-178. + +#include <boost/math/distributions/fwd.hpp> // TODO add skew_normal distribution to fwd.hpp! +#include <boost/math/special_functions/owens_t.hpp> // Owen's T function +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/normal.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/tools/tuple.hpp> +#include <boost/math/tools/roots.hpp> // Newton-Raphson +#include <boost/assert.hpp> +#include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder. + +#include <utility> +#include <algorithm> // std::lower_bound, std::distance + +namespace boost{ namespace math{ + + namespace detail + { + template <class RealType, class Policy> + inline bool check_skew_normal_shape( + const char* function, + RealType shape, + RealType* result, + const Policy& pol) + { + if(!(boost::math::isfinite)(shape)) + { + *result = + policies::raise_domain_error<RealType>(function, + "Shape parameter is %1%, but must be finite!", + shape, pol); + return false; + } + return true; + } + + } // namespace detail + + template <class RealType = double, class Policy = policies::policy<> > + class skew_normal_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + skew_normal_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0) + : location_(location), scale_(scale), shape_(shape) + { // Default is a 'standard' normal distribution N01. (shape=0 results in the normal distribution with no skew) + static const char* function = "boost::math::skew_normal_distribution<%1%>::skew_normal_distribution"; + + RealType result; + detail::check_scale(function, scale, &result, Policy()); + detail::check_location(function, location, &result, Policy()); + detail::check_skew_normal_shape(function, shape, &result, Policy()); + } + + RealType location()const + { + return location_; + } + + RealType scale()const + { + return scale_; + } + + RealType shape()const + { + return shape_; + } + + + private: + // + // Data members: + // + RealType location_; // distribution location. + RealType scale_; // distribution scale. + RealType shape_; // distribution shape. + }; // class skew_normal_distribution + + typedef skew_normal_distribution<double> skew_normal; + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const skew_normal_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>()); // - to + max value. + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const skew_normal_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. + + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } + + template <class RealType, class Policy> + inline RealType pdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::pdf(const skew_normal_distribution<%1%>&, %1%)"; + if((boost::math::isinf)(x)) + { + return 0; // pdf + and - infinity is zero. + } + // Below produces MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity()) + //{ // pdf + and - infinity is zero. + // return 0; + //} + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + { + return result; + } + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = pdf(std_normal, transformed_x) * cdf(std_normal, shape*transformed_x) * 2 / scale; + + return result; + } // pdf + + template <class RealType, class Policy> + inline RealType cdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::cdf(const skew_normal_distribution<%1%>&, %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + if(x < 0) return 0; // -infinity + return 1; // + infinity + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf +infinity is unity. + // return 1; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf -infinity is zero. + // return 0; + //} + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2); + + return result; + } // cdf + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) + { + const RealType scale = c.dist.scale(); + const RealType location = c.dist.location(); + const RealType shape = c.dist.shape(); + const RealType x = c.param; + + static const char* function = "boost::math::cdf(const complement(skew_normal_distribution<%1%>&), %1%)"; + + if((boost::math::isinf)(x)) + { + if(x < 0) return 1; // cdf complement -infinity is unity. + return 0; // cdf complement +infinity is zero + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf complement +infinity is zero. + // return 0; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf complement -infinity is unity. + // return 1; + //} + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + if(false == detail::check_x(function, x, &result, Policy())) + return result; + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2); + return result; + } // cdf complement + + template <class RealType, class Policy> + inline RealType location(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.location(); + } + + template <class RealType, class Policy> + inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.scale(); + } + + template <class RealType, class Policy> + inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.shape(); + } + + template <class RealType, class Policy> + inline RealType mean(const skew_normal_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // for ADL of std functions + + using namespace boost::math::constants; + + //const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); + + //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>(); + + return dist.location() + dist.scale() * dist.shape() / sqrt(pi<RealType>()+pi<RealType>()*dist.shape()*dist.shape()) * root_two<RealType>(); + } + + template <class RealType, class Policy> + inline RealType variance(const skew_normal_distribution<RealType, Policy>& dist) + { + using namespace boost::math::constants; + + const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); + //const RealType inv_delta2 = static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape()); + + RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta2); + //RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()/inv_delta2); + + return variance; + } + + namespace detail + { + /* + TODO No closed expression for mode, so use max of pdf. + */ + + template <class RealType, class Policy> + inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist) + { // mode. + static const char* function = "mode(skew_normal_distribution<%1%> const&)"; + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + RealType result; + if(!detail::check_scale( + function, + scale, &result, Policy()) + || + !detail::check_skew_normal_shape( + function, + shape, + &result, + Policy())) + return result; + + if( shape == 0 ) + { + return location; + } + + if( shape < 0 ) + { + skew_normal_distribution<RealType, Policy> D(0, 1, -shape); + result = mode_fallback(D); + result = location-scale*result; + return result; + } + + BOOST_MATH_STD_USING + + // 21 elements + static const RealType shapes[] = { + 0.0, + 1.000000000000000e-004, + 2.069138081114790e-004, + 4.281332398719396e-004, + 8.858667904100824e-004, + 1.832980710832436e-003, + 3.792690190732250e-003, + 7.847599703514606e-003, + 1.623776739188722e-002, + 3.359818286283781e-002, + 6.951927961775606e-002, + 1.438449888287663e-001, + 2.976351441631319e-001, + 6.158482110660261e-001, + 1.274274985703135e+000, + 2.636650898730361e+000, + 5.455594781168514e+000, + 1.128837891684688e+001, + 2.335721469090121e+001, + 4.832930238571753e+001, + 1.000000000000000e+002}; + + // 21 elements + static const RealType guess[] = { + 0.0, + 5.000050000525391e-005, + 1.500015000148736e-004, + 3.500035000350010e-004, + 7.500075000752560e-004, + 1.450014500145258e-003, + 3.050030500305390e-003, + 6.250062500624765e-003, + 1.295012950129504e-002, + 2.675026750267495e-002, + 5.525055250552491e-002, + 1.132511325113255e-001, + 2.249522495224952e-001, + 3.992539925399257e-001, + 5.353553535535358e-001, + 4.954549545495457e-001, + 3.524535245352451e-001, + 2.182521825218249e-001, + 1.256512565125654e-001, + 6.945069450694508e-002, + 3.735037350373460e-002 + }; + + const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); + + typedef typename std::iterator_traits<RealType*>::difference_type diff_type; + + const diff_type d = std::distance(shapes, result_ptr); + + BOOST_ASSERT(d > static_cast<diff_type>(0)); + + // refine + if(d < static_cast<diff_type>(21)) // shape smaller 100 + { + result = guess[d-static_cast<diff_type>(1)] + + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) + * (shape-shapes[d-static_cast<diff_type>(1)]); + } + else // shape greater 100 + { + result = 1e-4; + } + + skew_normal_distribution<RealType, Policy> helper(0, 1, shape); + + result = detail::generic_find_mode_01(helper, result, function); + + result = result*scale + location; + + return result; + } // mode_fallback + + + /* + * TODO No closed expression for mode, so use f'(x) = 0 + */ + template <class RealType, class Policy> + struct skew_normal_mode_functor + { + skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist) + : distribution(dist) + { + } + + boost::math::tuple<RealType, RealType> operator()(RealType const& x) + { + normal_distribution<RealType, Policy> std_normal; + const RealType shape = distribution.shape(); + const RealType pdf_x = pdf(distribution, x); + const RealType normpdf_x = pdf(std_normal, x); + const RealType normpdf_ax = pdf(std_normal, x*shape); + RealType fx = static_cast<RealType>(2)*shape*normpdf_ax*normpdf_x - x*pdf_x; + RealType dx = static_cast<RealType>(2)*shape*x*normpdf_x*normpdf_ax*(static_cast<RealType>(1) + shape*shape) + pdf_x + x*fx; + // return both function evaluation difference f(x) and 1st derivative f'(x). + return boost::math::make_tuple(fx, -dx); + } + private: + const boost::math::skew_normal_distribution<RealType, Policy> distribution; + }; + + } // namespace detail + + template <class RealType, class Policy> + inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::mode(const skew_normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + + if( shape == 0 ) + { + return location; + } + + if( shape < 0 ) + { + skew_normal_distribution<RealType, Policy> D(0, 1, -shape); + result = mode(D); + result = location-scale*result; + return result; + } + + // 21 elements + static const RealType shapes[] = { + 0.0, + 1.000000000000000e-004, + 2.069138081114790e-004, + 4.281332398719396e-004, + 8.858667904100824e-004, + 1.832980710832436e-003, + 3.792690190732250e-003, + 7.847599703514606e-003, + 1.623776739188722e-002, + 3.359818286283781e-002, + 6.951927961775606e-002, + 1.438449888287663e-001, + 2.976351441631319e-001, + 6.158482110660261e-001, + 1.274274985703135e+000, + 2.636650898730361e+000, + 5.455594781168514e+000, + 1.128837891684688e+001, + 2.335721469090121e+001, + 4.832930238571753e+001, + 1.000000000000000e+002}; + + // 21 elements + static const RealType guess[] = { + 0.0, + 5.000050000525391e-005, + 1.500015000148736e-004, + 3.500035000350010e-004, + 7.500075000752560e-004, + 1.450014500145258e-003, + 3.050030500305390e-003, + 6.250062500624765e-003, + 1.295012950129504e-002, + 2.675026750267495e-002, + 5.525055250552491e-002, + 1.132511325113255e-001, + 2.249522495224952e-001, + 3.992539925399257e-001, + 5.353553535535358e-001, + 4.954549545495457e-001, + 3.524535245352451e-001, + 2.182521825218249e-001, + 1.256512565125654e-001, + 6.945069450694508e-002, + 3.735037350373460e-002 + }; + + const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); + + typedef typename std::iterator_traits<RealType*>::difference_type diff_type; + + const diff_type d = std::distance(shapes, result_ptr); + + BOOST_ASSERT(d > static_cast<diff_type>(0)); + + // TODO: make the search bounds smarter, depending on the shape parameter + RealType search_min = 0; // below zero was caught above + RealType search_max = 0.55; // will never go above 0.55 + + // refine + if(d < static_cast<diff_type>(21)) // shape smaller 100 + { + // it is safe to assume that d > 0, because shape==0.0 is caught earlier + result = guess[d-static_cast<diff_type>(1)] + + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) + * (shape-shapes[d-static_cast<diff_type>(1)]); + } + else // shape greater 100 + { + result = 1e-4; + search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100 + } + + const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, + boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. + + skew_normal_distribution<RealType, Policy> helper(0, 1, shape); + + result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result, + search_min, search_max, get_digits, m); + + result = result*scale + location; + + return result; + } + + + + template <class RealType, class Policy> + inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // for ADL of std functions + using namespace boost::math::constants; + + static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2); + const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); + + return factor * pow(root_two_div_pi<RealType>() * delta, 3) / + pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5)); + } + + template <class RealType, class Policy> + inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist) + { + return kurtosis_excess(dist)+static_cast<RealType>(3); + } + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist) + { + using namespace boost::math::constants; + + static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2); + + const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); + + const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta2; + const RealType y = two_div_pi<RealType>() * delta2; + + return factor * y*y / (x*x); + } + + namespace detail + { + + template <class RealType, class Policy> + struct skew_normal_quantile_functor + { + skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p) + : distribution(dist), prob(p) + { + } + + boost::math::tuple<RealType, RealType> operator()(RealType const& x) + { + RealType c = cdf(distribution, x); + RealType fx = c - prob; // Difference cdf - value - to minimize. + RealType dx = pdf(distribution, x); // pdf is 1st derivative. + // return both function evaluation difference f(x) and 1st derivative f'(x). + return boost::math::make_tuple(fx, dx); + } + private: + const boost::math::skew_normal_distribution<RealType, Policy> distribution; + RealType prob; + }; + + } // namespace detail + + template <class RealType, class Policy> + inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::quantile(const skew_normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, p, &result, Policy())) + return result; + + // compute initial guess via Cornish-Fisher expansion + RealType x = -boost::math::erfc_inv(2 * p, Policy()) * constants::root_two<RealType>(); + + // avoid unnecessary computations if there is no skew + if(shape != 0) + { + const RealType skew = skewness(dist); + const RealType exk = kurtosis_excess(dist); + + x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6) + + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24) + - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36); + } // if(shape != 0) + + result = standard_deviation(dist)*x+mean(dist); + + // handle special case of non-skew normal distribution + if(shape == 0) + return result; + + // refine the result by numerically searching the root of (p-cdf) + + const RealType search_min = range(dist).first; + const RealType search_max = range(dist).second; + + const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, + boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. + + result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(dist, p), result, + search_min, search_max, get_digits, m); + + return result; + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) + { + const RealType scale = c.dist.scale(); + const RealType location = c.dist.location(); + const RealType shape = c.dist.shape(); + + static const char* function = "boost::math::quantile(const complement(skew_normal_distribution<%1%>&), %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + RealType q = c.param; + if(false == detail::check_probability(function, q, &result, Policy())) + return result; + + skew_normal_distribution<RealType, Policy> D(-location, scale, -shape); + + result = -quantile(D, q); + + return result; + } // quantile + + +} // 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_SKEW_NORMAL_HPP + + |