Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:39:44

0001 //  Copyright Benjamin Sobotta 2012
0002 
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 
0007 #ifndef BOOST_STATS_SKEW_NORMAL_HPP
0008 #define BOOST_STATS_SKEW_NORMAL_HPP
0009 
0010 // http://en.wikipedia.org/wiki/Skew_normal_distribution
0011 // http://azzalini.stat.unipd.it/SN/
0012 // Also:
0013 // Azzalini, A. (1985). "A class of distributions which includes the normal ones".
0014 // Scand. J. Statist. 12: 171-178.
0015 
0016 #include <boost/math/distributions/fwd.hpp> // TODO add skew_normal distribution to fwd.hpp!
0017 #include <boost/math/special_functions/owens_t.hpp> // Owen's T function
0018 #include <boost/math/distributions/complement.hpp>
0019 #include <boost/math/distributions/normal.hpp>
0020 #include <boost/math/distributions/detail/common_error_handling.hpp>
0021 #include <boost/math/constants/constants.hpp>
0022 #include <boost/math/tools/tuple.hpp>
0023 #include <boost/math/tools/roots.hpp> // Newton-Raphson
0024 #include <boost/math/tools/assert.hpp>
0025 #include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder.
0026 
0027 #include <utility>
0028 #include <algorithm> // std::lower_bound, std::distance
0029 
0030 namespace boost{ namespace math{
0031 
0032   namespace detail
0033   {
0034     template <class RealType, class Policy>
0035     inline bool check_skew_normal_shape(
0036       const char* function,
0037       RealType shape,
0038       RealType* result,
0039       const Policy& pol)
0040     {
0041       if(!(boost::math::isfinite)(shape))
0042       {
0043         *result =
0044           policies::raise_domain_error<RealType>(function,
0045           "Shape parameter is %1%, but must be finite!",
0046           shape, pol);
0047         return false;
0048       }
0049       return true;
0050     }
0051 
0052   } // namespace detail
0053 
0054   template <class RealType = double, class Policy = policies::policy<> >
0055   class skew_normal_distribution
0056   {
0057   public:
0058     typedef RealType value_type;
0059     typedef Policy policy_type;
0060 
0061     skew_normal_distribution(RealType l_location = 0, RealType l_scale = 1, RealType l_shape = 0)
0062       : location_(l_location), scale_(l_scale), shape_(l_shape)
0063     { // Default is a 'standard' normal distribution N01. (shape=0 results in the normal distribution with no skew)
0064       static const char* function = "boost::math::skew_normal_distribution<%1%>::skew_normal_distribution";
0065 
0066       RealType result;
0067       detail::check_scale(function, l_scale, &result, Policy());
0068       detail::check_location(function, l_location, &result, Policy());
0069       detail::check_skew_normal_shape(function, l_shape, &result, Policy());
0070     }
0071 
0072     RealType location()const
0073     {
0074       return location_;
0075     }
0076 
0077     RealType scale()const
0078     {
0079       return scale_;
0080     }
0081 
0082     RealType shape()const
0083     {
0084       return shape_;
0085     }
0086 
0087 
0088   private:
0089     //
0090     // Data members:
0091     //
0092     RealType location_;  // distribution location.
0093     RealType scale_;    // distribution scale.
0094     RealType shape_;    // distribution shape.
0095   }; // class skew_normal_distribution
0096 
0097   typedef skew_normal_distribution<double> skew_normal;
0098 
0099   #ifdef __cpp_deduction_guides
0100   template <class RealType>
0101   skew_normal_distribution(RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0102   template <class RealType>
0103   skew_normal_distribution(RealType,RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0104   template <class RealType>
0105   skew_normal_distribution(RealType,RealType,RealType)->skew_normal_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0106   #endif
0107 
0108   template <class RealType, class Policy>
0109   inline const std::pair<RealType, RealType> range(const skew_normal_distribution<RealType, Policy>& /*dist*/)
0110   { // Range of permissible values for random variable x.
0111     using boost::math::tools::max_value;
0112     return std::pair<RealType, RealType>(
0113        std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(),
0114        std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); // - to + max value.
0115   }
0116 
0117   template <class RealType, class Policy>
0118   inline const std::pair<RealType, RealType> support(const skew_normal_distribution<RealType, Policy>& /*dist*/)
0119   { // Range of supported values for random variable x.
0120     // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0121 
0122     using boost::math::tools::max_value;
0123     return std::pair<RealType, RealType>(-max_value<RealType>(),  max_value<RealType>()); // - to + max value.
0124   }
0125 
0126   template <class RealType, class Policy>
0127   inline RealType pdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x)
0128   {
0129     const RealType scale = dist.scale();
0130     const RealType location = dist.location();
0131     const RealType shape = dist.shape();
0132 
0133     static const char* function = "boost::math::pdf(const skew_normal_distribution<%1%>&, %1%)";
0134 
0135     RealType result = 0;
0136     if(false == detail::check_scale(function, scale, &result, Policy()))
0137     {
0138       return result;
0139     }
0140     if(false == detail::check_location(function, location, &result, Policy()))
0141     {
0142       return result;
0143     }
0144     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0145     {
0146       return result;
0147     }
0148     if((boost::math::isinf)(x))
0149     {
0150        return 0; // pdf + and - infinity is zero.
0151     }
0152     // Below produces MSVC 4127 warnings, so the above used instead.
0153     //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())
0154     //{ // pdf + and - infinity is zero.
0155     //  return 0;
0156     //}
0157     if(false == detail::check_x(function, x, &result, Policy()))
0158     {
0159       return result;
0160     }
0161 
0162     const RealType transformed_x = (x-location)/scale;
0163 
0164     normal_distribution<RealType, Policy> std_normal;
0165 
0166     result = pdf(std_normal, transformed_x) * cdf(std_normal, shape*transformed_x) * 2 / scale;
0167 
0168     return result;
0169   } // pdf
0170 
0171   template <class RealType, class Policy>
0172   inline RealType cdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x)
0173   {
0174     const RealType scale = dist.scale();
0175     const RealType location = dist.location();
0176     const RealType shape = dist.shape();
0177 
0178     static const char* function = "boost::math::cdf(const skew_normal_distribution<%1%>&, %1%)";
0179     RealType result = 0;
0180     if(false == detail::check_scale(function, scale, &result, Policy()))
0181     {
0182       return result;
0183     }
0184     if(false == detail::check_location(function, location, &result, Policy()))
0185     {
0186       return result;
0187     }
0188     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0189     {
0190       return result;
0191     }
0192     if((boost::math::isinf)(x))
0193     {
0194       if(x < 0) return 0; // -infinity
0195       return 1; // + infinity
0196     }
0197     // These produce MSVC 4127 warnings, so the above used instead.
0198     //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
0199     //{ // cdf +infinity is unity.
0200     //  return 1;
0201     //}
0202     //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
0203     //{ // cdf -infinity is zero.
0204     //  return 0;
0205     //}
0206     if(false == detail::check_x(function, x, &result, Policy()))
0207     {
0208       return result;
0209     }
0210 
0211     const RealType transformed_x = (x-location)/scale;
0212 
0213     normal_distribution<RealType, Policy> std_normal;
0214 
0215     result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2);
0216 
0217     return result;
0218   } // cdf
0219 
0220   template <class RealType, class Policy>
0221   inline RealType cdf(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c)
0222   {
0223     const RealType scale = c.dist.scale();
0224     const RealType location = c.dist.location();
0225     const RealType shape = c.dist.shape();
0226     const RealType x = c.param;
0227 
0228     static const char* function = "boost::math::cdf(const complement(skew_normal_distribution<%1%>&), %1%)";
0229 
0230     if((boost::math::isinf)(x))
0231     {
0232       if(x < 0) return 1; // cdf complement -infinity is unity.
0233       return 0; // cdf complement +infinity is zero
0234     }
0235     // These produce MSVC 4127 warnings, so the above used instead.
0236     //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
0237     //{ // cdf complement +infinity is zero.
0238     //  return 0;
0239     //}
0240     //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
0241     //{ // cdf complement -infinity is unity.
0242     //  return 1;
0243     //}
0244     RealType result = 0;
0245     if(false == detail::check_scale(function, scale, &result, Policy()))
0246       return result;
0247     if(false == detail::check_location(function, location, &result, Policy()))
0248       return result;
0249     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0250       return result;
0251     if(false == detail::check_x(function, x, &result, Policy()))
0252       return result;
0253 
0254     const RealType transformed_x = (x-location)/scale;
0255 
0256     normal_distribution<RealType, Policy> std_normal;
0257 
0258     result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2);
0259     return result;
0260   } // cdf complement
0261 
0262   template <class RealType, class Policy>
0263   inline RealType location(const skew_normal_distribution<RealType, Policy>& dist)
0264   {
0265     return dist.location();
0266   }
0267 
0268   template <class RealType, class Policy>
0269   inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist)
0270   {
0271     return dist.scale();
0272   }
0273 
0274   template <class RealType, class Policy>
0275   inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist)
0276   {
0277     return dist.shape();
0278   }
0279 
0280   template <class RealType, class Policy>
0281   inline RealType mean(const skew_normal_distribution<RealType, Policy>& dist)
0282   {
0283     BOOST_MATH_STD_USING  // for ADL of std functions
0284 
0285     using namespace boost::math::constants;
0286 
0287     //const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
0288 
0289     //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>();
0290 
0291     return dist.location() + dist.scale() * dist.shape() / sqrt(pi<RealType>()+pi<RealType>()*dist.shape()*dist.shape()) * root_two<RealType>();
0292   }
0293 
0294   template <class RealType, class Policy>
0295   inline RealType variance(const skew_normal_distribution<RealType, Policy>& dist)
0296   {
0297     using namespace boost::math::constants;
0298 
0299     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);
0300     //const RealType inv_delta2 = static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape());
0301 
0302     RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta2);
0303     //RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()/inv_delta2);
0304 
0305     return variance;
0306   }
0307 
0308   namespace detail
0309   {
0310     /*
0311       TODO No closed expression for mode, so use max of pdf.
0312     */
0313 
0314     template <class RealType, class Policy>
0315     inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist)
0316     { // mode.
0317         static const char* function = "mode(skew_normal_distribution<%1%> const&)";
0318         const RealType scale = dist.scale();
0319         const RealType location = dist.location();
0320         const RealType shape = dist.shape();
0321 
0322         RealType result;
0323         if(!detail::check_scale(
0324           function,
0325           scale, &result, Policy())
0326           ||
0327         !detail::check_skew_normal_shape(
0328           function,
0329           shape,
0330           &result,
0331           Policy()))
0332         return result;
0333 
0334         if( shape == 0 )
0335         {
0336           return location;
0337         }
0338 
0339         if( shape < 0 )
0340         {
0341           skew_normal_distribution<RealType, Policy> D(0, 1, -shape);
0342           result = mode_fallback(D);
0343           result = location-scale*result;
0344           return result;
0345         }
0346 
0347         BOOST_MATH_STD_USING
0348 
0349         // 21 elements
0350         static const RealType shapes[] = {
0351           0.0,
0352           1.000000000000000e-004,
0353           2.069138081114790e-004,
0354           4.281332398719396e-004,
0355           8.858667904100824e-004,
0356           1.832980710832436e-003,
0357           3.792690190732250e-003,
0358           7.847599703514606e-003,
0359           1.623776739188722e-002,
0360           3.359818286283781e-002,
0361           6.951927961775606e-002,
0362           1.438449888287663e-001,
0363           2.976351441631319e-001,
0364           6.158482110660261e-001,
0365           1.274274985703135e+000,
0366           2.636650898730361e+000,
0367           5.455594781168514e+000,
0368           1.128837891684688e+001,
0369           2.335721469090121e+001,
0370           4.832930238571753e+001,
0371           1.000000000000000e+002};
0372 
0373         // 21 elements
0374         static const RealType guess[] = {
0375           0.0,
0376           5.000050000525391e-005,
0377           1.500015000148736e-004,
0378           3.500035000350010e-004,
0379           7.500075000752560e-004,
0380           1.450014500145258e-003,
0381           3.050030500305390e-003,
0382           6.250062500624765e-003,
0383           1.295012950129504e-002,
0384           2.675026750267495e-002,
0385           5.525055250552491e-002,
0386           1.132511325113255e-001,
0387           2.249522495224952e-001,
0388           3.992539925399257e-001,
0389           5.353553535535358e-001,
0390           4.954549545495457e-001,
0391           3.524535245352451e-001,
0392           2.182521825218249e-001,
0393           1.256512565125654e-001,
0394           6.945069450694508e-002,
0395           3.735037350373460e-002
0396         };
0397 
0398         const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);
0399 
0400         typedef typename std::iterator_traits<RealType*>::difference_type diff_type;
0401 
0402         const diff_type d = std::distance(shapes, result_ptr);
0403 
0404         BOOST_MATH_ASSERT(d > static_cast<diff_type>(0));
0405 
0406         // refine
0407         if(d < static_cast<diff_type>(21)) // shape smaller 100
0408         {
0409           result = guess[d-static_cast<diff_type>(1)]
0410             + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])
0411             * (shape-shapes[d-static_cast<diff_type>(1)]);
0412         }
0413         else // shape greater 100
0414         {
0415           result = 1e-4;
0416         }
0417 
0418         skew_normal_distribution<RealType, Policy> helper(0, 1, shape);
0419 
0420         result = detail::generic_find_mode_01(helper, result, function);
0421 
0422         result = result*scale + location;
0423 
0424         return result;
0425     } // mode_fallback
0426 
0427 
0428     /*
0429      * TODO No closed expression for mode, so use f'(x) = 0
0430      */
0431     template <class RealType, class Policy>
0432     struct skew_normal_mode_functor
0433     {
0434       skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist)
0435         : distribution(dist)
0436       {
0437       }
0438 
0439       boost::math::tuple<RealType, RealType> operator()(RealType const& x)
0440       {
0441         normal_distribution<RealType, Policy> std_normal;
0442         const RealType shape = distribution.shape();
0443         const RealType pdf_x = pdf(distribution, x);
0444         const RealType normpdf_x = pdf(std_normal, x);
0445         const RealType normpdf_ax = pdf(std_normal, x*shape);
0446         RealType fx = static_cast<RealType>(2)*shape*normpdf_ax*normpdf_x - x*pdf_x;
0447         RealType dx = static_cast<RealType>(2)*shape*x*normpdf_x*normpdf_ax*(static_cast<RealType>(1) + shape*shape) + pdf_x + x*fx;
0448         // return both function evaluation difference f(x) and 1st derivative f'(x).
0449         return boost::math::make_tuple(fx, -dx);
0450       }
0451     private:
0452       const boost::math::skew_normal_distribution<RealType, Policy> distribution;
0453     };
0454 
0455   } // namespace detail
0456 
0457   template <class RealType, class Policy>
0458   inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist)
0459   {
0460     const RealType scale = dist.scale();
0461     const RealType location = dist.location();
0462     const RealType shape = dist.shape();
0463 
0464     static const char* function = "boost::math::mode(const skew_normal_distribution<%1%>&, %1%)";
0465 
0466     RealType result = 0;
0467     if(false == detail::check_scale(function, scale, &result, Policy()))
0468       return result;
0469     if(false == detail::check_location(function, location, &result, Policy()))
0470       return result;
0471     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0472       return result;
0473 
0474     if( shape == 0 )
0475     {
0476       return location;
0477     }
0478 
0479     if( shape < 0 )
0480     {
0481       skew_normal_distribution<RealType, Policy> D(0, 1, -shape);
0482       result = mode(D);
0483       result = location-scale*result;
0484       return result;
0485     }
0486 
0487     // 21 elements
0488     static const RealType shapes[] = {
0489       static_cast<RealType>(0.0),
0490       static_cast<RealType>(1.000000000000000e-004),
0491       static_cast<RealType>(2.069138081114790e-004),
0492       static_cast<RealType>(4.281332398719396e-004),
0493       static_cast<RealType>(8.858667904100824e-004),
0494       static_cast<RealType>(1.832980710832436e-003),
0495       static_cast<RealType>(3.792690190732250e-003),
0496       static_cast<RealType>(7.847599703514606e-003),
0497       static_cast<RealType>(1.623776739188722e-002),
0498       static_cast<RealType>(3.359818286283781e-002),
0499       static_cast<RealType>(6.951927961775606e-002),
0500       static_cast<RealType>(1.438449888287663e-001),
0501       static_cast<RealType>(2.976351441631319e-001),
0502       static_cast<RealType>(6.158482110660261e-001),
0503       static_cast<RealType>(1.274274985703135e+000),
0504       static_cast<RealType>(2.636650898730361e+000),
0505       static_cast<RealType>(5.455594781168514e+000),
0506       static_cast<RealType>(1.128837891684688e+001),
0507       static_cast<RealType>(2.335721469090121e+001),
0508       static_cast<RealType>(4.832930238571753e+001),
0509       static_cast<RealType>(1.000000000000000e+002)
0510     };
0511 
0512     // 21 elements
0513     static const RealType guess[] = {
0514       static_cast<RealType>(0.0),
0515       static_cast<RealType>(5.000050000525391e-005),
0516       static_cast<RealType>(1.500015000148736e-004),
0517       static_cast<RealType>(3.500035000350010e-004),
0518       static_cast<RealType>(7.500075000752560e-004),
0519       static_cast<RealType>(1.450014500145258e-003),
0520       static_cast<RealType>(3.050030500305390e-003),
0521       static_cast<RealType>(6.250062500624765e-003),
0522       static_cast<RealType>(1.295012950129504e-002),
0523       static_cast<RealType>(2.675026750267495e-002),
0524       static_cast<RealType>(5.525055250552491e-002),
0525       static_cast<RealType>(1.132511325113255e-001),
0526       static_cast<RealType>(2.249522495224952e-001),
0527       static_cast<RealType>(3.992539925399257e-001),
0528       static_cast<RealType>(5.353553535535358e-001),
0529       static_cast<RealType>(4.954549545495457e-001),
0530       static_cast<RealType>(3.524535245352451e-001),
0531       static_cast<RealType>(2.182521825218249e-001),
0532       static_cast<RealType>(1.256512565125654e-001),
0533       static_cast<RealType>(6.945069450694508e-002),
0534       static_cast<RealType>(3.735037350373460e-002)
0535     };
0536 
0537     const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape);
0538 
0539     typedef typename std::iterator_traits<RealType*>::difference_type diff_type;
0540 
0541     const diff_type d = std::distance(shapes, result_ptr);
0542 
0543     BOOST_MATH_ASSERT(d > static_cast<diff_type>(0));
0544 
0545     // TODO: make the search bounds smarter, depending on the shape parameter
0546     RealType search_min = 0; // below zero was caught above
0547     RealType search_max = 0.55f; // will never go above 0.55
0548 
0549     // refine
0550     if(d < static_cast<diff_type>(21)) // shape smaller 100
0551     {
0552       // it is safe to assume that d > 0, because shape==0.0 is caught earlier
0553       result = guess[d-static_cast<diff_type>(1)]
0554         + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)])
0555         * (shape-shapes[d-static_cast<diff_type>(1)]);
0556     }
0557     else // shape greater 100
0558     {
0559       result = 1e-4f;
0560       search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100
0561     }
0562 
0563     const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
0564     std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
0565 
0566     skew_normal_distribution<RealType, Policy> helper(0, 1, shape);
0567 
0568     result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result,
0569       search_min, search_max, get_digits, m);
0570 
0571     result = result*scale + location;
0572 
0573     return result;
0574   }
0575 
0576 
0577 
0578   template <class RealType, class Policy>
0579   inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist)
0580   {
0581     BOOST_MATH_STD_USING  // for ADL of std functions
0582     using namespace boost::math::constants;
0583 
0584     static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2);
0585     const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape());
0586 
0587     return static_cast<RealType>(factor * pow(root_two_div_pi<RealType>() * delta, 3) /
0588       pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5)));
0589   }
0590 
0591   template <class RealType, class Policy>
0592   inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist)
0593   {
0594     return kurtosis_excess(dist)+static_cast<RealType>(3);
0595   }
0596 
0597   template <class RealType, class Policy>
0598   inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist)
0599   {
0600     using namespace boost::math::constants;
0601 
0602     static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2);
0603 
0604     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);
0605 
0606     const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta2;
0607     const RealType y = two_div_pi<RealType>() * delta2;
0608 
0609     return factor * y*y / (x*x);
0610   }
0611 
0612   namespace detail
0613   {
0614 
0615     template <class RealType, class Policy>
0616     struct skew_normal_quantile_functor
0617     {
0618       skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p)
0619         : distribution(dist), prob(p)
0620       {
0621       }
0622 
0623       boost::math::tuple<RealType, RealType> operator()(RealType const& x)
0624       {
0625         RealType c = cdf(distribution, x);
0626         RealType fx = c - prob;  // Difference cdf - value - to minimize.
0627         RealType dx = pdf(distribution, x); // pdf is 1st derivative.
0628         // return both function evaluation difference f(x) and 1st derivative f'(x).
0629         return boost::math::make_tuple(fx, dx);
0630       }
0631     private:
0632       const boost::math::skew_normal_distribution<RealType, Policy> distribution;
0633       RealType prob;
0634     };
0635 
0636   } // namespace detail
0637 
0638   template <class RealType, class Policy>
0639   inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p)
0640   {
0641     const RealType scale = dist.scale();
0642     const RealType location = dist.location();
0643     const RealType shape = dist.shape();
0644 
0645     static const char* function = "boost::math::quantile(const skew_normal_distribution<%1%>&, %1%)";
0646 
0647     RealType result = 0;
0648     if(false == detail::check_scale(function, scale, &result, Policy()))
0649       return result;
0650     if(false == detail::check_location(function, location, &result, Policy()))
0651       return result;
0652     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0653       return result;
0654     if(false == detail::check_probability(function, p, &result, Policy()))
0655       return result;
0656 
0657     // Compute initial guess via Cornish-Fisher expansion.
0658     RealType x = -boost::math::erfc_inv(2 * p, Policy()) * constants::root_two<RealType>();
0659 
0660     // Avoid unnecessary computations if there is no skew.
0661     if(shape != 0)
0662     {
0663       const RealType skew = skewness(dist);
0664       const RealType exk = kurtosis_excess(dist);
0665 
0666       x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6)
0667       + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24)
0668       - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36);
0669     } // if(shape != 0)
0670 
0671     result = standard_deviation(dist)*x+mean(dist);
0672 
0673     // handle special case of non-skew normal distribution.
0674     if(shape == 0)
0675       return result;
0676 
0677     // refine the result by numerically searching the root of (p-cdf)
0678 
0679     const RealType search_min = support(dist).first;
0680     const RealType search_max = support(dist).second;
0681 
0682     const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
0683     std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
0684 
0685     result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(dist, p), result,
0686       search_min, search_max, get_digits, m);
0687 
0688     return result;
0689   } // quantile
0690 
0691   template <class RealType, class Policy>
0692   inline RealType quantile(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c)
0693   {
0694     const RealType scale = c.dist.scale();
0695     const RealType location = c.dist.location();
0696     const RealType shape = c.dist.shape();
0697 
0698     static const char* function = "boost::math::quantile(const complement(skew_normal_distribution<%1%>&), %1%)";
0699     RealType result = 0;
0700     if(false == detail::check_scale(function, scale, &result, Policy()))
0701       return result;
0702     if(false == detail::check_location(function, location, &result, Policy()))
0703       return result;
0704     if(false == detail::check_skew_normal_shape(function, shape, &result, Policy()))
0705       return result;
0706     RealType q = c.param;
0707     if(false == detail::check_probability(function, q, &result, Policy()))
0708       return result;
0709 
0710     skew_normal_distribution<RealType, Policy> D(-location, scale, -shape);
0711 
0712     result = -quantile(D, q);
0713 
0714     return result;
0715   } // quantile
0716 
0717 
0718 } // namespace math
0719 } // namespace boost
0720 
0721 // This include must be at the end, *after* the accessors
0722 // for this distribution have been defined, in order to
0723 // keep compilers that support two-phase lookup happy.
0724 #include <boost/math/distributions/detail/derived_accessors.hpp>
0725 
0726 #endif // BOOST_STATS_SKEW_NORMAL_HPP
0727 
0728