skew_normal.hpp 25 KB

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