1 // Copyright Benjamin Sobotta 2012 2 3 // Use, modification and distribution are subject to the 4 // Boost Software License, Version 1.0. (See accompanying file 5 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_STATS_SKEW_NORMAL_HPP 8 #define BOOST_STATS_SKEW_NORMAL_HPP 9 10 // http://en.wikipedia.org/wiki/Skew_normal_distribution 11 // 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> // TODO add skew_normal distribution to fwd.hpp! 17 #include <boost/math/special_functions/owens_t.hpp> // Owen's T function 18 #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-Raphson 24 #include <boost/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::distance 29 30 namespace boost{ namespace math{ 31 32 namespace detail 33 { 34 template <class RealType, class Policy> check_skew_normal_shape(const char * function,RealType shape,RealType * result,const Policy & pol)35 inline bool check_skew_normal_shape( 36 const char* function, 37 RealType shape, 38 RealType* result, 39 const Policy& pol) 40 { 41 if(!(boost::math::isfinite)(shape)) 42 { 43 *result = 44 policies::raise_domain_error<RealType>(function, 45 "Shape parameter is %1%, but must be finite!", 46 shape, pol); 47 return false; 48 } 49 return true; 50 } 51 52 } // namespace detail 53 54 template <class RealType = double, class Policy = policies::policy<> > 55 class skew_normal_distribution 56 { 57 public: 58 typedef RealType value_type; 59 typedef Policy policy_type; 60 skew_normal_distribution(RealType l_location=0,RealType l_scale=1,RealType l_shape=0)61 skew_normal_distribution(RealType l_location = 0, RealType l_scale = 1, RealType l_shape = 0) 62 : location_(l_location), scale_(l_scale), shape_(l_shape) 63 { // Default is a 'standard' normal distribution N01. (shape=0 results in the normal distribution with no skew) 64 static const char* function = "boost::math::skew_normal_distribution<%1%>::skew_normal_distribution"; 65 66 RealType result; 67 detail::check_scale(function, l_scale, &result, Policy()); 68 detail::check_location(function, l_location, &result, Policy()); 69 detail::check_skew_normal_shape(function, l_shape, &result, Policy()); 70 } 71 location() const72 RealType location()const 73 { 74 return location_; 75 } 76 scale() const77 RealType scale()const 78 { 79 return scale_; 80 } 81 shape() const82 RealType shape()const 83 { 84 return shape_; 85 } 86 87 88 private: 89 // 90 // Data members: 91 // 92 RealType location_; // distribution location. 93 RealType scale_; // distribution scale. 94 RealType shape_; // distribution shape. 95 }; // class skew_normal_distribution 96 97 typedef skew_normal_distribution<double> skew_normal; 98 99 template <class RealType, class Policy> range(const skew_normal_distribution<RealType,Policy> &)100 inline const std::pair<RealType, RealType> range(const skew_normal_distribution<RealType, Policy>& /*dist*/) 101 { // Range of permissible values for random variable x. 102 using boost::math::tools::max_value; 103 return std::pair<RealType, RealType>( 104 std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(), 105 std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); // - to + max value. 106 } 107 108 template <class RealType, class Policy> support(const skew_normal_distribution<RealType,Policy> &)109 inline const std::pair<RealType, RealType> support(const skew_normal_distribution<RealType, Policy>& /*dist*/) 110 { // Range of supported values for random variable x. 111 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 112 113 using boost::math::tools::max_value; 114 return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. 115 } 116 117 template <class RealType, class Policy> pdf(const skew_normal_distribution<RealType,Policy> & dist,const RealType & x)118 inline RealType pdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) 119 { 120 const RealType scale = dist.scale(); 121 const RealType location = dist.location(); 122 const RealType shape = dist.shape(); 123 124 static const char* function = "boost::math::pdf(const skew_normal_distribution<%1%>&, %1%)"; 125 if((boost::math::isinf)(x)) 126 { 127 return 0; // pdf + and - infinity is zero. 128 } 129 // Below produces MSVC 4127 warnings, so the above used instead. 130 //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity()) 131 //{ // pdf + and - infinity is zero. 132 // return 0; 133 //} 134 135 RealType result = 0; 136 if(false == detail::check_scale(function, scale, &result, Policy())) 137 { 138 return result; 139 } 140 if(false == detail::check_location(function, location, &result, Policy())) 141 { 142 return result; 143 } 144 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 145 { 146 return result; 147 } 148 if(false == detail::check_x(function, x, &result, Policy())) 149 { 150 return result; 151 } 152 153 const RealType transformed_x = (x-location)/scale; 154 155 normal_distribution<RealType, Policy> std_normal; 156 157 result = pdf(std_normal, transformed_x) * cdf(std_normal, shape*transformed_x) * 2 / scale; 158 159 return result; 160 } // pdf 161 162 template <class RealType, class Policy> cdf(const skew_normal_distribution<RealType,Policy> & dist,const RealType & x)163 inline RealType cdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) 164 { 165 const RealType scale = dist.scale(); 166 const RealType location = dist.location(); 167 const RealType shape = dist.shape(); 168 169 static const char* function = "boost::math::cdf(const skew_normal_distribution<%1%>&, %1%)"; 170 RealType result = 0; 171 if(false == detail::check_scale(function, scale, &result, Policy())) 172 { 173 return result; 174 } 175 if(false == detail::check_location(function, location, &result, Policy())) 176 { 177 return result; 178 } 179 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 180 { 181 return result; 182 } 183 if((boost::math::isinf)(x)) 184 { 185 if(x < 0) return 0; // -infinity 186 return 1; // + infinity 187 } 188 // These produce MSVC 4127 warnings, so the above used instead. 189 //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) 190 //{ // cdf +infinity is unity. 191 // return 1; 192 //} 193 //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) 194 //{ // cdf -infinity is zero. 195 // return 0; 196 //} 197 if(false == detail::check_x(function, x, &result, Policy())) 198 { 199 return result; 200 } 201 202 const RealType transformed_x = (x-location)/scale; 203 204 normal_distribution<RealType, Policy> std_normal; 205 206 result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2); 207 208 return result; 209 } // cdf 210 211 template <class RealType, class Policy> cdf(const complemented2_type<skew_normal_distribution<RealType,Policy>,RealType> & c)212 inline RealType cdf(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) 213 { 214 const RealType scale = c.dist.scale(); 215 const RealType location = c.dist.location(); 216 const RealType shape = c.dist.shape(); 217 const RealType x = c.param; 218 219 static const char* function = "boost::math::cdf(const complement(skew_normal_distribution<%1%>&), %1%)"; 220 221 if((boost::math::isinf)(x)) 222 { 223 if(x < 0) return 1; // cdf complement -infinity is unity. 224 return 0; // cdf complement +infinity is zero 225 } 226 // These produce MSVC 4127 warnings, so the above used instead. 227 //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) 228 //{ // cdf complement +infinity is zero. 229 // return 0; 230 //} 231 //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) 232 //{ // cdf complement -infinity is unity. 233 // return 1; 234 //} 235 RealType result = 0; 236 if(false == detail::check_scale(function, scale, &result, Policy())) 237 return result; 238 if(false == detail::check_location(function, location, &result, Policy())) 239 return result; 240 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 241 return result; 242 if(false == detail::check_x(function, x, &result, Policy())) 243 return result; 244 245 const RealType transformed_x = (x-location)/scale; 246 247 normal_distribution<RealType, Policy> std_normal; 248 249 result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2); 250 return result; 251 } // cdf complement 252 253 template <class RealType, class Policy> location(const skew_normal_distribution<RealType,Policy> & dist)254 inline RealType location(const skew_normal_distribution<RealType, Policy>& dist) 255 { 256 return dist.location(); 257 } 258 259 template <class RealType, class Policy> scale(const skew_normal_distribution<RealType,Policy> & dist)260 inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist) 261 { 262 return dist.scale(); 263 } 264 265 template <class RealType, class Policy> shape(const skew_normal_distribution<RealType,Policy> & dist)266 inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist) 267 { 268 return dist.shape(); 269 } 270 271 template <class RealType, class Policy> mean(const skew_normal_distribution<RealType,Policy> & dist)272 inline RealType mean(const skew_normal_distribution<RealType, Policy>& dist) 273 { 274 BOOST_MATH_STD_USING // for ADL of std functions 275 276 using namespace boost::math::constants; 277 278 //const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); 279 280 //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>(); 281 282 return dist.location() + dist.scale() * dist.shape() / sqrt(pi<RealType>()+pi<RealType>()*dist.shape()*dist.shape()) * root_two<RealType>(); 283 } 284 285 template <class RealType, class Policy> variance(const skew_normal_distribution<RealType,Policy> & dist)286 inline RealType variance(const skew_normal_distribution<RealType, Policy>& dist) 287 { 288 using namespace boost::math::constants; 289 290 const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); 291 //const RealType inv_delta2 = static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape()); 292 293 RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta2); 294 //RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()/inv_delta2); 295 296 return variance; 297 } 298 299 namespace detail 300 { 301 /* 302 TODO No closed expression for mode, so use max of pdf. 303 */ 304 305 template <class RealType, class Policy> mode_fallback(const skew_normal_distribution<RealType,Policy> & dist)306 inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist) 307 { // mode. 308 static const char* function = "mode(skew_normal_distribution<%1%> const&)"; 309 const RealType scale = dist.scale(); 310 const RealType location = dist.location(); 311 const RealType shape = dist.shape(); 312 313 RealType result; 314 if(!detail::check_scale( 315 function, 316 scale, &result, Policy()) 317 || 318 !detail::check_skew_normal_shape( 319 function, 320 shape, 321 &result, 322 Policy())) 323 return result; 324 325 if( shape == 0 ) 326 { 327 return location; 328 } 329 330 if( shape < 0 ) 331 { 332 skew_normal_distribution<RealType, Policy> D(0, 1, -shape); 333 result = mode_fallback(D); 334 result = location-scale*result; 335 return result; 336 } 337 338 BOOST_MATH_STD_USING 339 340 // 21 elements 341 static const RealType shapes[] = { 342 0.0, 343 1.000000000000000e-004, 344 2.069138081114790e-004, 345 4.281332398719396e-004, 346 8.858667904100824e-004, 347 1.832980710832436e-003, 348 3.792690190732250e-003, 349 7.847599703514606e-003, 350 1.623776739188722e-002, 351 3.359818286283781e-002, 352 6.951927961775606e-002, 353 1.438449888287663e-001, 354 2.976351441631319e-001, 355 6.158482110660261e-001, 356 1.274274985703135e+000, 357 2.636650898730361e+000, 358 5.455594781168514e+000, 359 1.128837891684688e+001, 360 2.335721469090121e+001, 361 4.832930238571753e+001, 362 1.000000000000000e+002}; 363 364 // 21 elements 365 static const RealType guess[] = { 366 0.0, 367 5.000050000525391e-005, 368 1.500015000148736e-004, 369 3.500035000350010e-004, 370 7.500075000752560e-004, 371 1.450014500145258e-003, 372 3.050030500305390e-003, 373 6.250062500624765e-003, 374 1.295012950129504e-002, 375 2.675026750267495e-002, 376 5.525055250552491e-002, 377 1.132511325113255e-001, 378 2.249522495224952e-001, 379 3.992539925399257e-001, 380 5.353553535535358e-001, 381 4.954549545495457e-001, 382 3.524535245352451e-001, 383 2.182521825218249e-001, 384 1.256512565125654e-001, 385 6.945069450694508e-002, 386 3.735037350373460e-002 387 }; 388 389 const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); 390 391 typedef typename std::iterator_traits<RealType*>::difference_type diff_type; 392 393 const diff_type d = std::distance(shapes, result_ptr); 394 395 BOOST_ASSERT(d > static_cast<diff_type>(0)); 396 397 // refine 398 if(d < static_cast<diff_type>(21)) // shape smaller 100 399 { 400 result = guess[d-static_cast<diff_type>(1)] 401 + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) 402 * (shape-shapes[d-static_cast<diff_type>(1)]); 403 } 404 else // shape greater 100 405 { 406 result = 1e-4; 407 } 408 409 skew_normal_distribution<RealType, Policy> helper(0, 1, shape); 410 411 result = detail::generic_find_mode_01(helper, result, function); 412 413 result = result*scale + location; 414 415 return result; 416 } // mode_fallback 417 418 419 /* 420 * TODO No closed expression for mode, so use f'(x) = 0 421 */ 422 template <class RealType, class Policy> 423 struct skew_normal_mode_functor 424 { skew_normal_mode_functorboost::math::detail::skew_normal_mode_functor425 skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist) 426 : distribution(dist) 427 { 428 } 429 operator ()boost::math::detail::skew_normal_mode_functor430 boost::math::tuple<RealType, RealType> operator()(RealType const& x) 431 { 432 normal_distribution<RealType, Policy> std_normal; 433 const RealType shape = distribution.shape(); 434 const RealType pdf_x = pdf(distribution, x); 435 const RealType normpdf_x = pdf(std_normal, x); 436 const RealType normpdf_ax = pdf(std_normal, x*shape); 437 RealType fx = static_cast<RealType>(2)*shape*normpdf_ax*normpdf_x - x*pdf_x; 438 RealType dx = static_cast<RealType>(2)*shape*x*normpdf_x*normpdf_ax*(static_cast<RealType>(1) + shape*shape) + pdf_x + x*fx; 439 // return both function evaluation difference f(x) and 1st derivative f'(x). 440 return boost::math::make_tuple(fx, -dx); 441 } 442 private: 443 const boost::math::skew_normal_distribution<RealType, Policy> distribution; 444 }; 445 446 } // namespace detail 447 448 template <class RealType, class Policy> mode(const skew_normal_distribution<RealType,Policy> & dist)449 inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist) 450 { 451 const RealType scale = dist.scale(); 452 const RealType location = dist.location(); 453 const RealType shape = dist.shape(); 454 455 static const char* function = "boost::math::mode(const skew_normal_distribution<%1%>&, %1%)"; 456 457 RealType result = 0; 458 if(false == detail::check_scale(function, scale, &result, Policy())) 459 return result; 460 if(false == detail::check_location(function, location, &result, Policy())) 461 return result; 462 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 463 return result; 464 465 if( shape == 0 ) 466 { 467 return location; 468 } 469 470 if( shape < 0 ) 471 { 472 skew_normal_distribution<RealType, Policy> D(0, 1, -shape); 473 result = mode(D); 474 result = location-scale*result; 475 return result; 476 } 477 478 // 21 elements 479 static const RealType shapes[] = { 480 0.0, 481 static_cast<RealType>(1.000000000000000e-004), 482 static_cast<RealType>(2.069138081114790e-004), 483 static_cast<RealType>(4.281332398719396e-004), 484 static_cast<RealType>(8.858667904100824e-004), 485 static_cast<RealType>(1.832980710832436e-003), 486 static_cast<RealType>(3.792690190732250e-003), 487 static_cast<RealType>(7.847599703514606e-003), 488 static_cast<RealType>(1.623776739188722e-002), 489 static_cast<RealType>(3.359818286283781e-002), 490 static_cast<RealType>(6.951927961775606e-002), 491 static_cast<RealType>(1.438449888287663e-001), 492 static_cast<RealType>(2.976351441631319e-001), 493 static_cast<RealType>(6.158482110660261e-001), 494 static_cast<RealType>(1.274274985703135e+000), 495 static_cast<RealType>(2.636650898730361e+000), 496 static_cast<RealType>(5.455594781168514e+000), 497 static_cast<RealType>(1.128837891684688e+001), 498 static_cast<RealType>(2.335721469090121e+001), 499 static_cast<RealType>(4.832930238571753e+001), 500 static_cast<RealType>(1.000000000000000e+002) 501 }; 502 503 // 21 elements 504 static const RealType guess[] = { 505 0.0, 506 static_cast<RealType>(5.000050000525391e-005), 507 static_cast<RealType>(1.500015000148736e-004), 508 static_cast<RealType>(3.500035000350010e-004), 509 static_cast<RealType>(7.500075000752560e-004), 510 static_cast<RealType>(1.450014500145258e-003), 511 static_cast<RealType>(3.050030500305390e-003), 512 static_cast<RealType>(6.250062500624765e-003), 513 static_cast<RealType>(1.295012950129504e-002), 514 static_cast<RealType>(2.675026750267495e-002), 515 static_cast<RealType>(5.525055250552491e-002), 516 static_cast<RealType>(1.132511325113255e-001), 517 static_cast<RealType>(2.249522495224952e-001), 518 static_cast<RealType>(3.992539925399257e-001), 519 static_cast<RealType>(5.353553535535358e-001), 520 static_cast<RealType>(4.954549545495457e-001), 521 static_cast<RealType>(3.524535245352451e-001), 522 static_cast<RealType>(2.182521825218249e-001), 523 static_cast<RealType>(1.256512565125654e-001), 524 static_cast<RealType>(6.945069450694508e-002), 525 static_cast<RealType>(3.735037350373460e-002) 526 }; 527 528 const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); 529 530 typedef typename std::iterator_traits<RealType*>::difference_type diff_type; 531 532 const diff_type d = std::distance(shapes, result_ptr); 533 534 BOOST_ASSERT(d > static_cast<diff_type>(0)); 535 536 // TODO: make the search bounds smarter, depending on the shape parameter 537 RealType search_min = 0; // below zero was caught above 538 RealType search_max = 0.55f; // will never go above 0.55 539 540 // refine 541 if(d < static_cast<diff_type>(21)) // shape smaller 100 542 { 543 // it is safe to assume that d > 0, because shape==0.0 is caught earlier 544 result = guess[d-static_cast<diff_type>(1)] 545 + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) 546 * (shape-shapes[d-static_cast<diff_type>(1)]); 547 } 548 else // shape greater 100 549 { 550 result = 1e-4f; 551 search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100 552 } 553 554 const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, 555 boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. 556 557 skew_normal_distribution<RealType, Policy> helper(0, 1, shape); 558 559 result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result, 560 search_min, search_max, get_digits, m); 561 562 result = result*scale + location; 563 564 return result; 565 } 566 567 568 569 template <class RealType, class Policy> skewness(const skew_normal_distribution<RealType,Policy> & dist)570 inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist) 571 { 572 BOOST_MATH_STD_USING // for ADL of std functions 573 using namespace boost::math::constants; 574 575 static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2); 576 const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); 577 578 return factor * pow(root_two_div_pi<RealType>() * delta, 3) / 579 pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5)); 580 } 581 582 template <class RealType, class Policy> kurtosis(const skew_normal_distribution<RealType,Policy> & dist)583 inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist) 584 { 585 return kurtosis_excess(dist)+static_cast<RealType>(3); 586 } 587 588 template <class RealType, class Policy> kurtosis_excess(const skew_normal_distribution<RealType,Policy> & dist)589 inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist) 590 { 591 using namespace boost::math::constants; 592 593 static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2); 594 595 const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); 596 597 const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta2; 598 const RealType y = two_div_pi<RealType>() * delta2; 599 600 return factor * y*y / (x*x); 601 } 602 603 namespace detail 604 { 605 606 template <class RealType, class Policy> 607 struct skew_normal_quantile_functor 608 { skew_normal_quantile_functorboost::math::detail::skew_normal_quantile_functor609 skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p) 610 : distribution(dist), prob(p) 611 { 612 } 613 operator ()boost::math::detail::skew_normal_quantile_functor614 boost::math::tuple<RealType, RealType> operator()(RealType const& x) 615 { 616 RealType c = cdf(distribution, x); 617 RealType fx = c - prob; // Difference cdf - value - to minimize. 618 RealType dx = pdf(distribution, x); // pdf is 1st derivative. 619 // return both function evaluation difference f(x) and 1st derivative f'(x). 620 return boost::math::make_tuple(fx, dx); 621 } 622 private: 623 const boost::math::skew_normal_distribution<RealType, Policy> distribution; 624 RealType prob; 625 }; 626 627 } // namespace detail 628 629 template <class RealType, class Policy> quantile(const skew_normal_distribution<RealType,Policy> & dist,const RealType & p)630 inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p) 631 { 632 const RealType scale = dist.scale(); 633 const RealType location = dist.location(); 634 const RealType shape = dist.shape(); 635 636 static const char* function = "boost::math::quantile(const skew_normal_distribution<%1%>&, %1%)"; 637 638 RealType result = 0; 639 if(false == detail::check_scale(function, scale, &result, Policy())) 640 return result; 641 if(false == detail::check_location(function, location, &result, Policy())) 642 return result; 643 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 644 return result; 645 if(false == detail::check_probability(function, p, &result, Policy())) 646 return result; 647 648 // Compute initial guess via Cornish-Fisher expansion. 649 RealType x = -boost::math::erfc_inv(2 * p, Policy()) * constants::root_two<RealType>(); 650 651 // Avoid unnecessary computations if there is no skew. 652 if(shape != 0) 653 { 654 const RealType skew = skewness(dist); 655 const RealType exk = kurtosis_excess(dist); 656 657 x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6) 658 + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24) 659 - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36); 660 } // if(shape != 0) 661 662 result = standard_deviation(dist)*x+mean(dist); 663 664 // handle special case of non-skew normal distribution. 665 if(shape == 0) 666 return result; 667 668 // refine the result by numerically searching the root of (p-cdf) 669 670 const RealType search_min = range(dist).first; 671 const RealType search_max = range(dist).second; 672 673 const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, 674 boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. 675 676 result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(dist, p), result, 677 search_min, search_max, get_digits, m); 678 679 return result; 680 } // quantile 681 682 template <class RealType, class Policy> quantile(const complemented2_type<skew_normal_distribution<RealType,Policy>,RealType> & c)683 inline RealType quantile(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) 684 { 685 const RealType scale = c.dist.scale(); 686 const RealType location = c.dist.location(); 687 const RealType shape = c.dist.shape(); 688 689 static const char* function = "boost::math::quantile(const complement(skew_normal_distribution<%1%>&), %1%)"; 690 RealType result = 0; 691 if(false == detail::check_scale(function, scale, &result, Policy())) 692 return result; 693 if(false == detail::check_location(function, location, &result, Policy())) 694 return result; 695 if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) 696 return result; 697 RealType q = c.param; 698 if(false == detail::check_probability(function, q, &result, Policy())) 699 return result; 700 701 skew_normal_distribution<RealType, Policy> D(-location, scale, -shape); 702 703 result = -quantile(D, q); 704 705 return result; 706 } // quantile 707 708 709 } // namespace math 710 } // namespace boost 711 712 // This include must be at the end, *after* the accessors 713 // for this distribution have been defined, in order to 714 // keep compilers that support two-phase lookup happy. 715 #include <boost/math/distributions/detail/derived_accessors.hpp> 716 717 #endif // BOOST_STATS_SKEW_NORMAL_HPP 718 719 720