]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 purturb 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* symetric 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 purturbation 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 eqaution 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 distribtion | |
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 expantions, 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 |