]>
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 | |
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 | ||
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 |