]>
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 beta.hpp | |
8 | // after it has defined it's definitions. This inverts the incomplete | |
9 | // beta functions ibeta and ibetac on the first parameters "a" | |
10 | // and "b" using a generic root finding algorithm (TOMS Algorithm 748). | |
11 | // | |
12 | ||
13 | #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB | |
14 | #define BOOST_MATH_SP_DETAIL_BETA_INV_AB | |
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 beta_inv_ab_t | |
27 | { | |
28 | beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {} | |
29 | T operator()(T a) | |
30 | { | |
31 | return invert ? | |
32 | p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) | |
33 | : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p; | |
34 | } | |
35 | private: | |
36 | T b, z, p; | |
37 | bool invert, swap_ab; | |
38 | }; | |
39 | ||
40 | template <class T, class Policy> | |
41 | T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol) | |
42 | { | |
43 | BOOST_MATH_STD_USING | |
44 | // mean: | |
45 | T m = n * (sfc) / sf; | |
46 | T t = sqrt(n * (sfc)); | |
47 | // standard deviation: | |
48 | T sigma = t / sf; | |
49 | // skewness | |
50 | T sk = (1 + sfc) / t; | |
51 | // kurtosis: | |
52 | T k = (6 - sf * (5+sfc)) / (n * (sfc)); | |
53 | // Get the inverse of a std normal distribution: | |
54 | T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>(); | |
55 | // Set the sign: | |
56 | if(p < 0.5) | |
57 | x = -x; | |
58 | T x2 = x * x; | |
59 | // w is correction term due to skewness | |
60 | T w = x + sk * (x2 - 1) / 6; | |
61 | // | |
62 | // Add on correction due to kurtosis. | |
63 | // | |
64 | if(n >= 10) | |
65 | w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36; | |
66 | ||
67 | w = m + sigma * w; | |
68 | if(w < tools::min_value<T>()) | |
69 | return tools::min_value<T>(); | |
70 | return w; | |
71 | } | |
72 | ||
73 | template <class T, class Policy> | |
74 | T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol) | |
75 | { | |
76 | BOOST_MATH_STD_USING // for ADL of std lib math functions | |
77 | // | |
78 | // Special cases first: | |
79 | // | |
80 | BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab); | |
81 | if(p == 0) | |
82 | { | |
83 | return swap_ab ? tools::min_value<T>() : tools::max_value<T>(); | |
84 | } | |
85 | if(q == 0) | |
86 | { | |
87 | return swap_ab ? tools::max_value<T>() : tools::min_value<T>(); | |
88 | } | |
89 | // | |
90 | // Function object, this is the functor whose root | |
91 | // we have to solve: | |
92 | // | |
93 | beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab); | |
94 | // | |
95 | // Tolerance: full precision. | |
96 | // | |
97 | tools::eps_tolerance<T> tol(policies::digits<T, Policy>()); | |
98 | // | |
99 | // Now figure out a starting guess for what a may be, | |
100 | // we'll start out with a value that'll put p or q | |
101 | // right bang in the middle of their range, the functions | |
102 | // are quite sensitive so we should need too many steps | |
103 | // to bracket the root from there: | |
104 | // | |
105 | T guess = 0; | |
106 | T factor = 5; | |
107 | // | |
108 | // Convert variables to parameters of a negative binomial distribution: | |
109 | // | |
110 | T n = b; | |
111 | T sf = swap_ab ? z : 1-z; | |
112 | T sfc = swap_ab ? 1-z : z; | |
113 | T u = swap_ab ? p : q; | |
114 | T v = swap_ab ? q : p; | |
115 | if(u <= pow(sf, n)) | |
116 | { | |
117 | // | |
118 | // Result is less than 1, negative binomial approximation | |
119 | // is useless.... | |
120 | // | |
121 | if((p < q) != swap_ab) | |
122 | { | |
123 | guess = (std::min)(T(b * 2), T(1)); | |
124 | } | |
125 | else | |
126 | { | |
127 | guess = (std::min)(T(b / 2), T(1)); | |
128 | } | |
129 | } | |
130 | if(n * n * n * u * sf > 0.005) | |
131 | guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol); | |
132 | ||
133 | if(guess < 10) | |
134 | { | |
135 | // | |
136 | // Negative binomial approximation not accurate in this area: | |
137 | // | |
138 | if((p < q) != swap_ab) | |
139 | { | |
140 | guess = (std::min)(T(b * 2), T(10)); | |
141 | } | |
142 | else | |
143 | { | |
144 | guess = (std::min)(T(b / 2), T(10)); | |
145 | } | |
146 | } | |
147 | else | |
148 | factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); | |
149 | BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); | |
150 | // | |
151 | // Max iterations permitted: | |
152 | // | |
153 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
154 | std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol); | |
155 | if(max_iter >= policies::get_max_root_iterations<Policy>()) | |
156 | return policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol); | |
157 | return (r.first + r.second) / 2; | |
158 | } | |
159 | ||
160 | } // namespace detail | |
161 | ||
162 | template <class RT1, class RT2, class RT3, class Policy> | |
163 | typename tools::promote_args<RT1, RT2, RT3>::type | |
164 | ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol) | |
165 | { | |
166 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
167 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
168 | typedef typename policies::normalise< | |
169 | Policy, | |
170 | policies::promote_float<false>, | |
171 | policies::promote_double<false>, | |
172 | policies::discrete_quantile<>, | |
173 | policies::assert_undefined<> >::type forwarding_policy; | |
174 | ||
175 | static const char* function = "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)"; | |
176 | if(p == 0) | |
177 | { | |
178 | return policies::raise_overflow_error<result_type>(function, 0, Policy()); | |
179 | } | |
180 | if(p == 1) | |
181 | { | |
182 | return tools::min_value<result_type>(); | |
183 | } | |
184 | ||
185 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
186 | detail::ibeta_inv_ab_imp( | |
187 | static_cast<value_type>(b), | |
188 | static_cast<value_type>(x), | |
189 | static_cast<value_type>(p), | |
190 | static_cast<value_type>(1 - static_cast<value_type>(p)), | |
191 | false, pol), | |
192 | function); | |
193 | } | |
194 | ||
195 | template <class RT1, class RT2, class RT3, class Policy> | |
196 | typename tools::promote_args<RT1, RT2, RT3>::type | |
197 | ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol) | |
198 | { | |
199 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
200 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
201 | typedef typename policies::normalise< | |
202 | Policy, | |
203 | policies::promote_float<false>, | |
204 | policies::promote_double<false>, | |
205 | policies::discrete_quantile<>, | |
206 | policies::assert_undefined<> >::type forwarding_policy; | |
207 | ||
208 | static const char* function = "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)"; | |
209 | if(q == 1) | |
210 | { | |
211 | return policies::raise_overflow_error<result_type>(function, 0, Policy()); | |
212 | } | |
213 | if(q == 0) | |
214 | { | |
215 | return tools::min_value<result_type>(); | |
216 | } | |
217 | ||
218 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
219 | detail::ibeta_inv_ab_imp( | |
220 | static_cast<value_type>(b), | |
221 | static_cast<value_type>(x), | |
222 | static_cast<value_type>(1 - static_cast<value_type>(q)), | |
223 | static_cast<value_type>(q), | |
224 | false, pol), | |
225 | function); | |
226 | } | |
227 | ||
228 | template <class RT1, class RT2, class RT3, class Policy> | |
229 | typename tools::promote_args<RT1, RT2, RT3>::type | |
230 | ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol) | |
231 | { | |
232 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
233 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
234 | typedef typename policies::normalise< | |
235 | Policy, | |
236 | policies::promote_float<false>, | |
237 | policies::promote_double<false>, | |
238 | policies::discrete_quantile<>, | |
239 | policies::assert_undefined<> >::type forwarding_policy; | |
240 | ||
241 | static const char* function = "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)"; | |
242 | if(p == 0) | |
243 | { | |
244 | return tools::min_value<result_type>(); | |
245 | } | |
246 | if(p == 1) | |
247 | { | |
248 | return policies::raise_overflow_error<result_type>(function, 0, Policy()); | |
249 | } | |
250 | ||
251 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
252 | detail::ibeta_inv_ab_imp( | |
253 | static_cast<value_type>(a), | |
254 | static_cast<value_type>(x), | |
255 | static_cast<value_type>(p), | |
256 | static_cast<value_type>(1 - static_cast<value_type>(p)), | |
257 | true, pol), | |
258 | function); | |
259 | } | |
260 | ||
261 | template <class RT1, class RT2, class RT3, class Policy> | |
262 | typename tools::promote_args<RT1, RT2, RT3>::type | |
263 | ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol) | |
264 | { | |
265 | static const char* function = "boost::math::ibeta_invb<%1%>(%1%, %1%, %1%)"; | |
266 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
267 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
268 | typedef typename policies::normalise< | |
269 | Policy, | |
270 | policies::promote_float<false>, | |
271 | policies::promote_double<false>, | |
272 | policies::discrete_quantile<>, | |
273 | policies::assert_undefined<> >::type forwarding_policy; | |
274 | ||
275 | if(q == 1) | |
276 | { | |
277 | return tools::min_value<result_type>(); | |
278 | } | |
279 | if(q == 0) | |
280 | { | |
281 | return policies::raise_overflow_error<result_type>(function, 0, Policy()); | |
282 | } | |
283 | ||
284 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( | |
285 | detail::ibeta_inv_ab_imp( | |
286 | static_cast<value_type>(a), | |
287 | static_cast<value_type>(x), | |
288 | static_cast<value_type>(1 - static_cast<value_type>(q)), | |
289 | static_cast<value_type>(q), | |
290 | true, pol), | |
291 | function); | |
292 | } | |
293 | ||
294 | template <class RT1, class RT2, class RT3> | |
295 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
296 | ibeta_inva(RT1 b, RT2 x, RT3 p) | |
297 | { | |
298 | return boost::math::ibeta_inva(b, x, p, policies::policy<>()); | |
299 | } | |
300 | ||
301 | template <class RT1, class RT2, class RT3> | |
302 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
303 | ibetac_inva(RT1 b, RT2 x, RT3 q) | |
304 | { | |
305 | return boost::math::ibetac_inva(b, x, q, policies::policy<>()); | |
306 | } | |
307 | ||
308 | template <class RT1, class RT2, class RT3> | |
309 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
310 | ibeta_invb(RT1 a, RT2 x, RT3 p) | |
311 | { | |
312 | return boost::math::ibeta_invb(a, x, p, policies::policy<>()); | |
313 | } | |
314 | ||
315 | template <class RT1, class RT2, class RT3> | |
316 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
317 | ibetac_invb(RT1 a, RT2 x, RT3 q) | |
318 | { | |
319 | return boost::math::ibetac_invb(a, x, q, policies::policy<>()); | |
320 | } | |
321 | ||
322 | } // namespace math | |
323 | } // namespace boost | |
324 | ||
325 | #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB | |
326 | ||
327 | ||
328 |