ceph/src/boost/boost/math/distributions/skew_normal.hpp

   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>
  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
  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
  72     RealType location()const
  73     {
  74       return location_;
  75     }
  76
  77     RealType scale()const
  78     {
  79       return scale_;
  80     }
  81
  82     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>
 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>
 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>
 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
 126     RealType result = 0;
 127     if(false == detail::check_scale(function, scale, &result, Policy()))
 128     {
 129       return result;
 130     }
 131     if(false == detail::check_location(function, location, &result, Policy()))
 132     {
 133       return result;
 134     }
 135     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
 136     {
 137       return result;
 138     }
 139     if((boost::math::isinf)(x))
 140     {
 141        return 0; // pdf + and - infinity is zero.
 142     }
 143     // Below produces MSVC 4127 warnings, so the above used instead.
 144     //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())
 145     //{ // pdf + and - infinity is zero.
 146     //  return 0;
 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>
 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>
 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>
 254   inline RealType location(const skew_normal_distribution<RealType, Policy>& dist)
 255   {
 256     return dist.location();
 257   }
 258
 259   template <class RealType, class Policy>
 260   inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist)
 261   {
 262     return dist.scale();
 263   }
 264
 265   template <class RealType, class Policy>
 266   inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist)
 267   {
 268     return dist.shape();
 269   }
 270
 271   template <class RealType, class Policy>
 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>
 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>
 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     {
 425       skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist)
 426         : distribution(dist)
 427       {
 428       }
 429
 430       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>
 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>
 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>
 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>
 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     {
 609       skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p)
 610         : distribution(dist), prob(p)
 611       {
 612       }
 613
 614       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>
 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>
 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