1 | // Copyright John Maddock 2006. |
2 | // Copyright Paul A. Bristow 2007 |
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_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP |
8 | #define BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP |
9 | |
10 | #ifdef _MSC_VER |
11 | #pragma once |
12 | #endif |
13 | |
14 | #include <boost/math/special_functions/beta.hpp> |
15 | #include <boost/math/special_functions/erf.hpp> |
16 | #include <boost/math/tools/roots.hpp> |
17 | #include <boost/math/special_functions/detail/t_distribution_inv.hpp> |
18 | |
19 | namespace boost{ namespace math{ namespace detail{ |
20 | |
21 | // |
22 | // Helper object used by root finding |
23 | // code to convert eta to x. |
24 | // |
25 | template <class T> |
26 | struct temme_root_finder |
27 | { |
28 | temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {} |
29 | |
30 | boost::math::tuple<T, T> operator()(T x) |
31 | { |
32 | BOOST_MATH_STD_USING // ADL of std names |
33 | |
34 | T y = 1 - x; |
35 | if(y == 0) |
36 | { |
37 | T big = tools::max_value<T>() / 4; |
38 | return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big)); |
39 | } |
40 | if(x == 0) |
41 | { |
42 | T big = tools::max_value<T>() / 4; |
43 | return boost::math::make_tuple(static_cast<T>(-big), big); |
44 | } |
45 | T f = log(x) + a * log(y) + t; |
46 | T f1 = (1 / x) - (a / (y)); |
47 | return boost::math::make_tuple(f, f1); |
48 | } |
49 | private: |
50 | T t, a; |
51 | }; |
52 | // |
53 | // See: |
54 | // "Asymptotic Inversion of the Incomplete Beta Function" |
55 | // N.M. Temme |
56 | // Journal of Computation and Applied Mathematics 41 (1992) 145-157. |
57 | // Section 2. |
58 | // |
59 | template <class T, class Policy> |
60 | T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol) |
61 | { |
62 | BOOST_MATH_STD_USING // ADL of std names |
63 | |
64 | const T r2 = sqrt(T(2)); |
65 | // |
66 | // get the first approximation for eta from the inverse |
67 | // error function (Eq: 2.9 and 2.10). |
68 | // |
69 | T eta0 = boost::math::erfc_inv(2 * z, pol); |
70 | eta0 /= -sqrt(a / 2); |
71 | |
72 | T terms[4] = { eta0 }; |
73 | T workspace[7]; |
74 | // |
75 | // calculate powers: |
76 | // |
77 | T B = b - a; |
78 | T B_2 = B * B; |
79 | T B_3 = B_2 * B; |
80 | // |
81 | // Calculate correction terms: |
82 | // |
83 | |
84 | // See eq following 2.15: |
85 | workspace[0] = -B * r2 / 2; |
86 | workspace[1] = (1 - 2 * B) / 8; |
87 | workspace[2] = -(B * r2 / 48); |
88 | workspace[3] = T(-1) / 192; |
89 | workspace[4] = -B * r2 / 3840; |
90 | terms[1] = tools::evaluate_polynomial(workspace, eta0, 5); |
91 | // Eq Following 2.17: |
92 | workspace[0] = B * r2 * (3 * B - 2) / 12; |
93 | workspace[1] = (20 * B_2 - 12 * B + 1) / 128; |
94 | workspace[2] = B * r2 * (20 * B - 1) / 960; |
95 | workspace[3] = (16 * B_2 + 30 * B - 15) / 4608; |
96 | workspace[4] = B * r2 * (21 * B + 32) / 53760; |
97 | workspace[5] = (-32 * B_2 + 63) / 368640; |
98 | workspace[6] = -B * r2 * (120 * B + 17) / 25804480; |
99 | terms[2] = tools::evaluate_polynomial(workspace, eta0, 7); |
100 | // Eq Following 2.17: |
101 | workspace[0] = B * r2 * (-75 * B_2 + 80 * B - 16) / 480; |
102 | workspace[1] = (-1080 * B_3 + 868 * B_2 - 90 * B - 45) / 9216; |
103 | workspace[2] = B * r2 * (-1190 * B_2 + 84 * B + 373) / 53760; |
104 | workspace[3] = (-2240 * B_3 - 2508 * B_2 + 2100 * B - 165) / 368640; |
105 | terms[3] = tools::evaluate_polynomial(workspace, eta0, 4); |
106 | // |
107 | // Bring them together to get a final estimate for eta: |
108 | // |
109 | T eta = tools::evaluate_polynomial(terms, T(1/a), 4); |
110 | // |
111 | // now we need to convert eta to x, by solving the appropriate |
112 | // quadratic equation: |
113 | // |
114 | T eta_2 = eta * eta; |
115 | T c = -exp(-eta_2 / 2); |
116 | T x; |
117 | if(eta_2 == 0) |
118 | x = 0.5; |
119 | else |
120 | x = (1 + eta * sqrt((1 + c) / eta_2)) / 2; |
121 | |
122 | BOOST_ASSERT(x >= 0); |
123 | BOOST_ASSERT(x <= 1); |
124 | BOOST_ASSERT(eta * (x - 0.5) >= 0); |
125 | #ifdef BOOST_INSTRUMENT |
126 | std::cout << "Estimating x with Temme method 1: " << x << std::endl; |
127 | #endif |
128 | return x; |
129 | } |
130 | // |
131 | // See: |
132 | // "Asymptotic Inversion of the Incomplete Beta Function" |
133 | // N.M. Temme |
134 | // Journal of Computation and Applied Mathematics 41 (1992) 145-157. |
135 | // Section 3. |
136 | // |
137 | template <class T, class Policy> |
138 | T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta, const Policy& pol) |
139 | { |
140 | BOOST_MATH_STD_USING // ADL of std names |
141 | |
142 | // |
143 | // Get first estimate for eta, see Eq 3.9 and 3.10, |
144 | // but note there is a typo in Eq 3.10: |
145 | // |
146 | T eta0 = boost::math::erfc_inv(2 * z, pol); |
147 | eta0 /= -sqrt(r / 2); |
148 | |
149 | T s = sin(theta); |
150 | T c = cos(theta); |
151 | // |
152 | // Now we need to perturb eta0 to get eta, which we do by |
153 | // evaluating the polynomial in 1/r at the bottom of page 151, |
154 | // to do this we first need the error terms e1, e2 e3 |
155 | // which we'll fill into the array "terms". Since these |
156 | // terms are themselves polynomials, we'll need another |
157 | // array "workspace" to calculate those... |
158 | // |
159 | T terms[4] = { eta0 }; |
160 | T workspace[6]; |
161 | // |
162 | // some powers of sin(theta)cos(theta) that we'll need later: |
163 | // |
164 | T sc = s * c; |
165 | T sc_2 = sc * sc; |
166 | T sc_3 = sc_2 * sc; |
167 | T sc_4 = sc_2 * sc_2; |
168 | T sc_5 = sc_2 * sc_3; |
169 | T sc_6 = sc_3 * sc_3; |
170 | T sc_7 = sc_4 * sc_3; |
171 | // |
172 | // Calculate e1 and put it in terms[1], see the middle of page 151: |
173 | // |
174 | workspace[0] = (2 * s * s - 1) / (3 * s * c); |
175 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co1[] = { -1, -5, 5 }; |
176 | workspace[1] = -tools::evaluate_even_polynomial(co1, s, 3) / (36 * sc_2); |
177 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co2[] = { 1, 21, -69, 46 }; |
178 | workspace[2] = tools::evaluate_even_polynomial(co2, s, 4) / (1620 * sc_3); |
179 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co3[] = { 7, -2, 33, -62, 31 }; |
180 | workspace[3] = -tools::evaluate_even_polynomial(co3, s, 5) / (6480 * sc_4); |
181 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co4[] = { 25, -52, -17, 88, -115, 46 }; |
182 | workspace[4] = tools::evaluate_even_polynomial(co4, s, 6) / (90720 * sc_5); |
183 | terms[1] = tools::evaluate_polynomial(workspace, eta0, 5); |
184 | // |
185 | // Now evaluate e2 and put it in terms[2]: |
186 | // |
187 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co5[] = { 7, 12, -78, 52 }; |
188 | workspace[0] = -tools::evaluate_even_polynomial(co5, s, 4) / (405 * sc_3); |
189 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co6[] = { -7, 2, 183, -370, 185 }; |
190 | workspace[1] = tools::evaluate_even_polynomial(co6, s, 5) / (2592 * sc_4); |
191 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co7[] = { -533, 776, -1835, 10240, -13525, 5410 }; |
192 | workspace[2] = -tools::evaluate_even_polynomial(co7, s, 6) / (204120 * sc_5); |
193 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co8[] = { -1579, 3747, -3372, -15821, 45588, -45213, 15071 }; |
194 | workspace[3] = -tools::evaluate_even_polynomial(co8, s, 7) / (2099520 * sc_6); |
195 | terms[2] = tools::evaluate_polynomial(workspace, eta0, 4); |
196 | // |
197 | // And e3, and put it in terms[3]: |
198 | // |
199 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co9[] = {449, -1259, -769, 6686, -9260, 3704 }; |
200 | workspace[0] = tools::evaluate_even_polynomial(co9, s, 6) / (102060 * sc_5); |
201 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co10[] = { 63149, -151557, 140052, -727469, 2239932, -2251437, 750479 }; |
202 | workspace[1] = -tools::evaluate_even_polynomial(co10, s, 7) / (20995200 * sc_6); |
203 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co11[] = { 29233, -78755, 105222, 146879, -1602610, 3195183, -2554139, 729754 }; |
204 | workspace[2] = tools::evaluate_even_polynomial(co11, s, 8) / (36741600 * sc_7); |
205 | terms[3] = tools::evaluate_polynomial(workspace, eta0, 3); |
206 | // |
207 | // Bring the correction terms together to evaluate eta, |
208 | // this is the last equation on page 151: |
209 | // |
210 | T eta = tools::evaluate_polynomial(terms, T(1/r), 4); |
211 | // |
212 | // Now that we have eta we need to back solve for x, |
213 | // we seek the value of x that gives eta in Eq 3.2. |
214 | // The two methods used are described in section 5. |
215 | // |
216 | // Begin by defining a few variables we'll need later: |
217 | // |
218 | T x; |
219 | T s_2 = s * s; |
220 | T c_2 = c * c; |
221 | T alpha = c / s; |
222 | alpha *= alpha; |
223 | T lu = (-(eta * eta) / (2 * s_2) + log(s_2) + c_2 * log(c_2) / s_2); |
224 | // |
225 | // Temme doesn't specify what value to switch on here, |
226 | // but this seems to work pretty well: |
227 | // |
228 | if(fabs(eta) < 0.7) |
229 | { |
230 | // |
231 | // Small eta use the expansion Temme gives in the second equation |
232 | // of section 5, it's a polynomial in eta: |
233 | // |
234 | workspace[0] = s * s; |
235 | workspace[1] = s * c; |
236 | workspace[2] = (1 - 2 * workspace[0]) / 3; |
237 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co12[] = { 1, -13, 13 }; |
238 | workspace[3] = tools::evaluate_polynomial(co12, workspace[0], 3) / (36 * s * c); |
239 | static const BOOST_MATH_INT_TABLE_TYPE(T, int) co13[] = { 1, 21, -69, 46 }; |
240 | workspace[4] = tools::evaluate_polynomial(co13, workspace[0], 4) / (270 * workspace[0] * c * c); |
241 | x = tools::evaluate_polynomial(workspace, eta, 5); |
242 | #ifdef BOOST_INSTRUMENT |
243 | std::cout << "Estimating x with Temme method 2 (small eta): " << x << std::endl; |
244 | #endif |
245 | } |
246 | else |
247 | { |
248 | // |
249 | // If eta is large we need to solve Eq 3.2 more directly, |
250 | // begin by getting an initial approximation for x from |
251 | // the last equation on page 155, this is a polynomial in u: |
252 | // |
253 | T u = exp(lu); |
254 | workspace[0] = u; |
255 | workspace[1] = alpha; |
256 | workspace[2] = 0; |
257 | workspace[3] = 3 * alpha * (3 * alpha + 1) / 6; |
258 | workspace[4] = 4 * alpha * (4 * alpha + 1) * (4 * alpha + 2) / 24; |
259 | workspace[5] = 5 * alpha * (5 * alpha + 1) * (5 * alpha + 2) * (5 * alpha + 3) / 120; |
260 | x = tools::evaluate_polynomial(workspace, u, 6); |
261 | // |
262 | // At this point we may or may not have the right answer, Eq-3.2 has |
263 | // two solutions for x for any given eta, however the mapping in 3.2 |
264 | // is 1:1 with the sign of eta and x-sin^2(theta) being the same. |
265 | // So we can check if we have the right root of 3.2, and if not |
266 | // switch x for 1-x. This transformation is motivated by the fact |
267 | // that the distribution is *almost* symmetric so 1-x will be in the right |
268 | // ball park for the solution: |
269 | // |
270 | if((x - s_2) * eta < 0) |
271 | x = 1 - x; |
272 | #ifdef BOOST_INSTRUMENT |
273 | std::cout << "Estimating x with Temme method 2 (large eta): " << x << std::endl; |
274 | #endif |
275 | } |
276 | // |
277 | // The final step is a few Newton-Raphson iterations to |
278 | // clean up our approximation for x, this is pretty cheap |
279 | // in general, and very cheap compared to an incomplete beta |
280 | // evaluation. The limits set on x come from the observation |
281 | // that the sign of eta and x-sin^2(theta) are the same. |
282 | // |
283 | T lower, upper; |
284 | if(eta < 0) |
285 | { |
286 | lower = 0; |
287 | upper = s_2; |
288 | } |
289 | else |
290 | { |
291 | lower = s_2; |
292 | upper = 1; |
293 | } |
294 | // |
295 | // If our initial approximation is out of bounds then bisect: |
296 | // |
297 | if((x < lower) || (x > upper)) |
298 | x = (lower+upper) / 2; |
299 | // |
300 | // And iterate: |
301 | // |
302 | x = tools::newton_raphson_iterate( |
303 | temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2); |
304 | |
305 | return x; |
306 | } |
307 | // |
308 | // See: |
309 | // "Asymptotic Inversion of the Incomplete Beta Function" |
310 | // N.M. Temme |
311 | // Journal of Computation and Applied Mathematics 41 (1992) 145-157. |
312 | // Section 4. |
313 | // |
314 | template <class T, class Policy> |
315 | T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol) |
316 | { |
317 | BOOST_MATH_STD_USING // ADL of std names |
318 | |
319 | // |
320 | // Begin by getting an initial approximation for the quantity |
321 | // eta from the dominant part of the incomplete beta: |
322 | // |
323 | T eta0; |
324 | if(p < q) |
325 | eta0 = boost::math::gamma_q_inv(b, p, pol); |
326 | else |
327 | eta0 = boost::math::gamma_p_inv(b, q, pol); |
328 | eta0 /= a; |
329 | // |
330 | // Define the variables and powers we'll need later on: |
331 | // |
332 | T mu = b / a; |
333 | T w = sqrt(1 + mu); |
334 | T w_2 = w * w; |
335 | T w_3 = w_2 * w; |
336 | T w_4 = w_2 * w_2; |
337 | T w_5 = w_3 * w_2; |
338 | T w_6 = w_3 * w_3; |
339 | T w_7 = w_4 * w_3; |
340 | T w_8 = w_4 * w_4; |
341 | T w_9 = w_5 * w_4; |
342 | T w_10 = w_5 * w_5; |
343 | T d = eta0 - mu; |
344 | T d_2 = d * d; |
345 | T d_3 = d_2 * d; |
346 | T d_4 = d_2 * d_2; |
347 | T w1 = w + 1; |
348 | T w1_2 = w1 * w1; |
349 | T w1_3 = w1 * w1_2; |
350 | T w1_4 = w1_2 * w1_2; |
351 | // |
352 | // Now we need to compute the perturbation error terms that |
353 | // convert eta0 to eta, these are all polynomials of polynomials. |
354 | // Probably these should be re-written to use tabulated data |
355 | // (see examples above), but it's less of a win in this case as we |
356 | // need to calculate the individual powers for the denominator terms |
357 | // anyway, so we might as well use them for the numerator-polynomials |
358 | // as well.... |
359 | // |
360 | // Refer to p154-p155 for the details of these expansions: |
361 | // |
362 | T e1 = (w + 2) * (w - 1) / (3 * w); |
363 | e1 += (w_3 + 9 * w_2 + 21 * w + 5) * d / (36 * w_2 * w1); |
364 | e1 -= (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2 / (1620 * w1_2 * w_3); |
365 | e1 -= (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3 / (6480 * w1_3 * w_4); |
366 | e1 -= (75 * w_6 + 202 * w_5 + 188 * w_4 - 888 * w_3 - 1345 * w_2 + 118 * w + 138) * d_4 / (272160 * w1_4 * w_5); |
367 | |
368 | T e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) / (1620 * w1 * w_3); |
369 | e2 -= (35 * w_6 - 154 * w_5 - 623 * w_4 - 1636 * w_3 - 3983 * w_2 - 3514 * w - 925) * d / (12960 * w1_2 * w_4); |
370 | e2 -= (2132 * w_7 + 7915 * w_6 + 16821 * w_5 + 35066 * w_4 + 87490 * w_3 + 141183 * w_2 + 95993 * w + 21640) * d_2 / (816480 * w_5 * w1_3); |
371 | e2 -= (11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4 - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497) * d_3 / (14696640 * w1_4 * w_6); |
372 | |
373 | T e3 = -((3592 * w_7 + 8375 * w_6 - 1323 * w_5 - 29198 * w_4 - 89578 * w_3 - 154413 * w_2 - 116063 * w - 29632) * (w - 1)) / (816480 * w_5 * w1_2); |
374 | e3 -= (442043 * w_9 + 2054169 * w_8 + 3803094 * w_7 + 3470754 * w_6 + 2141568 * w_5 - 2393568 * w_4 - 19904934 * w_3 - 34714674 * w_2 - 23128299 * w - 5253353) * d / (146966400 * w_6 * w1_3); |
375 | e3 -= (116932 * w_10 + 819281 * w_9 + 2378172 * w_8 + 4341330 * w_7 + 6806004 * w_6 + 10622748 * w_5 + 18739500 * w_4 + 30651894 * w_3 + 30869976 * w_2 + 15431867 * w + 2919016) * d_2 / (146966400 * w1_4 * w_7); |
376 | // |
377 | // Combine eta0 and the error terms to compute eta (Second equation p155): |
378 | // |
379 | T eta = eta0 + e1 / a + e2 / (a * a) + e3 / (a * a * a); |
380 | // |
381 | // Now we need to solve Eq 4.2 to obtain x. For any given value of |
382 | // eta there are two solutions to this equation, and since the distribution |
383 | // may be very skewed, these are not related by x ~ 1-x we used when |
384 | // implementing section 3 above. However we know that: |
385 | // |
386 | // cross < x <= 1 ; iff eta < mu |
387 | // x == cross ; iff eta == mu |
388 | // 0 <= x < cross ; iff eta > mu |
389 | // |
390 | // Where cross == 1 / (1 + mu) |
391 | // Many thanks to Prof Temme for clarifying this point. |
392 | // |
393 | // Therefore we'll just jump straight into Newton iterations |
394 | // to solve Eq 4.2 using these bounds, and simple bisection |
395 | // as the first guess, in practice this converges pretty quickly |
396 | // and we only need a few digits correct anyway: |
397 | // |
398 | if(eta <= 0) |
399 | eta = tools::min_value<T>(); |
400 | T u = eta - mu * log(eta) + (1 + mu) * log(1 + mu) - mu; |
401 | T cross = 1 / (1 + mu); |
402 | T lower = eta < mu ? cross : 0; |
403 | T upper = eta < mu ? 1 : cross; |
404 | T x = (lower + upper) / 2; |
405 | x = tools::newton_raphson_iterate( |
406 | temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2); |
407 | #ifdef BOOST_INSTRUMENT |
408 | std::cout << "Estimating x with Temme method 3: " << x << std::endl; |
409 | #endif |
410 | return x; |
411 | } |
412 | |
413 | template <class T, class Policy> |
414 | struct ibeta_roots |
415 | { |
416 | ibeta_roots(T _a, T _b, T t, bool inv = false) |
417 | : a(_a), b(_b), target(t), invert(inv) {} |
418 | |
419 | boost::math::tuple<T, T, T> operator()(T x) |
420 | { |
421 | BOOST_MATH_STD_USING // ADL of std names |
422 | |
423 | BOOST_FPU_EXCEPTION_GUARD |
424 | |
425 | T f1; |
426 | T y = 1 - x; |
427 | T f = ibeta_imp(a, b, x, Policy(), invert, true, &f1) - target; |
428 | if(invert) |
429 | f1 = -f1; |
430 | if(y == 0) |
431 | y = tools::min_value<T>() * 64; |
432 | if(x == 0) |
433 | x = tools::min_value<T>() * 64; |
434 | |
435 | T f2 = f1 * (-y * a + (b - 2) * x + 1); |
436 | if(fabs(f2) < y * x * tools::max_value<T>()) |
437 | f2 /= (y * x); |
438 | if(invert) |
439 | f2 = -f2; |
440 | |
441 | // make sure we don't have a zero derivative: |
442 | if(f1 == 0) |
443 | f1 = (invert ? -1 : 1) * tools::min_value<T>() * 64; |
444 | |
445 | return boost::math::make_tuple(f, f1, f2); |
446 | } |
447 | private: |
448 | T a, b, target; |
449 | bool invert; |
450 | }; |
451 | |
452 | template <class T, class Policy> |
453 | T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) |
454 | { |
455 | BOOST_MATH_STD_USING // For ADL of math functions. |
456 | |
457 | // |
458 | // The flag invert is set to true if we swap a for b and p for q, |
459 | // in which case the result has to be subtracted from 1: |
460 | // |
461 | bool invert = false; |
462 | // |
463 | // Handle trivial cases first: |
464 | // |
465 | if(q == 0) |
466 | { |
467 | if(py) *py = 0; |
468 | return 1; |
469 | } |
470 | else if(p == 0) |
471 | { |
472 | if(py) *py = 1; |
473 | return 0; |
474 | } |
475 | else if(a == 1) |
476 | { |
477 | if(b == 1) |
478 | { |
479 | if(py) *py = 1 - p; |
480 | return p; |
481 | } |
482 | // Change things around so we can handle as b == 1 special case below: |
483 | std::swap(a, b); |
484 | std::swap(p, q); |
485 | invert = true; |
486 | } |
487 | // |
488 | // Depending upon which approximation method we use, we may end up |
489 | // calculating either x or y initially (where y = 1-x): |
490 | // |
491 | T x = 0; // Set to a safe zero to avoid a |
492 | // MSVC 2005 warning C4701: potentially uninitialized local variable 'x' used |
493 | // But code inspection appears to ensure that x IS assigned whatever the code path. |
494 | T y; |
495 | |
496 | // For some of the methods we can put tighter bounds |
497 | // on the result than simply [0,1]: |
498 | // |
499 | T lower = 0; |
500 | T upper = 1; |
501 | // |
502 | // Student's T with b = 0.5 gets handled as a special case, swap |
503 | // around if the arguments are in the "wrong" order: |
504 | // |
505 | if(a == 0.5f) |
506 | { |
507 | if(b == 0.5f) |
508 | { |
509 | x = sin(p * constants::half_pi<T>()); |
510 | x *= x; |
511 | if(py) |
512 | { |
513 | *py = sin(q * constants::half_pi<T>()); |
514 | *py *= *py; |
515 | } |
516 | return x; |
517 | } |
518 | else if(b > 0.5f) |
519 | { |
520 | std::swap(a, b); |
521 | std::swap(p, q); |
522 | invert = !invert; |
523 | } |
524 | } |
525 | // |
526 | // Select calculation method for the initial estimate: |
527 | // |
528 | if((b == 0.5f) && (a >= 0.5f) && (p != 1)) |
529 | { |
530 | // |
531 | // We have a Student's T distribution: |
532 | x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol); |
533 | } |
534 | else if(b == 1) |
535 | { |
536 | if(p < q) |
537 | { |
538 | if(a > 1) |
539 | { |
540 | x = pow(p, 1 / a); |
541 | y = -boost::math::expm1(log(p) / a, pol); |
542 | } |
543 | else |
544 | { |
545 | x = pow(p, 1 / a); |
546 | y = 1 - x; |
547 | } |
548 | } |
549 | else |
550 | { |
551 | x = exp(boost::math::log1p(-q, pol) / a); |
552 | y = -boost::math::expm1(boost::math::log1p(-q, pol) / a, pol); |
553 | } |
554 | if(invert) |
555 | std::swap(x, y); |
556 | if(py) |
557 | *py = y; |
558 | return x; |
559 | } |
560 | else if(a + b > 5) |
561 | { |
562 | // |
563 | // When a+b is large then we can use one of Prof Temme's |
564 | // asymptotic expansions, begin by swapping things around |
565 | // so that p < 0.5, we do this to avoid cancellations errors |
566 | // when p is large. |
567 | // |
568 | if(p > 0.5) |
569 | { |
570 | std::swap(a, b); |
571 | std::swap(p, q); |
572 | invert = !invert; |
573 | } |
574 | T minv = (std::min)(a, b); |
575 | T maxv = (std::max)(a, b); |
576 | if((sqrt(minv) > (maxv - minv)) && (minv > 5)) |
577 | { |
578 | // |
579 | // When a and b differ by a small amount |
580 | // the curve is quite symmetrical and we can use an error |
581 | // function to approximate the inverse. This is the cheapest |
582 | // of the three Temme expansions, and the calculated value |
583 | // for x will never be much larger than p, so we don't have |
584 | // to worry about cancellation as long as p is small. |
585 | // |
586 | x = temme_method_1_ibeta_inverse(a, b, p, pol); |
587 | y = 1 - x; |
588 | } |
589 | else |
590 | { |
591 | T r = a + b; |
592 | T theta = asin(sqrt(a / r)); |
593 | T lambda = minv / r; |
594 | if((lambda >= 0.2) && (lambda <= 0.8) && (r >= 10)) |
595 | { |
596 | // |
597 | // The second error function case is the next cheapest |
598 | // to use, it brakes down when the result is likely to be |
599 | // very small, if a+b is also small, but we can use a |
600 | // cheaper expansion there in any case. As before x won't |
601 | // be much larger than p, so as long as p is small we should |
602 | // be free of cancellation error. |
603 | // |
604 | T ppa = pow(p, 1/a); |
605 | if((ppa < 0.0025) && (a + b < 200)) |
606 | { |
607 | x = ppa * pow(a * boost::math::beta(a, b, pol), 1/a); |
608 | } |
609 | else |
610 | x = temme_method_2_ibeta_inverse(a, b, p, r, theta, pol); |
611 | y = 1 - x; |
612 | } |
613 | else |
614 | { |
615 | // |
616 | // If we get here then a and b are very different in magnitude |
617 | // and we need to use the third of Temme's methods which |
618 | // involves inverting the incomplete gamma. This is much more |
619 | // expensive than the other methods. We also can only use this |
620 | // method when a > b, which can lead to cancellation errors |
621 | // if we really want y (as we will when x is close to 1), so |
622 | // a different expansion is used in that case. |
623 | // |
624 | if(a < b) |
625 | { |
626 | std::swap(a, b); |
627 | std::swap(p, q); |
628 | invert = !invert; |
629 | } |
630 | // |
631 | // Try and compute the easy way first: |
632 | // |
633 | T bet = 0; |
634 | if(b < 2) |
635 | bet = boost::math::beta(a, b, pol); |
636 | if(bet != 0) |
637 | { |
638 | y = pow(b * q * bet, 1/b); |
639 | x = 1 - y; |
640 | } |
641 | else |
642 | y = 1; |
643 | if(y > 1e-5) |
644 | { |
645 | x = temme_method_3_ibeta_inverse(a, b, p, q, pol); |
646 | y = 1 - x; |
647 | } |
648 | } |
649 | } |
650 | } |
651 | else if((a < 1) && (b < 1)) |
652 | { |
653 | // |
654 | // Both a and b less than 1, |
655 | // there is a point of inflection at xs: |
656 | // |
657 | T xs = (1 - a) / (2 - a - b); |
658 | // |
659 | // Now we need to ensure that we start our iteration from the |
660 | // right side of the inflection point: |
661 | // |
662 | T fs = boost::math::ibeta(a, b, xs, pol) - p; |
663 | if(fabs(fs) / p < tools::epsilon<T>() * 3) |
664 | { |
665 | // The result is at the point of inflection, best just return it: |
666 | *py = invert ? xs : 1 - xs; |
667 | return invert ? 1-xs : xs; |
668 | } |
669 | if(fs < 0) |
670 | { |
671 | std::swap(a, b); |
672 | std::swap(p, q); |
673 | invert = !invert; |
674 | xs = 1 - xs; |
675 | } |
676 | T xg = pow(a * p * boost::math::beta(a, b, pol), 1/a); |
677 | x = xg / (1 + xg); |
678 | y = 1 / (1 + xg); |
679 | // |
680 | // And finally we know that our result is below the inflection |
681 | // point, so set an upper limit on our search: |
682 | // |
683 | if(x > xs) |
684 | x = xs; |
685 | upper = xs; |
686 | } |
687 | else if((a > 1) && (b > 1)) |
688 | { |
689 | // |
690 | // Small a and b, both greater than 1, |
691 | // there is a point of inflection at xs, |
692 | // and it's complement is xs2, we must always |
693 | // start our iteration from the right side of the |
694 | // point of inflection. |
695 | // |
696 | T xs = (a - 1) / (a + b - 2); |
697 | T xs2 = (b - 1) / (a + b - 2); |
698 | T ps = boost::math::ibeta(a, b, xs, pol) - p; |
699 | |
700 | if(ps < 0) |
701 | { |
702 | std::swap(a, b); |
703 | std::swap(p, q); |
704 | std::swap(xs, xs2); |
705 | invert = !invert; |
706 | } |
707 | // |
708 | // Estimate x and y, using expm1 to get a good estimate |
709 | // for y when it's very small: |
710 | // |
711 | T lx = log(p * a * boost::math::beta(a, b, pol)) / a; |
712 | x = exp(lx); |
713 | y = x < 0.9 ? T(1 - x) : (T)(-boost::math::expm1(lx, pol)); |
714 | |
715 | if((b < a) && (x < 0.2)) |
716 | { |
717 | // |
718 | // Under a limited range of circumstances we can improve |
719 | // our estimate for x, frankly it's clear if this has much effect! |
720 | // |
721 | T ap1 = a - 1; |
722 | T bm1 = b - 1; |
723 | T a_2 = a * a; |
724 | T a_3 = a * a_2; |
725 | T b_2 = b * b; |
726 | T terms[5] = { 0, 1 }; |
727 | terms[2] = bm1 / ap1; |
728 | ap1 *= ap1; |
729 | terms[3] = bm1 * (3 * a * b + 5 * b + a_2 - a - 4) / (2 * (a + 2) * ap1); |
730 | ap1 *= (a + 1); |
731 | terms[4] = bm1 * (33 * a * b_2 + 31 * b_2 + 8 * a_2 * b_2 - 30 * a * b - 47 * b + 11 * a_2 * b + 6 * a_3 * b + 18 + 4 * a - a_3 + a_2 * a_2 - 10 * a_2) |
732 | / (3 * (a + 3) * (a + 2) * ap1); |
733 | x = tools::evaluate_polynomial(terms, x, 5); |
734 | } |
735 | // |
736 | // And finally we know that our result is below the inflection |
737 | // point, so set an upper limit on our search: |
738 | // |
739 | if(x > xs) |
740 | x = xs; |
741 | upper = xs; |
742 | } |
743 | else /*if((a <= 1) != (b <= 1))*/ |
744 | { |
745 | // |
746 | // If all else fails we get here, only one of a and b |
747 | // is above 1, and a+b is small. Start by swapping |
748 | // things around so that we have a concave curve with b > a |
749 | // and no points of inflection in [0,1]. As long as we expect |
750 | // x to be small then we can use the simple (and cheap) power |
751 | // term to estimate x, but when we expect x to be large then |
752 | // this greatly underestimates x and leaves us trying to |
753 | // iterate "round the corner" which may take almost forever... |
754 | // |
755 | // We could use Temme's inverse gamma function case in that case, |
756 | // this works really rather well (albeit expensively) even though |
757 | // strictly speaking we're outside it's defined range. |
758 | // |
759 | // However it's expensive to compute, and an alternative approach |
760 | // which models the curve as a distorted quarter circle is much |
761 | // cheaper to compute, and still keeps the number of iterations |
762 | // required down to a reasonable level. With thanks to Prof Temme |
763 | // for this suggestion. |
764 | // |
765 | if(b < a) |
766 | { |
767 | std::swap(a, b); |
768 | std::swap(p, q); |
769 | invert = !invert; |
770 | } |
771 | if(pow(p, 1/a) < 0.5) |
772 | { |
773 | x = pow(p * a * boost::math::beta(a, b, pol), 1 / a); |
774 | if(x == 0) |
775 | x = boost::math::tools::min_value<T>(); |
776 | y = 1 - x; |
777 | } |
778 | else /*if(pow(q, 1/b) < 0.1)*/ |
779 | { |
780 | // model a distorted quarter circle: |
781 | y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b); |
782 | if(y == 0) |
783 | y = boost::math::tools::min_value<T>(); |
784 | x = 1 - y; |
785 | } |
786 | } |
787 | |
788 | // |
789 | // Now we have a guess for x (and for y) we can set things up for |
790 | // iteration. If x > 0.5 it pays to swap things round: |
791 | // |
792 | if(x > 0.5) |
793 | { |
794 | std::swap(a, b); |
795 | std::swap(p, q); |
796 | std::swap(x, y); |
797 | invert = !invert; |
798 | T l = 1 - upper; |
799 | T u = 1 - lower; |
800 | lower = l; |
801 | upper = u; |
802 | } |
803 | // |
804 | // lower bound for our search: |
805 | // |
806 | // We're not interested in denormalised answers as these tend to |
807 | // these tend to take up lots of iterations, given that we can't get |
808 | // accurate derivatives in this area (they tend to be infinite). |
809 | // |
810 | if(lower == 0) |
811 | { |
812 | if(invert && (py == 0)) |
813 | { |
814 | // |
815 | // We're not interested in answers smaller than machine epsilon: |
816 | // |
817 | lower = boost::math::tools::epsilon<T>(); |
818 | if(x < lower) |
819 | x = lower; |
820 | } |
821 | else |
822 | lower = boost::math::tools::min_value<T>(); |
823 | if(x < lower) |
824 | x = lower; |
825 | } |
826 | // |
827 | // Figure out how many digits to iterate towards: |
828 | // |
829 | int digits = boost::math::policies::digits<T, Policy>() / 2; |
830 | if((x < 1e-50) && ((a < 1) || (b < 1))) |
831 | { |
832 | // |
833 | // If we're in a region where the first derivative is very |
834 | // large, then we have to take care that the root-finder |
835 | // doesn't terminate prematurely. We'll bump the precision |
836 | // up to avoid this, but we have to take care not to set the |
837 | // precision too high or the last few iterations will just |
838 | // thrash around and convergence may be slow in this case. |
839 | // Try 3/4 of machine epsilon: |
840 | // |
841 | digits *= 3; |
842 | digits /= 2; |
843 | } |
844 | // |
845 | // Now iterate, we can use either p or q as the target here |
846 | // depending on which is smaller: |
847 | // |
848 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); |
849 | x = boost::math::tools::halley_iterate( |
850 | boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter); |
851 | policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)" , max_iter, pol); |
852 | // |
853 | // We don't really want these asserts here, but they are useful for sanity |
854 | // checking that we have the limits right, uncomment if you suspect bugs *only*. |
855 | // |
856 | //BOOST_ASSERT(x != upper); |
857 | //BOOST_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>())); |
858 | // |
859 | // Tidy up, if we "lower" was too high then zero is the best answer we have: |
860 | // |
861 | if(x == lower) |
862 | x = 0; |
863 | if(py) |
864 | *py = invert ? x : 1 - x; |
865 | return invert ? 1-x : x; |
866 | } |
867 | |
868 | } // namespace detail |
869 | |
870 | template <class T1, class T2, class T3, class T4, class Policy> |
871 | inline typename tools::promote_args<T1, T2, T3, T4>::type |
872 | ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy& pol) |
873 | { |
874 | static const char* function = "boost::math::ibeta_inv<%1%>(%1%,%1%,%1%)" ; |
875 | BOOST_FPU_EXCEPTION_GUARD |
876 | typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type; |
877 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
878 | typedef typename policies::normalise< |
879 | Policy, |
880 | policies::promote_float<false>, |
881 | policies::promote_double<false>, |
882 | policies::discrete_quantile<>, |
883 | policies::assert_undefined<> >::type forwarding_policy; |
884 | |
885 | if(a <= 0) |
886 | return policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%)." , a, pol); |
887 | if(b <= 0) |
888 | return policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%)." , b, pol); |
889 | if((p < 0) || (p > 1)) |
890 | return policies::raise_domain_error<result_type>(function, "Argument p outside the range [0,1] in the incomplete beta function inverse (got p=%1%)." , p, pol); |
891 | |
892 | value_type rx, ry; |
893 | |
894 | rx = detail::ibeta_inv_imp( |
895 | static_cast<value_type>(a), |
896 | static_cast<value_type>(b), |
897 | static_cast<value_type>(p), |
898 | static_cast<value_type>(1 - p), |
899 | forwarding_policy(), &ry); |
900 | |
901 | if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function); |
902 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function); |
903 | } |
904 | |
905 | template <class T1, class T2, class T3, class T4> |
906 | inline typename tools::promote_args<T1, T2, T3, T4>::type |
907 | ibeta_inv(T1 a, T2 b, T3 p, T4* py) |
908 | { |
909 | return ibeta_inv(a, b, p, py, policies::policy<>()); |
910 | } |
911 | |
912 | template <class T1, class T2, class T3> |
913 | inline typename tools::promote_args<T1, T2, T3>::type |
914 | ibeta_inv(T1 a, T2 b, T3 p) |
915 | { |
916 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; |
917 | return ibeta_inv(a, b, p, static_cast<result_type*>(0), policies::policy<>()); |
918 | } |
919 | |
920 | template <class T1, class T2, class T3, class Policy> |
921 | inline typename tools::promote_args<T1, T2, T3>::type |
922 | ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol) |
923 | { |
924 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; |
925 | return ibeta_inv(a, b, p, static_cast<result_type*>(0), pol); |
926 | } |
927 | |
928 | template <class T1, class T2, class T3, class T4, class Policy> |
929 | inline typename tools::promote_args<T1, T2, T3, T4>::type |
930 | ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy& pol) |
931 | { |
932 | static const char* function = "boost::math::ibetac_inv<%1%>(%1%,%1%,%1%)" ; |
933 | BOOST_FPU_EXCEPTION_GUARD |
934 | typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type; |
935 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
936 | typedef typename policies::normalise< |
937 | Policy, |
938 | policies::promote_float<false>, |
939 | policies::promote_double<false>, |
940 | policies::discrete_quantile<>, |
941 | policies::assert_undefined<> >::type forwarding_policy; |
942 | |
943 | if(a <= 0) |
944 | return policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%)." , a, pol); |
945 | if(b <= 0) |
946 | return policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%)." , b, pol); |
947 | if((q < 0) || (q > 1)) |
948 | return policies::raise_domain_error<result_type>(function, "Argument q outside the range [0,1] in the incomplete beta function inverse (got q=%1%)." , q, pol); |
949 | |
950 | value_type rx, ry; |
951 | |
952 | rx = detail::ibeta_inv_imp( |
953 | static_cast<value_type>(a), |
954 | static_cast<value_type>(b), |
955 | static_cast<value_type>(1 - q), |
956 | static_cast<value_type>(q), |
957 | forwarding_policy(), &ry); |
958 | |
959 | if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function); |
960 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function); |
961 | } |
962 | |
963 | template <class T1, class T2, class T3, class T4> |
964 | inline typename tools::promote_args<T1, T2, T3, T4>::type |
965 | ibetac_inv(T1 a, T2 b, T3 q, T4* py) |
966 | { |
967 | return ibetac_inv(a, b, q, py, policies::policy<>()); |
968 | } |
969 | |
970 | template <class RT1, class RT2, class RT3> |
971 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
972 | ibetac_inv(RT1 a, RT2 b, RT3 q) |
973 | { |
974 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; |
975 | return ibetac_inv(a, b, q, static_cast<result_type*>(0), policies::policy<>()); |
976 | } |
977 | |
978 | template <class RT1, class RT2, class RT3, class Policy> |
979 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
980 | ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol) |
981 | { |
982 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; |
983 | return ibetac_inv(a, b, q, static_cast<result_type*>(0), pol); |
984 | } |
985 | |
986 | } // namespace math |
987 | } // namespace boost |
988 | |
989 | #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP |
990 | |
991 | |
992 | |
993 | |
994 | |