761 lines · plain
1// Copyright Benjamin Sobotta 20122 3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_STATS_SKEW_NORMAL_HPP8#define BOOST_STATS_SKEW_NORMAL_HPP9 10// http://en.wikipedia.org/wiki/Skew_normal_distribution11// http://azzalini.stat.unipd.it/SN/12// Also:13// Azzalini, A. (1985). "A class of distributions which includes the normal ones".14// Scand. J. Statist. 12: 171-178.15 16#include <boost/math/distributions/fwd.hpp>17#include <boost/math/special_functions/owens_t.hpp> // Owen's T function18#include <boost/math/distributions/complement.hpp>19#include <boost/math/distributions/normal.hpp>20#include <boost/math/distributions/detail/common_error_handling.hpp>21#include <boost/math/constants/constants.hpp>22#include <boost/math/tools/tuple.hpp>23#include <boost/math/tools/roots.hpp> // Newton-Raphson24#include <boost/math/tools/assert.hpp>25#include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder.26 27#include <utility>28#include <algorithm> // std::lower_bound, std::distance29 30#ifdef BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS31extern std::uintmax_t global_iter_count;32#endif33 34namespace boost{ namespace math{35 36 namespace detail37 {38 template <class RealType, class Policy>39 inline bool check_skew_normal_shape(40 const char* function,41 RealType shape,42 RealType* result,43 const Policy& pol)44 {45 if(!(boost::math::isfinite)(shape))46 {47 *result =48 policies::raise_domain_error<RealType>(function,49 "Shape parameter is %1%, but must be finite!",50 shape, pol);51 return false;52 }53 return true;54 }55 56 } // namespace detail57 58 template <class RealType = double, class Policy = policies::policy<> >59 class skew_normal_distribution60 {61 public:62 typedef RealType value_type;63 typedef Policy policy_type;64 65 skew_normal_distribution(RealType l_location = 0, RealType l_scale = 1, RealType l_shape = 0)66 : location_(l_location), scale_(l_scale), shape_(l_shape)67 { // Default is a 'standard' normal distribution N01. (shape=0 results in the normal distribution with no skew)68 static const char* function = "boost::math::skew_normal_distribution<%1%>::skew_normal_distribution";69 70 RealType result;71 detail::check_scale(function, l_scale, &result, Policy());72 detail::check_location(function, l_location, &result, Policy());73 detail::check_skew_normal_shape(function, l_shape, &result, Policy());74 }75 76 RealType location()const77 {78 return location_;79 }80 81 RealType scale()const82 {83 return scale_;84 }85 86 RealType shape()const87 {88 return shape_;89 }90 91 92 private:93 //94 // Data members:95 //96 RealType location_; // distribution location.97 RealType scale_; // distribution scale.98 RealType shape_; // distribution shape.99 }; // class skew_normal_distribution100 101 typedef skew_normal_distribution<double> skew_normal;102 103 #ifdef __cpp_deduction_guides104 template <class RealType>105 skew_normal_distribution(RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;106 template <class RealType>107 skew_normal_distribution(RealType,RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;108 template <class RealType>109 skew_normal_distribution(RealType,RealType,RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;110 #endif111 112 template <class RealType, class Policy>113 inline const std::pair<RealType, RealType> range(const skew_normal_distribution<RealType, Policy>& /*dist*/)114 { // Range of permissible values for random variable x.115 using boost::math::tools::max_value;116 return std::pair<RealType, RealType>(117 std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(),118 std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); // - to + max value.119 }120 121 template <class RealType, class Policy>122 inline const std::pair<RealType, RealType> support(const skew_normal_distribution<RealType, Policy>& /*dist*/)123 { // Range of supported values for random variable x.124 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.125 126 using boost::math::tools::max_value;127 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value.128 }129 130 template <class RealType, class Policy>131 inline RealType pdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x)132 {133 const RealType scale = dist.scale();134 const RealType location = dist.location();135 const RealType shape = dist.shape();136 137 static const char* function = "boost::math::pdf(const skew_normal_distribution<%1%>&, %1%)";138 139 RealType result = 0;140 if(false == detail::check_scale(function, scale, &result, Policy()))141 {142 return result;143 }144 if(false == detail::check_location(function, location, &result, Policy()))145 {146 return result;147 }148 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))149 {150 return result;151 }152 if((boost::math::isinf)(x))153 {154 return 0; // pdf + and - infinity is zero.155 }156 // Below produces MSVC 4127 warnings, so the above used instead.157 //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())158 //{ // pdf + and - infinity is zero.159 // return 0;160 //}161 if(false == detail::check_x(function, x, &result, Policy()))162 {163 return result;164 }165 166 const RealType transformed_x = (x-location)/scale;167 168 normal_distribution<RealType, Policy> std_normal;169 170 result = pdf(std_normal, transformed_x) * cdf(std_normal, shape*transformed_x) * 2 / scale;171 172 return result;173 } // pdf174 175 template <class RealType, class Policy>176 inline RealType cdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x)177 {178 const RealType scale = dist.scale();179 const RealType location = dist.location();180 const RealType shape = dist.shape();181 182 static const char* function = "boost::math::cdf(const skew_normal_distribution<%1%>&, %1%)";183 RealType result = 0;184 if(false == detail::check_scale(function, scale, &result, Policy()))185 {186 return result;187 }188 if(false == detail::check_location(function, location, &result, Policy()))189 {190 return result;191 }192 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))193 {194 return result;195 }196 if((boost::math::isinf)(x))197 {198 if(x < 0) return 0; // -infinity199 return 1; // + infinity200 }201 // These produce MSVC 4127 warnings, so the above used instead.202 //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())203 //{ // cdf +infinity is unity.204 // return 1;205 //}206 //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())207 //{ // cdf -infinity is zero.208 // return 0;209 //}210 if(false == detail::check_x(function, x, &result, Policy()))211 {212 return result;213 }214 215 const RealType transformed_x = (x-location)/scale;216 217 normal_distribution<RealType, Policy> std_normal;218 219 result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2);220 221 return result;222 } // cdf223 224 template <class RealType, class Policy>225 inline RealType cdf(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c)226 {227 const RealType scale = c.dist.scale();228 const RealType location = c.dist.location();229 const RealType shape = c.dist.shape();230 const RealType x = c.param;231 232 static const char* function = "boost::math::cdf(const complement(skew_normal_distribution<%1%>&), %1%)";233 234 if((boost::math::isinf)(x))235 {236 if(x < 0) return 1; // cdf complement -infinity is unity.237 return 0; // cdf complement +infinity is zero238 }239 // These produce MSVC 4127 warnings, so the above used instead.240 //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())241 //{ // cdf complement +infinity is zero.242 // return 0;243 //}244 //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())245 //{ // cdf complement -infinity is unity.246 // return 1;247 //}248 RealType result = 0;249 if(false == detail::check_scale(function, scale, &result, Policy()))250 return result;251 if(false == detail::check_location(function, location, &result, Policy()))252 return result;253 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))254 return result;255 if(false == detail::check_x(function, x, &result, Policy()))256 return result;257 258 const RealType transformed_x = (x-location)/scale;259 260 normal_distribution<RealType, Policy> std_normal;261 262 result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2);263 return result;264 } // cdf complement265 266 template <class RealType, class Policy>267 inline RealType location(const skew_normal_distribution<RealType, Policy>& dist)268 {269 return dist.location();270 }271 272 template <class RealType, class Policy>273 inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist)274 {275 return dist.scale();276 }277 278 template <class RealType, class Policy>279 inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist)280 {281 return dist.shape();282 }283 284 template <class RealType, class Policy>285 inline RealType mean(const skew_normal_distribution<RealType, Policy>& dist)286 {287 BOOST_MATH_STD_USING // for ADL of std functions288 289 using namespace boost::math::constants;290 291 //const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());292 293 //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>();294 295 return dist.location() + dist.scale() * dist.shape() / sqrt(pi<RealType>()+pi<RealType>()*dist.shape()*dist.shape()) * root_two<RealType>();296 }297 298 template <class RealType, class Policy>299 inline RealType variance(const skew_normal_distribution<RealType, Policy>& dist)300 {301 using namespace boost::math::constants;302 303 const RealType delta2 = dist.shape() != 0 ? static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())) : static_cast<RealType>(0);304 //const RealType inv_delta2 = static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape());305 306 RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta2);307 //RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()/inv_delta2);308 309 return variance;310 }311 312 namespace detail313 {314 /*315 TODO No closed expression for mode, so use max of pdf.316 */317 318 template <class RealType, class Policy>319 inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist)320 { // mode.321 static const char* function = "mode(skew_normal_distribution<%1%> const&)";322 const RealType scale = dist.scale();323 const RealType location = dist.location();324 const RealType shape = dist.shape();325 326 RealType result;327 if(!detail::check_scale(328 function,329 scale, &result, Policy())330 ||331 !detail::check_skew_normal_shape(332 function,333 shape,334 &result,335 Policy()))336 return result;337 338 if( shape == 0 )339 {340 return location;341 }342 343 if( shape < 0 )344 {345 skew_normal_distribution<RealType, Policy> D(0, 1, -shape);346 result = mode_fallback(D);347 result = location-scale*result;348 return result;349 }350 351 BOOST_MATH_STD_USING352 353 // 21 elements354 static const RealType shapes[] = {355 0.0,356 1.000000000000000e-004,357 2.069138081114790e-004,358 4.281332398719396e-004,359 8.858667904100824e-004,360 1.832980710832436e-003,361 3.792690190732250e-003,362 7.847599703514606e-003,363 1.623776739188722e-002,364 3.359818286283781e-002,365 6.951927961775606e-002,366 1.438449888287663e-001,367 2.976351441631319e-001,368 6.158482110660261e-001,369 1.274274985703135e+000,370 2.636650898730361e+000,371 5.455594781168514e+000,372 1.128837891684688e+001,373 2.335721469090121e+001,374 4.832930238571753e+001,375 1.000000000000000e+002};376 377 // 21 elements378 static const RealType guess[] = {379 0.0,380 5.000050000525391e-005,381 1.500015000148736e-004,382 3.500035000350010e-004,383 7.500075000752560e-004,384 1.450014500145258e-003,385 3.050030500305390e-003,386 6.250062500624765e-003,387 1.295012950129504e-002,388 2.675026750267495e-002,389 5.525055250552491e-002,390 1.132511325113255e-001,391 2.249522495224952e-001,392 3.992539925399257e-001,393 5.353553535535358e-001,394 4.954549545495457e-001,395 3.524535245352451e-001,396 2.182521825218249e-001,397 1.256512565125654e-001,398 6.945069450694508e-002,399 3.735037350373460e-002400 };401 402 const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);403 404 typedef typename std::iterator_traits<RealType*>::difference_type diff_type;405 406 const diff_type d = std::distance(shapes, result_ptr);407 408 BOOST_MATH_ASSERT(d > static_cast<diff_type>(0));409 410 // refine411 if(d < static_cast<diff_type>(21)) // shape smaller 100412 {413 result = guess[d-static_cast<diff_type>(1)]414 + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])415 * (shape-shapes[d-static_cast<diff_type>(1)]);416 }417 else // shape greater 100418 {419 result = 1e-4;420 }421 422 skew_normal_distribution<RealType, Policy> helper(0, 1, shape);423 424 result = detail::generic_find_mode_01(helper, result, function);425 426 result = result*scale + location;427 428 return result;429 } // mode_fallback430 431 432 /*433 * TODO No closed expression for mode, so use f'(x) = 0434 */435 template <class RealType, class Policy>436 struct skew_normal_mode_functor437 {438 skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist)439 : distribution(dist)440 {441 }442 443 boost::math::tuple<RealType, RealType> operator()(RealType const& x)444 {445 normal_distribution<RealType, Policy> std_normal;446 const RealType shape = distribution.shape();447 const RealType pdf_x = pdf(distribution, x);448 const RealType normpdf_x = pdf(std_normal, x);449 const RealType normpdf_ax = pdf(std_normal, x*shape);450 RealType fx = static_cast<RealType>(2)*shape*normpdf_ax*normpdf_x - x*pdf_x;451 RealType dx = static_cast<RealType>(2)*shape*x*normpdf_x*normpdf_ax*(static_cast<RealType>(1) + shape*shape) + pdf_x + x*fx;452 // return both function evaluation difference f(x) and 1st derivative f'(x).453 return boost::math::make_tuple(fx, -dx);454 }455 private:456 const boost::math::skew_normal_distribution<RealType, Policy> distribution;457 };458 459 } // namespace detail460 461 template <class RealType, class Policy>462 inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist)463 {464 const RealType scale = dist.scale();465 const RealType location = dist.location();466 const RealType shape = dist.shape();467 468 static const char* function = "boost::math::mode(const skew_normal_distribution<%1%>&, %1%)";469 470 RealType result = 0;471 if(false == detail::check_scale(function, scale, &result, Policy()))472 return result;473 if(false == detail::check_location(function, location, &result, Policy()))474 return result;475 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))476 return result;477 478 if( shape == 0 )479 {480 return location;481 }482 483 if( shape < 0 )484 {485 skew_normal_distribution<RealType, Policy> D(0, 1, -shape);486 result = mode(D);487 result = location-scale*result;488 return result;489 }490 491 // 21 elements492 static const RealType shapes[] = {493 static_cast<RealType>(0.0),494 static_cast<RealType>(1.000000000000000e-004),495 static_cast<RealType>(2.069138081114790e-004),496 static_cast<RealType>(4.281332398719396e-004),497 static_cast<RealType>(8.858667904100824e-004),498 static_cast<RealType>(1.832980710832436e-003),499 static_cast<RealType>(3.792690190732250e-003),500 static_cast<RealType>(7.847599703514606e-003),501 static_cast<RealType>(1.623776739188722e-002),502 static_cast<RealType>(3.359818286283781e-002),503 static_cast<RealType>(6.951927961775606e-002),504 static_cast<RealType>(1.438449888287663e-001),505 static_cast<RealType>(2.976351441631319e-001),506 static_cast<RealType>(6.158482110660261e-001),507 static_cast<RealType>(1.274274985703135e+000),508 static_cast<RealType>(2.636650898730361e+000),509 static_cast<RealType>(5.455594781168514e+000),510 static_cast<RealType>(1.128837891684688e+001),511 static_cast<RealType>(2.335721469090121e+001),512 static_cast<RealType>(4.832930238571753e+001),513 static_cast<RealType>(1.000000000000000e+002)514 };515 516 // 21 elements517 static const RealType guess[] = {518 static_cast<RealType>(0.0),519 static_cast<RealType>(5.000050000525391e-005),520 static_cast<RealType>(1.500015000148736e-004),521 static_cast<RealType>(3.500035000350010e-004),522 static_cast<RealType>(7.500075000752560e-004),523 static_cast<RealType>(1.450014500145258e-003),524 static_cast<RealType>(3.050030500305390e-003),525 static_cast<RealType>(6.250062500624765e-003),526 static_cast<RealType>(1.295012950129504e-002),527 static_cast<RealType>(2.675026750267495e-002),528 static_cast<RealType>(5.525055250552491e-002),529 static_cast<RealType>(1.132511325113255e-001),530 static_cast<RealType>(2.249522495224952e-001),531 static_cast<RealType>(3.992539925399257e-001),532 static_cast<RealType>(5.353553535535358e-001),533 static_cast<RealType>(4.954549545495457e-001),534 static_cast<RealType>(3.524535245352451e-001),535 static_cast<RealType>(2.182521825218249e-001),536 static_cast<RealType>(1.256512565125654e-001),537 static_cast<RealType>(6.945069450694508e-002),538 static_cast<RealType>(3.735037350373460e-002)539 };540 541 const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);542 543 typedef typename std::iterator_traits<RealType*>::difference_type diff_type;544 545 const diff_type d = std::distance(shapes, result_ptr);546 547 BOOST_MATH_ASSERT(d > static_cast<diff_type>(0));548 549 // TODO: make the search bounds smarter, depending on the shape parameter550 RealType search_min = 0; // below zero was caught above551 RealType search_max = 0.55f; // will never go above 0.55552 553 // refine554 if(d < static_cast<diff_type>(21)) // shape smaller 100555 {556 // it is safe to assume that d > 0, because shape==0.0 is caught earlier557 result = guess[d-static_cast<diff_type>(1)]558 + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])559 * (shape-shapes[d-static_cast<diff_type>(1)]);560 }561 else // shape greater 100562 {563 result = 1e-4f;564 search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100565 }566 567 const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,568 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.569 570 skew_normal_distribution<RealType, Policy> helper(0, 1, shape);571 572 result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result,573 search_min, search_max, get_digits, max_iter);574 if (max_iter >= policies::get_max_root_iterations<Policy>())575 {576 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE577 " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE578 }579 580 result = result*scale + location;581 582 return result;583 }584 585 586 587 template <class RealType, class Policy>588 inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist)589 {590 BOOST_MATH_STD_USING // for ADL of std functions591 using namespace boost::math::constants;592 593 static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2);594 const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());595 596 return static_cast<RealType>(factor * pow(root_two_div_pi<RealType>() * delta, 3) /597 pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5)));598 }599 600 template <class RealType, class Policy>601 inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist)602 {603 return kurtosis_excess(dist)+static_cast<RealType>(3);604 }605 606 template <class RealType, class Policy>607 inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist)608 {609 using namespace boost::math::constants;610 611 static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2);612 613 const RealType delta2 = dist.shape() != 0 ? static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())) : static_cast<RealType>(0);614 615 const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta2;616 const RealType y = two_div_pi<RealType>() * delta2;617 618 return factor * y*y / (x*x);619 }620 621 template <class RealType, class Policy>622 inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p)623 {624 const RealType scale = dist.scale();625 const RealType location = dist.location();626 const RealType shape = dist.shape();627 628 static const char* function = "boost::math::quantile(const skew_normal_distribution<%1%>&, %1%)";629 630 RealType result = 0;631 if(false == detail::check_scale(function, scale, &result, Policy()))632 return result;633 if(false == detail::check_location(function, location, &result, Policy()))634 return result;635 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))636 return result;637 if(false == detail::check_probability(function, p, &result, Policy()))638 return result;639 640 // Compute initial guess via Cornish-Fisher expansion.641 RealType x = -boost::math::erfc_inv(2 * p, Policy()) * constants::root_two<RealType>();642 643 // Avoid unnecessary computations if there is no skew.644 if(shape != 0)645 {646 const RealType skew = skewness(dist);647 const RealType exk = kurtosis_excess(dist);648 649 x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6)650 + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24)651 - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36);652 } // if(shape != 0)653 654 result = standard_deviation(dist)*x+mean(dist);655 656 // handle special case of non-skew normal distribution.657 if(shape == 0)658 return result;659 660 // refine the result by numerically searching the root of (p-cdf)661 662 const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,663 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.664 665 if (result == 0)666 result = tools::min_value<RealType>(); // we need to be one side of zero or the other for the root finder to work.667 668 auto fun = [&, dist, p](const RealType& x)->RealType { return cdf(dist, x) - p; };669 670 RealType f_result = fun(result);671 672 if (f_result == 0)673 return result;674 675 if (f_result * result > 0)676 {677 // If the root is in the direction of zero, we need to check that we're the correct side of it:678 RealType f_zero = fun(static_cast<RealType>(0));679 if (f_zero * f_result > 0)680 {681 // we're the wrong side of zero:682 result = -result;683 f_result = fun(result);684 }685 }686 687 RealType scaling_factor = 1.25;688 if (f_result * result > 0)689 {690 // We're heading towards zero... it's a long way down so use a larger scaling factor:691 scaling_factor = 16;692 }693 694 auto p_result = tools::bracket_and_solve_root(fun, result, scaling_factor, true, tools::eps_tolerance<RealType>(get_digits), max_iter, Policy());695 696#ifdef BOOST_MATH_INSTRUMENT_SKEW_NORMAL_ITERATIONS697 global_iter_count += max_iter;698#endif699 700 result = (p_result.first + p_result.second) / 2;701 702 //703 // Try one last Newton step, just to close up the interval:704 //705 RealType step = fun(result) / pdf(dist, result);706 707 if (result - step <= p_result.first)708 result = p_result.first;709 else if (result - step >= p_result.second)710 result = p_result.second;711 else712 result -= step;713 714 if (max_iter >= policies::get_max_root_iterations<Policy>())715 {716 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time: either there is no answer to quantile" // LCOV_EXCL_LINE717 " or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE718 }719 720 return result;721 } // quantile722 723 template <class RealType, class Policy>724 inline RealType quantile(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c)725 {726 const RealType scale = c.dist.scale();727 const RealType location = c.dist.location();728 const RealType shape = c.dist.shape();729 730 static const char* function = "boost::math::quantile(const complement(skew_normal_distribution<%1%>&), %1%)";731 RealType result = 0;732 if(false == detail::check_scale(function, scale, &result, Policy()))733 return result;734 if(false == detail::check_location(function, location, &result, Policy()))735 return result;736 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))737 return result;738 RealType q = c.param;739 if(false == detail::check_probability(function, q, &result, Policy()))740 return result;741 742 skew_normal_distribution<RealType, Policy> D(-location, scale, -shape);743 744 result = -quantile(D, q);745 746 return result;747 } // quantile748 749 750} // namespace math751} // namespace boost752 753// This include must be at the end, *after* the accessors754// for this distribution have been defined, in order to755// keep compilers that support two-phase lookup happy.756#include <boost/math/distributions/detail/derived_accessors.hpp>757 758#endif // BOOST_STATS_SKEW_NORMAL_HPP759 760 761