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