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