]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
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 | ||
6 | // | |
7 | // This is not a complete header file, it is included by gamma.hpp | |
8 | // after it has defined it's definitions. This inverts the incomplete | |
9 | // gamma functions P and Q on the first parameter "a" using a generic | |
10 | // root finding algorithm (TOMS Algorithm 748). | |
11 | // | |
12 | ||
13 | #ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA | |
14 | #define BOOST_MATH_SP_DETAIL_GAMMA_INVA | |
15 | ||
16 | #ifdef _MSC_VER | |
17 | #pragma once | |
18 | #endif | |
19 | ||
20 | #include <boost/math/tools/toms748_solve.hpp> | |
21 | #include <boost/cstdint.hpp> | |
22 | ||
23 | namespace boost{ namespace math{ namespace detail{ | |
24 | ||
25 | template <class T, class Policy> | |
26 | struct gamma_inva_t | |
27 | { | |
28 | gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {} | |
29 | T operator()(T a) | |
30 | { | |
31 | return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p; | |
32 | } | |
33 | private: | |
34 | T z, p; | |
35 | bool invert; | |
36 | }; | |
37 | ||
38 | template <class T, class Policy> | |
39 | T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol) | |
40 | { | |
41 | BOOST_MATH_STD_USING | |
42 | // mean: | |
43 | T m = lambda; | |
44 | // standard deviation: | |
45 | T sigma = sqrt(lambda); | |
46 | // skewness | |
47 | T sk = 1 / sigma; | |
48 | // kurtosis: | |
49 | // T k = 1/lambda; | |
50 | // Get the inverse of a std normal distribution: | |
51 | T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>(); | |
52 | // Set the sign: | |
53 | if(p < 0.5) | |
54 | x = -x; | |
55 | T x2 = x * x; | |
56 | // w is correction term due to skewness | |
57 | T w = x + sk * (x2 - 1) / 6; | |
58 | /* | |
59 | // Add on correction due to kurtosis. | |
60 | // Disabled for now, seems to make things worse? | |
61 | // | |
62 | if(lambda >= 10) | |
63 | w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36; | |
64 | */ | |
65 | w = m + sigma * w; | |
66 | return w > tools::min_value<T>() ? w : tools::min_value<T>(); | |
67 | } | |
68 | ||
69 | template <class T, class Policy> | |
70 | T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol) | |
71 | { | |
72 | BOOST_MATH_STD_USING // for ADL of std lib math functions | |
73 | // | |
74 | // Special cases first: | |
75 | // | |
76 | if(p == 0) | |
77 | { | |
78 | return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy()); | |
79 | } | |
80 | if(q == 0) | |
81 | { | |
82 | return tools::min_value<T>(); | |
83 | } | |
84 | // | |
85 | // Function object, this is the functor whose root | |
86 | // we have to solve: | |
87 | // | |
88 | gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true); | |
89 | // | |
90 | // Tolerance: full precision. | |
91 | // | |
92 | tools::eps_tolerance<T> tol(policies::digits<T, Policy>()); | |
93 | // | |
94 | // Now figure out a starting guess for what a may be, | |
95 | // we'll start out with a value that'll put p or q | |
96 | // right bang in the middle of their range, the functions | |
97 | // are quite sensitive so we should need too many steps | |
98 | // to bracket the root from there: | |
99 | // | |
100 | T guess; | |
101 | T factor = 8; | |
102 | if(z >= 1) | |
103 | { | |
104 | // | |
105 | // We can use the relationship between the incomplete | |
106 | // gamma function and the poisson distribution to | |
107 | // calculate an approximate inverse, for large z | |
108 | // this is actually pretty accurate, but it fails badly | |
109 | // when z is very small. Also set our step-factor according | |
110 | // to how accurate we think the result is likely to be: | |
111 | // | |
112 | guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol); | |
113 | if(z > 5) | |
114 | { | |
115 | if(z > 1000) | |
116 | factor = 1.01f; | |
117 | else if(z > 50) | |
118 | factor = 1.1f; | |
119 | else if(guess > 10) | |
120 | factor = 1.25f; | |
121 | else | |
122 | factor = 2; | |
123 | if(guess < 1.1) | |
124 | factor = 8; | |
125 | } | |
126 | } | |
127 | else if(z > 0.5) | |
128 | { | |
129 | guess = z * 1.2f; | |
130 | } | |
131 | else | |
132 | { | |
133 | guess = -0.4f / log(z); | |
134 | } | |
135 | // | |
136 | // Max iterations permitted: | |
137 | // | |
138 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
139 | // | |
140 | // Use our generic derivative-free root finding procedure. | |
141 | // We could use Newton steps here, taking the PDF of the | |
142 | // Poisson distribution as our derivative, but that's | |
143 | // even worse performance-wise than the generic method :-( | |
144 | // | |
145 | std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol); | |
146 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
147 | return policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol); | |
148 | return (r.first + r.second) / 2; | |
149 | } | |
150 | ||
151 | } // namespace detail | |
152 | ||
153 | template <class T1, class T2, class Policy> | |
154 | inline typename tools::promote_args<T1, T2>::type | |
155 | gamma_p_inva(T1 x, T2 p, const Policy& pol) | |
156 | { | |
157 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
158 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
159 | typedef typename policies::normalise< | |
160 | Policy, | |
161 | policies::promote_float<false>, | |
162 | policies::promote_double<false>, | |
163 | policies::discrete_quantile<>, | |
164 | policies::assert_undefined<> >::type forwarding_policy; | |
165 | ||
166 | if(p == 0) | |
167 | { | |
168 | policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy()); | |
169 | } | |
170 | if(p == 1) | |
171 | { | |
172 | return tools::min_value<result_type>(); | |
173 | } | |
174 | ||
175 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
176 | detail::gamma_inva_imp( | |
177 | static_cast<value_type>(x), | |
178 | static_cast<value_type>(p), | |
179 | static_cast<value_type>(1 - static_cast<value_type>(p)), | |
180 | pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)"); | |
181 | } | |
182 | ||
183 | template <class T1, class T2, class Policy> | |
184 | inline typename tools::promote_args<T1, T2>::type | |
185 | gamma_q_inva(T1 x, T2 q, const Policy& pol) | |
186 | { | |
187 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
188 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
189 | typedef typename policies::normalise< | |
190 | Policy, | |
191 | policies::promote_float<false>, | |
192 | policies::promote_double<false>, | |
193 | policies::discrete_quantile<>, | |
194 | policies::assert_undefined<> >::type forwarding_policy; | |
195 | ||
196 | if(q == 1) | |
197 | { | |
198 | policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", 0, Policy()); | |
199 | } | |
200 | if(q == 0) | |
201 | { | |
202 | return tools::min_value<result_type>(); | |
203 | } | |
204 | ||
205 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
206 | detail::gamma_inva_imp( | |
207 | static_cast<value_type>(x), | |
208 | static_cast<value_type>(1 - static_cast<value_type>(q)), | |
209 | static_cast<value_type>(q), | |
210 | pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)"); | |
211 | } | |
212 | ||
213 | template <class T1, class T2> | |
214 | inline typename tools::promote_args<T1, T2>::type | |
215 | gamma_p_inva(T1 x, T2 p) | |
216 | { | |
217 | return boost::math::gamma_p_inva(x, p, policies::policy<>()); | |
218 | } | |
219 | ||
220 | template <class T1, class T2> | |
221 | inline typename tools::promote_args<T1, T2>::type | |
222 | gamma_q_inva(T1 x, T2 q) | |
223 | { | |
224 | return boost::math::gamma_q_inva(x, q, policies::policy<>()); | |
225 | } | |
226 | ||
227 | } // namespace math | |
228 | } // namespace boost | |
229 | ||
230 | #endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA | |
231 | ||
232 | ||
233 |