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
19namespace boost{ namespace math{ namespace detail{
20
21//
22// Helper object used by root finding
23// code to convert eta to x.
24//
25template <class T>
26struct 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 }
49private:
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//
59template <class T, class Policy>
60T 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//
137template <class T, class Policy>
138T 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//
314template <class T, class Policy>
315T 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
413template <class T, class Policy>
414struct 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 }
447private:
448 T a, b, target;
449 bool invert;
450};
451
452template <class T, class Policy>
453T 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
870template <class T1, class T2, class T3, class T4, class Policy>
871inline 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
905template <class T1, class T2, class T3, class T4>
906inline 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
912template <class T1, class T2, class T3>
913inline 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
920template <class T1, class T2, class T3, class Policy>
921inline 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
928template <class T1, class T2, class T3, class T4, class Policy>
929inline 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
963template <class T1, class T2, class T3, class T4>
964inline 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
970template <class RT1, class RT2, class RT3>
971inline 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
978template <class RT1, class RT2, class RT3, class Policy>
979inline 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

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