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
23namespace boost{ namespace math{ namespace detail{
24
25template <class T, class Policy>
26struct 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 }
33private:
34 T z, p;
35 bool invert;
36};
37
38template <class T, class Policy>
39T 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
69template <class T, class Policy>
70T 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
153template <class T1, class T2, class Policy>
154inline 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
183template <class T1, class T2, class Policy>
184inline 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
213template <class T1, class T2>
214inline 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
220template <class T1, class T2>
221inline 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

source code of include/boost/math/special_functions/detail/gamma_inva.hpp