| // Copyright 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 l_location = 0, RealType l_scale = 1, RealType l_shape = 0) |
| : location_(l_location), scale_(l_scale), shape_(l_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, l_scale, &result, Policy()); |
| detail::check_location(function, l_location, &result, Policy()); |
| detail::check_skew_normal_shape(function, l_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>( |
| std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(), |
| std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : 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, |
| static_cast<RealType>(1.000000000000000e-004), |
| static_cast<RealType>(2.069138081114790e-004), |
| static_cast<RealType>(4.281332398719396e-004), |
| static_cast<RealType>(8.858667904100824e-004), |
| static_cast<RealType>(1.832980710832436e-003), |
| static_cast<RealType>(3.792690190732250e-003), |
| static_cast<RealType>(7.847599703514606e-003), |
| static_cast<RealType>(1.623776739188722e-002), |
| static_cast<RealType>(3.359818286283781e-002), |
| static_cast<RealType>(6.951927961775606e-002), |
| static_cast<RealType>(1.438449888287663e-001), |
| static_cast<RealType>(2.976351441631319e-001), |
| static_cast<RealType>(6.158482110660261e-001), |
| static_cast<RealType>(1.274274985703135e+000), |
| static_cast<RealType>(2.636650898730361e+000), |
| static_cast<RealType>(5.455594781168514e+000), |
| static_cast<RealType>(1.128837891684688e+001), |
| static_cast<RealType>(2.335721469090121e+001), |
| static_cast<RealType>(4.832930238571753e+001), |
| static_cast<RealType>(1.000000000000000e+002) |
| }; |
| |
| // 21 elements |
| static const RealType guess[] = { |
| 0.0, |
| static_cast<RealType>(5.000050000525391e-005), |
| static_cast<RealType>(1.500015000148736e-004), |
| static_cast<RealType>(3.500035000350010e-004), |
| static_cast<RealType>(7.500075000752560e-004), |
| static_cast<RealType>(1.450014500145258e-003), |
| static_cast<RealType>(3.050030500305390e-003), |
| static_cast<RealType>(6.250062500624765e-003), |
| static_cast<RealType>(1.295012950129504e-002), |
| static_cast<RealType>(2.675026750267495e-002), |
| static_cast<RealType>(5.525055250552491e-002), |
| static_cast<RealType>(1.132511325113255e-001), |
| static_cast<RealType>(2.249522495224952e-001), |
| static_cast<RealType>(3.992539925399257e-001), |
| static_cast<RealType>(5.353553535535358e-001), |
| static_cast<RealType>(4.954549545495457e-001), |
| static_cast<RealType>(3.524535245352451e-001), |
| static_cast<RealType>(2.182521825218249e-001), |
| static_cast<RealType>(1.256512565125654e-001), |
| static_cast<RealType>(6.945069450694508e-002), |
| static_cast<RealType>(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.55f; // 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-4f; |
| 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 |
| |
| |
