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

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