]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SPECIAL_BETA_HPP | |
7 | #define BOOST_MATH_SPECIAL_BETA_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/tools/config.hpp> | |
15 | #include <boost/math/special_functions/gamma.hpp> | |
16 | #include <boost/math/special_functions/binomial.hpp> | |
17 | #include <boost/math/special_functions/factorials.hpp> | |
18 | #include <boost/math/special_functions/erf.hpp> | |
19 | #include <boost/math/special_functions/log1p.hpp> | |
20 | #include <boost/math/special_functions/expm1.hpp> | |
21 | #include <boost/math/special_functions/trunc.hpp> | |
22 | #include <boost/math/tools/roots.hpp> | |
23 | #include <boost/static_assert.hpp> | |
24 | #include <boost/config/no_tr1/cmath.hpp> | |
25 | ||
26 | namespace boost{ namespace math{ | |
27 | ||
28 | namespace detail{ | |
29 | ||
30 | // | |
31 | // Implementation of Beta(a,b) using the Lanczos approximation: | |
32 | // | |
33 | template <class T, class Lanczos, class Policy> | |
34 | T beta_imp(T a, T b, const Lanczos&, const Policy& pol) | |
35 | { | |
36 | BOOST_MATH_STD_USING // for ADL of std names | |
37 | ||
38 | if(a <= 0) | |
39 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); | |
40 | if(b <= 0) | |
41 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); | |
42 | ||
43 | T result; | |
44 | ||
45 | T prefix = 1; | |
46 | T c = a + b; | |
47 | ||
48 | // Special cases: | |
49 | if((c == a) && (b < tools::epsilon<T>())) | |
50 | return 1 / b; | |
51 | else if((c == b) && (a < tools::epsilon<T>())) | |
52 | return 1 / a; | |
53 | if(b == 1) | |
54 | return 1/a; | |
55 | else if(a == 1) | |
56 | return 1/b; | |
57 | else if(c < tools::epsilon<T>()) | |
58 | { | |
59 | result = c / a; | |
60 | result /= b; | |
61 | return result; | |
62 | } | |
63 | ||
64 | /* | |
65 | // | |
66 | // This code appears to be no longer necessary: it was | |
67 | // used to offset errors introduced from the Lanczos | |
68 | // approximation, but the current Lanczos approximations | |
69 | // are sufficiently accurate for all z that we can ditch | |
70 | // this. It remains in the file for future reference... | |
71 | // | |
72 | // If a or b are less than 1, shift to greater than 1: | |
73 | if(a < 1) | |
74 | { | |
75 | prefix *= c / a; | |
76 | c += 1; | |
77 | a += 1; | |
78 | } | |
79 | if(b < 1) | |
80 | { | |
81 | prefix *= c / b; | |
82 | c += 1; | |
83 | b += 1; | |
84 | } | |
85 | */ | |
86 | ||
87 | if(a < b) | |
88 | std::swap(a, b); | |
89 | ||
90 | // Lanczos calculation: | |
91 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
92 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
93 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
94 | result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c)); | |
95 | T ambh = a - 0.5f - b; | |
96 | if((fabs(b * ambh) < (cgh * 100)) && (a > 100)) | |
97 | { | |
98 | // Special case where the base of the power term is close to 1 | |
99 | // compute (1+x)^y instead: | |
100 | result *= exp(ambh * boost::math::log1p(-b / cgh, pol)); | |
101 | } | |
102 | else | |
103 | { | |
104 | result *= pow(agh / cgh, a - T(0.5) - b); | |
105 | } | |
106 | if(cgh > 1e10f) | |
107 | // this avoids possible overflow, but appears to be marginally less accurate: | |
108 | result *= pow((agh / cgh) * (bgh / cgh), b); | |
109 | else | |
110 | result *= pow((agh * bgh) / (cgh * cgh), b); | |
111 | result *= sqrt(boost::math::constants::e<T>() / bgh); | |
112 | ||
113 | // If a and b were originally less than 1 we need to scale the result: | |
114 | result *= prefix; | |
115 | ||
116 | return result; | |
117 | } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&) | |
118 | ||
119 | // | |
120 | // Generic implementation of Beta(a,b) without Lanczos approximation support | |
121 | // (Caution this is slow!!!): | |
122 | // | |
123 | template <class T, class Policy> | |
124 | T beta_imp(T a, T b, const lanczos::undefined_lanczos& /* l */, const Policy& pol) | |
125 | { | |
126 | BOOST_MATH_STD_USING | |
127 | ||
128 | if(a <= 0) | |
129 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); | |
130 | if(b <= 0) | |
131 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); | |
132 | ||
133 | T result; | |
134 | ||
135 | T prefix = 1; | |
136 | T c = a + b; | |
137 | ||
138 | // special cases: | |
139 | if((c == a) && (b < tools::epsilon<T>())) | |
140 | return boost::math::tgamma(b, pol); | |
141 | else if((c == b) && (a < tools::epsilon<T>())) | |
142 | return boost::math::tgamma(a, pol); | |
143 | if(b == 1) | |
144 | return 1/a; | |
145 | else if(a == 1) | |
146 | return 1/b; | |
147 | ||
148 | // shift to a and b > 1 if required: | |
149 | if(a < 1) | |
150 | { | |
151 | prefix *= c / a; | |
152 | c += 1; | |
153 | a += 1; | |
154 | } | |
155 | if(b < 1) | |
156 | { | |
157 | prefix *= c / b; | |
158 | c += 1; | |
159 | b += 1; | |
160 | } | |
161 | if(a < b) | |
162 | std::swap(a, b); | |
163 | ||
164 | // set integration limits: | |
165 | T la = (std::max)(T(10), a); | |
166 | T lb = (std::max)(T(10), b); | |
167 | T lc = (std::max)(T(10), T(a+b)); | |
168 | ||
169 | // calculate the fraction parts: | |
170 | T sa = detail::lower_gamma_series(a, la, pol) / a; | |
171 | sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>()); | |
172 | T sb = detail::lower_gamma_series(b, lb, pol) / b; | |
173 | sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>()); | |
174 | T sc = detail::lower_gamma_series(c, lc, pol) / c; | |
175 | sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>()); | |
176 | ||
177 | // and the exponent part: | |
178 | result = exp(lc - la - lb) * pow(la/lc, a) * pow(lb/lc, b); | |
179 | ||
180 | // and combine: | |
181 | result *= sa * sb / sc; | |
182 | ||
183 | // if a and b were originally less than 1 we need to scale the result: | |
184 | result *= prefix; | |
185 | ||
186 | return result; | |
187 | } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l) | |
188 | ||
189 | ||
190 | // | |
191 | // Compute the leading power terms in the incomplete Beta: | |
192 | // | |
193 | // (x^a)(y^b)/Beta(a,b) when normalised, and | |
194 | // (x^a)(y^b) otherwise. | |
195 | // | |
196 | // Almost all of the error in the incomplete beta comes from this | |
197 | // function: particularly when a and b are large. Computing large | |
198 | // powers are *hard* though, and using logarithms just leads to | |
199 | // horrendous cancellation errors. | |
200 | // | |
201 | template <class T, class Lanczos, class Policy> | |
202 | T ibeta_power_terms(T a, | |
203 | T b, | |
204 | T x, | |
205 | T y, | |
206 | const Lanczos&, | |
207 | bool normalised, | |
208 | const Policy& pol, | |
209 | T prefix = 1, | |
210 | const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") | |
211 | { | |
212 | BOOST_MATH_STD_USING | |
213 | ||
214 | if(!normalised) | |
215 | { | |
216 | // can we do better here? | |
217 | return pow(x, a) * pow(y, b); | |
218 | } | |
219 | ||
220 | T result; | |
221 | ||
222 | T c = a + b; | |
223 | ||
224 | // combine power terms with Lanczos approximation: | |
225 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
226 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
227 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
228 | result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b)); | |
229 | result *= prefix; | |
230 | // combine with the leftover terms from the Lanczos approximation: | |
231 | result *= sqrt(bgh / boost::math::constants::e<T>()); | |
232 | result *= sqrt(agh / cgh); | |
233 | ||
234 | // l1 and l2 are the base of the exponents minus one: | |
235 | T l1 = (x * b - y * agh) / agh; | |
236 | T l2 = (y * a - x * bgh) / bgh; | |
237 | if(((std::min)(fabs(l1), fabs(l2)) < 0.2)) | |
238 | { | |
239 | // when the base of the exponent is very near 1 we get really | |
240 | // gross errors unless extra care is taken: | |
241 | if((l1 * l2 > 0) || ((std::min)(a, b) < 1)) | |
242 | { | |
243 | // | |
244 | // This first branch handles the simple cases where either: | |
245 | // | |
246 | // * The two power terms both go in the same direction | |
247 | // (towards zero or towards infinity). In this case if either | |
248 | // term overflows or underflows, then the product of the two must | |
249 | // do so also. | |
250 | // *Alternatively if one exponent is less than one, then we | |
251 | // can't productively use it to eliminate overflow or underflow | |
252 | // from the other term. Problems with spurious overflow/underflow | |
253 | // can't be ruled out in this case, but it is *very* unlikely | |
254 | // since one of the power terms will evaluate to a number close to 1. | |
255 | // | |
256 | if(fabs(l1) < 0.1) | |
257 | { | |
258 | result *= exp(a * boost::math::log1p(l1, pol)); | |
259 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
260 | } | |
261 | else | |
262 | { | |
263 | result *= pow((x * cgh) / agh, a); | |
264 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
265 | } | |
266 | if(fabs(l2) < 0.1) | |
267 | { | |
268 | result *= exp(b * boost::math::log1p(l2, pol)); | |
269 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
270 | } | |
271 | else | |
272 | { | |
273 | result *= pow((y * cgh) / bgh, b); | |
274 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
275 | } | |
276 | } | |
277 | else if((std::max)(fabs(l1), fabs(l2)) < 0.5) | |
278 | { | |
279 | // | |
280 | // Both exponents are near one and both the exponents are | |
281 | // greater than one and further these two | |
282 | // power terms tend in opposite directions (one towards zero, | |
283 | // the other towards infinity), so we have to combine the terms | |
284 | // to avoid any risk of overflow or underflow. | |
285 | // | |
286 | // We do this by moving one power term inside the other, we have: | |
287 | // | |
288 | // (1 + l1)^a * (1 + l2)^b | |
289 | // = ((1 + l1)*(1 + l2)^(b/a))^a | |
290 | // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1 | |
291 | // = exp((b/a) * log(1 + l2)) - 1 | |
292 | // | |
293 | // The tricky bit is deciding which term to move inside :-) | |
294 | // By preference we move the larger term inside, so that the | |
295 | // size of the largest exponent is reduced. However, that can | |
296 | // only be done as long as l3 (see above) is also small. | |
297 | // | |
298 | bool small_a = a < b; | |
299 | T ratio = b / a; | |
300 | if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1))) | |
301 | { | |
302 | T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol); | |
303 | l3 = l1 + l3 + l3 * l1; | |
304 | l3 = a * boost::math::log1p(l3, pol); | |
305 | result *= exp(l3); | |
306 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
307 | } | |
308 | else | |
309 | { | |
310 | T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol); | |
311 | l3 = l2 + l3 + l3 * l2; | |
312 | l3 = b * boost::math::log1p(l3, pol); | |
313 | result *= exp(l3); | |
314 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
315 | } | |
316 | } | |
317 | else if(fabs(l1) < fabs(l2)) | |
318 | { | |
319 | // First base near 1 only: | |
320 | T l = a * boost::math::log1p(l1, pol) | |
321 | + b * log((y * cgh) / bgh); | |
322 | if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) | |
323 | { | |
324 | l += log(result); | |
325 | if(l >= tools::log_max_value<T>()) | |
326 | return policies::raise_overflow_error<T>(function, 0, pol); | |
327 | result = exp(l); | |
328 | } | |
329 | else | |
330 | result *= exp(l); | |
331 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
332 | } | |
333 | else | |
334 | { | |
335 | // Second base near 1 only: | |
336 | T l = b * boost::math::log1p(l2, pol) | |
337 | + a * log((x * cgh) / agh); | |
338 | if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) | |
339 | { | |
340 | l += log(result); | |
341 | if(l >= tools::log_max_value<T>()) | |
342 | return policies::raise_overflow_error<T>(function, 0, pol); | |
343 | result = exp(l); | |
344 | } | |
345 | else | |
346 | result *= exp(l); | |
347 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
348 | } | |
349 | } | |
350 | else | |
351 | { | |
352 | // general case: | |
353 | T b1 = (x * cgh) / agh; | |
354 | T b2 = (y * cgh) / bgh; | |
355 | l1 = a * log(b1); | |
356 | l2 = b * log(b2); | |
357 | BOOST_MATH_INSTRUMENT_VARIABLE(b1); | |
358 | BOOST_MATH_INSTRUMENT_VARIABLE(b2); | |
359 | BOOST_MATH_INSTRUMENT_VARIABLE(l1); | |
360 | BOOST_MATH_INSTRUMENT_VARIABLE(l2); | |
361 | if((l1 >= tools::log_max_value<T>()) | |
362 | || (l1 <= tools::log_min_value<T>()) | |
363 | || (l2 >= tools::log_max_value<T>()) | |
364 | || (l2 <= tools::log_min_value<T>()) | |
365 | ) | |
366 | { | |
367 | // Oops, under/overflow, sidestep if we can: | |
368 | if(a < b) | |
369 | { | |
370 | T p1 = pow(b2, b / a); | |
371 | T l3 = a * (log(b1) + log(p1)); | |
372 | if((l3 < tools::log_max_value<T>()) | |
373 | && (l3 > tools::log_min_value<T>())) | |
374 | { | |
375 | result *= pow(p1 * b1, a); | |
376 | } | |
377 | else | |
378 | { | |
379 | l2 += l1 + log(result); | |
380 | if(l2 >= tools::log_max_value<T>()) | |
381 | return policies::raise_overflow_error<T>(function, 0, pol); | |
382 | result = exp(l2); | |
383 | } | |
384 | } | |
385 | else | |
386 | { | |
387 | T p1 = pow(b1, a / b); | |
388 | T l3 = (log(p1) + log(b2)) * b; | |
389 | if((l3 < tools::log_max_value<T>()) | |
390 | && (l3 > tools::log_min_value<T>())) | |
391 | { | |
392 | result *= pow(p1 * b2, b); | |
393 | } | |
394 | else | |
395 | { | |
396 | l2 += l1 + log(result); | |
397 | if(l2 >= tools::log_max_value<T>()) | |
398 | return policies::raise_overflow_error<T>(function, 0, pol); | |
399 | result = exp(l2); | |
400 | } | |
401 | } | |
402 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
403 | } | |
404 | else | |
405 | { | |
406 | // finally the normal case: | |
407 | result *= pow(b1, a) * pow(b2, b); | |
408 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
409 | } | |
410 | } | |
411 | ||
412 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
413 | ||
414 | return result; | |
415 | } | |
416 | // | |
417 | // Compute the leading power terms in the incomplete Beta: | |
418 | // | |
419 | // (x^a)(y^b)/Beta(a,b) when normalised, and | |
420 | // (x^a)(y^b) otherwise. | |
421 | // | |
422 | // Almost all of the error in the incomplete beta comes from this | |
423 | // function: particularly when a and b are large. Computing large | |
424 | // powers are *hard* though, and using logarithms just leads to | |
425 | // horrendous cancellation errors. | |
426 | // | |
427 | // This version is generic, slow, and does not use the Lanczos approximation. | |
428 | // | |
429 | template <class T, class Policy> | |
430 | T ibeta_power_terms(T a, | |
431 | T b, | |
432 | T x, | |
433 | T y, | |
434 | const boost::math::lanczos::undefined_lanczos&, | |
435 | bool normalised, | |
436 | const Policy& pol, | |
437 | T prefix = 1, | |
438 | const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") | |
439 | { | |
440 | BOOST_MATH_STD_USING | |
441 | ||
442 | if(!normalised) | |
443 | { | |
444 | return pow(x, a) * pow(y, b); | |
445 | } | |
446 | ||
447 | T result= 0; // assignment here silences warnings later | |
448 | ||
449 | T c = a + b; | |
450 | ||
451 | // integration limits for the gamma functions: | |
452 | //T la = (std::max)(T(10), a); | |
453 | //T lb = (std::max)(T(10), b); | |
454 | //T lc = (std::max)(T(10), a+b); | |
455 | T la = a + 5; | |
456 | T lb = b + 5; | |
457 | T lc = a + b + 5; | |
458 | // gamma function partials: | |
459 | T sa = detail::lower_gamma_series(a, la, pol) / a; | |
460 | sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>()); | |
461 | T sb = detail::lower_gamma_series(b, lb, pol) / b; | |
462 | sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>()); | |
463 | T sc = detail::lower_gamma_series(c, lc, pol) / c; | |
464 | sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>()); | |
465 | // gamma function powers combined with incomplete beta powers: | |
466 | ||
467 | T b1 = (x * lc) / la; | |
468 | T b2 = (y * lc) / lb; | |
469 | T e1 = -5; // lc - la - lb; | |
470 | T lb1 = a * log(b1); | |
471 | T lb2 = b * log(b2); | |
472 | ||
473 | if((lb1 >= tools::log_max_value<T>()) | |
474 | || (lb1 <= tools::log_min_value<T>()) | |
475 | || (lb2 >= tools::log_max_value<T>()) | |
476 | || (lb2 <= tools::log_min_value<T>()) | |
477 | || (e1 >= tools::log_max_value<T>()) | |
478 | || (e1 <= tools::log_min_value<T>()) | |
479 | ) | |
480 | { | |
481 | result = exp(lb1 + lb2 - e1 + log(prefix)); | |
482 | } | |
483 | else | |
484 | { | |
485 | T p1, p2; | |
486 | p1 = (x * b - y * la) / la; | |
487 | if(fabs(p1) < 0.5f) | |
488 | p1 = exp(a * boost::math::log1p(p1, pol)); | |
489 | else | |
490 | p1 = pow(b1, a); | |
491 | p2 = (y * a - x * lb) / lb; | |
492 | if(fabs(p2) < 0.5f) | |
493 | p2 = exp(b * boost::math::log1p(p2, pol)); | |
494 | else | |
495 | p2 = pow(b2, b); | |
496 | T p3 = exp(e1); | |
497 | result = prefix * p1 * (p2 / p3); | |
498 | } | |
499 | // and combine with the remaining gamma function components: | |
500 | result /= sa * sb / sc; | |
501 | ||
502 | return result; | |
503 | } | |
504 | // | |
505 | // Series approximation to the incomplete beta: | |
506 | // | |
507 | template <class T> | |
508 | struct ibeta_series_t | |
509 | { | |
510 | typedef T result_type; | |
511 | ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {} | |
512 | T operator()() | |
513 | { | |
514 | T r = result / apn; | |
515 | apn += 1; | |
516 | result *= poch * x / n; | |
517 | ++n; | |
518 | poch += 1; | |
519 | return r; | |
520 | } | |
521 | private: | |
522 | T result, x, apn, poch; | |
523 | int n; | |
524 | }; | |
525 | ||
526 | template <class T, class Lanczos, class Policy> | |
527 | T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol) | |
528 | { | |
529 | BOOST_MATH_STD_USING | |
530 | ||
531 | T result; | |
532 | ||
533 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
534 | ||
535 | if(normalised) | |
536 | { | |
537 | T c = a + b; | |
538 | ||
539 | // incomplete beta power term, combined with the Lanczos approximation: | |
540 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
541 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
542 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
543 | result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b)); | |
544 | ||
545 | T l1 = log(cgh / bgh) * (b - 0.5f); | |
546 | T l2 = log(x * cgh / agh) * a; | |
547 | // | |
548 | // Check for over/underflow in the power terms: | |
549 | // | |
550 | if((l1 > tools::log_min_value<T>()) | |
551 | && (l1 < tools::log_max_value<T>()) | |
552 | && (l2 > tools::log_min_value<T>()) | |
553 | && (l2 < tools::log_max_value<T>())) | |
554 | { | |
555 | if(a * b < bgh * 10) | |
556 | result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol)); | |
557 | else | |
558 | result *= pow(cgh / bgh, b - 0.5f); | |
559 | result *= pow(x * cgh / agh, a); | |
560 | result *= sqrt(agh / boost::math::constants::e<T>()); | |
561 | ||
562 | if(p_derivative) | |
563 | { | |
564 | *p_derivative = result * pow(y, b); | |
565 | BOOST_ASSERT(*p_derivative >= 0); | |
566 | } | |
567 | } | |
568 | else | |
569 | { | |
570 | // | |
571 | // Oh dear, we need logs, and this *will* cancel: | |
572 | // | |
573 | result = log(result) + l1 + l2 + (log(agh) - 1) / 2; | |
574 | if(p_derivative) | |
575 | *p_derivative = exp(result + b * log(y)); | |
576 | result = exp(result); | |
577 | } | |
578 | } | |
579 | else | |
580 | { | |
581 | // Non-normalised, just compute the power: | |
582 | result = pow(x, a); | |
583 | } | |
584 | if(result < tools::min_value<T>()) | |
585 | return s0; // Safeguard: series can't cope with denorms. | |
586 | ibeta_series_t<T> s(a, b, x, result); | |
587 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
588 | result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); | |
589 | policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol); | |
590 | return result; | |
591 | } | |
592 | // | |
593 | // Incomplete Beta series again, this time without Lanczos support: | |
594 | // | |
595 | template <class T, class Policy> | |
596 | T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol) | |
597 | { | |
598 | BOOST_MATH_STD_USING | |
599 | ||
600 | T result; | |
601 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
602 | ||
603 | if(normalised) | |
604 | { | |
605 | T c = a + b; | |
606 | ||
607 | // figure out integration limits for the gamma function: | |
608 | //T la = (std::max)(T(10), a); | |
609 | //T lb = (std::max)(T(10), b); | |
610 | //T lc = (std::max)(T(10), a+b); | |
611 | T la = a + 5; | |
612 | T lb = b + 5; | |
613 | T lc = a + b + 5; | |
614 | ||
615 | // calculate the gamma parts: | |
616 | T sa = detail::lower_gamma_series(a, la, pol) / a; | |
617 | sa += detail::upper_gamma_fraction(a, la, ::boost::math::policies::get_epsilon<T, Policy>()); | |
618 | T sb = detail::lower_gamma_series(b, lb, pol) / b; | |
619 | sb += detail::upper_gamma_fraction(b, lb, ::boost::math::policies::get_epsilon<T, Policy>()); | |
620 | T sc = detail::lower_gamma_series(c, lc, pol) / c; | |
621 | sc += detail::upper_gamma_fraction(c, lc, ::boost::math::policies::get_epsilon<T, Policy>()); | |
622 | ||
623 | // and their combined power-terms: | |
624 | T b1 = (x * lc) / la; | |
625 | T b2 = lc/lb; | |
626 | T e1 = lc - la - lb; | |
627 | T lb1 = a * log(b1); | |
628 | T lb2 = b * log(b2); | |
629 | ||
630 | if((lb1 >= tools::log_max_value<T>()) | |
631 | || (lb1 <= tools::log_min_value<T>()) | |
632 | || (lb2 >= tools::log_max_value<T>()) | |
633 | || (lb2 <= tools::log_min_value<T>()) | |
634 | || (e1 >= tools::log_max_value<T>()) | |
635 | || (e1 <= tools::log_min_value<T>()) ) | |
636 | { | |
637 | T p = lb1 + lb2 - e1; | |
638 | result = exp(p); | |
639 | } | |
640 | else | |
641 | { | |
642 | result = pow(b1, a); | |
643 | if(a * b < lb * 10) | |
644 | result *= exp(b * boost::math::log1p(a / lb, pol)); | |
645 | else | |
646 | result *= pow(b2, b); | |
647 | result /= exp(e1); | |
648 | } | |
649 | // and combine the results: | |
650 | result /= sa * sb / sc; | |
651 | ||
652 | if(p_derivative) | |
653 | { | |
654 | *p_derivative = result * pow(y, b); | |
655 | BOOST_ASSERT(*p_derivative >= 0); | |
656 | } | |
657 | } | |
658 | else | |
659 | { | |
660 | // Non-normalised, just compute the power: | |
661 | result = pow(x, a); | |
662 | } | |
663 | if(result < tools::min_value<T>()) | |
664 | return s0; // Safeguard: series can't cope with denorms. | |
665 | ibeta_series_t<T> s(a, b, x, result); | |
666 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
667 | result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); | |
668 | policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol); | |
669 | return result; | |
670 | } | |
671 | ||
672 | // | |
673 | // Continued fraction for the incomplete beta: | |
674 | // | |
675 | template <class T> | |
676 | struct ibeta_fraction2_t | |
677 | { | |
678 | typedef std::pair<T, T> result_type; | |
679 | ||
680 | ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {} | |
681 | ||
682 | result_type operator()() | |
683 | { | |
684 | T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x; | |
685 | T denom = (a + 2 * m - 1); | |
686 | aN /= denom * denom; | |
687 | ||
688 | T bN = static_cast<T>(m); | |
689 | bN += (m * (b - m) * x) / (a + 2*m - 1); | |
690 | bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1); | |
691 | ||
692 | ++m; | |
693 | ||
694 | return std::make_pair(aN, bN); | |
695 | } | |
696 | ||
697 | private: | |
698 | T a, b, x, y; | |
699 | int m; | |
700 | }; | |
701 | // | |
702 | // Evaluate the incomplete beta via the continued fraction representation: | |
703 | // | |
704 | template <class T, class Policy> | |
705 | inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative) | |
706 | { | |
707 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
708 | BOOST_MATH_STD_USING | |
709 | T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); | |
710 | if(p_derivative) | |
711 | { | |
712 | *p_derivative = result; | |
713 | BOOST_ASSERT(*p_derivative >= 0); | |
714 | } | |
715 | if(result == 0) | |
716 | return result; | |
717 | ||
718 | ibeta_fraction2_t<T> f(a, b, x, y); | |
719 | T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>()); | |
720 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
721 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
722 | return result / fract; | |
723 | } | |
724 | // | |
725 | // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x): | |
726 | // | |
727 | template <class T, class Policy> | |
728 | T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative) | |
729 | { | |
730 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
731 | ||
732 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
733 | ||
734 | T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); | |
735 | if(p_derivative) | |
736 | { | |
737 | *p_derivative = prefix; | |
738 | BOOST_ASSERT(*p_derivative >= 0); | |
739 | } | |
740 | prefix /= a; | |
741 | if(prefix == 0) | |
742 | return prefix; | |
743 | T sum = 1; | |
744 | T term = 1; | |
745 | // series summation from 0 to k-1: | |
746 | for(int i = 0; i < k-1; ++i) | |
747 | { | |
748 | term *= (a+b+i) * x / (a+i+1); | |
749 | sum += term; | |
750 | } | |
751 | prefix *= sum; | |
752 | ||
753 | return prefix; | |
754 | } | |
755 | // | |
756 | // This function is only needed for the non-regular incomplete beta, | |
757 | // it computes the delta in: | |
758 | // beta(a,b,x) = prefix + delta * beta(a+k,b,x) | |
759 | // it is currently only called for small k. | |
760 | // | |
761 | template <class T> | |
762 | inline T rising_factorial_ratio(T a, T b, int k) | |
763 | { | |
764 | // calculate: | |
765 | // (a)(a+1)(a+2)...(a+k-1) | |
766 | // _______________________ | |
767 | // (b)(b+1)(b+2)...(b+k-1) | |
768 | ||
769 | // This is only called with small k, for large k | |
770 | // it is grossly inefficient, do not use outside it's | |
771 | // intended purpose!!! | |
772 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
773 | if(k == 0) | |
774 | return 1; | |
775 | T result = 1; | |
776 | for(int i = 0; i < k; ++i) | |
777 | result *= (a+i) / (b+i); | |
778 | return result; | |
779 | } | |
780 | // | |
781 | // Routine for a > 15, b < 1 | |
782 | // | |
783 | // Begin by figuring out how large our table of Pn's should be, | |
784 | // quoted accuracies are "guestimates" based on empiracal observation. | |
785 | // Note that the table size should never exceed the size of our | |
786 | // tables of factorials. | |
787 | // | |
788 | template <class T> | |
789 | struct Pn_size | |
790 | { | |
791 | // This is likely to be enough for ~35-50 digit accuracy | |
792 | // but it's hard to quantify exactly: | |
793 | BOOST_STATIC_CONSTANT(unsigned, value = 50); | |
794 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<T>::value >= 100); | |
795 | }; | |
796 | template <> | |
797 | struct Pn_size<float> | |
798 | { | |
799 | BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy | |
800 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30); | |
801 | }; | |
802 | template <> | |
803 | struct Pn_size<double> | |
804 | { | |
805 | BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy | |
806 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60); | |
807 | }; | |
808 | template <> | |
809 | struct Pn_size<long double> | |
810 | { | |
811 | BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy | |
812 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100); | |
813 | }; | |
814 | ||
815 | template <class T, class Policy> | |
816 | T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised) | |
817 | { | |
818 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
819 | BOOST_MATH_STD_USING | |
820 | // | |
821 | // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6. | |
822 | // | |
823 | // Some values we'll need later, these are Eq 9.1: | |
824 | // | |
825 | T bm1 = b - 1; | |
826 | T t = a + bm1 / 2; | |
827 | T lx, u; | |
828 | if(y < 0.35) | |
829 | lx = boost::math::log1p(-y, pol); | |
830 | else | |
831 | lx = log(x); | |
832 | u = -t * lx; | |
833 | // and from from 9.2: | |
834 | T prefix; | |
835 | T h = regularised_gamma_prefix(b, u, pol, lanczos_type()); | |
836 | if(h <= tools::min_value<T>()) | |
837 | return s0; | |
838 | if(normalised) | |
839 | { | |
840 | prefix = h / boost::math::tgamma_delta_ratio(a, b, pol); | |
841 | prefix /= pow(t, b); | |
842 | } | |
843 | else | |
844 | { | |
845 | prefix = full_igamma_prefix(b, u, pol) / pow(t, b); | |
846 | } | |
847 | prefix *= mult; | |
848 | // | |
849 | // now we need the quantity Pn, unfortunatately this is computed | |
850 | // recursively, and requires a full history of all the previous values | |
851 | // so no choice but to declare a big table and hope it's big enough... | |
852 | // | |
853 | T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3. | |
854 | // | |
855 | // Now an initial value for J, see 9.6: | |
856 | // | |
857 | T j = boost::math::gamma_q(b, u, pol) / h; | |
858 | // | |
859 | // Now we can start to pull things together and evaluate the sum in Eq 9: | |
860 | // | |
861 | T sum = s0 + prefix * j; // Value at N = 0 | |
862 | // some variables we'll need: | |
863 | unsigned tnp1 = 1; // 2*N+1 | |
864 | T lx2 = lx / 2; | |
865 | lx2 *= lx2; | |
866 | T lxp = 1; | |
867 | T t4 = 4 * t * t; | |
868 | T b2n = b; | |
869 | ||
870 | for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n) | |
871 | { | |
872 | /* | |
873 | // debugging code, enable this if you want to determine whether | |
874 | // the table of Pn's is large enough... | |
875 | // | |
876 | static int max_count = 2; | |
877 | if(n > max_count) | |
878 | { | |
879 | max_count = n; | |
880 | std::cerr << "Max iterations in BGRAT was " << n << std::endl; | |
881 | } | |
882 | */ | |
883 | // | |
884 | // begin by evaluating the next Pn from Eq 9.4: | |
885 | // | |
886 | tnp1 += 2; | |
887 | p[n] = 0; | |
888 | T mbn = b - n; | |
889 | unsigned tmp1 = 3; | |
890 | for(unsigned m = 1; m < n; ++m) | |
891 | { | |
892 | mbn = m * b - n; | |
893 | p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1); | |
894 | tmp1 += 2; | |
895 | } | |
896 | p[n] /= n; | |
897 | p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1); | |
898 | // | |
899 | // Now we want Jn from Jn-1 using Eq 9.6: | |
900 | // | |
901 | j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4; | |
902 | lxp *= lx2; | |
903 | b2n += 2; | |
904 | // | |
905 | // pull it together with Eq 9: | |
906 | // | |
907 | T r = prefix * p[n] * j; | |
908 | sum += r; | |
909 | if(r > 1) | |
910 | { | |
911 | if(fabs(r) < fabs(tools::epsilon<T>() * sum)) | |
912 | break; | |
913 | } | |
914 | else | |
915 | { | |
916 | if(fabs(r / tools::epsilon<T>()) < fabs(sum)) | |
917 | break; | |
918 | } | |
919 | } | |
920 | return sum; | |
921 | } // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised) | |
922 | ||
923 | // | |
924 | // For integer arguments we can relate the incomplete beta to the | |
925 | // complement of the binomial distribution cdf and use this finite sum. | |
926 | // | |
927 | template <class T> | |
928 | T binomial_ccdf(T n, T k, T x, T y) | |
929 | { | |
930 | BOOST_MATH_STD_USING // ADL of std names | |
931 | ||
932 | T result = pow(x, n); | |
933 | ||
934 | if(result > tools::min_value<T>()) | |
935 | { | |
936 | T term = result; | |
937 | for(unsigned i = itrunc(T(n - 1)); i > k; --i) | |
938 | { | |
939 | term *= ((i + 1) * y) / ((n - i) * x); | |
940 | result += term; | |
941 | } | |
942 | } | |
943 | else | |
944 | { | |
945 | // First term underflows so we need to start at the mode of the | |
946 | // distribution and work outwards: | |
947 | int start = itrunc(n * x); | |
948 | if(start <= k + 1) | |
949 | start = itrunc(k + 2); | |
950 | result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start)); | |
951 | if(result == 0) | |
952 | { | |
953 | // OK, starting slightly above the mode didn't work, | |
954 | // we'll have to sum the terms the old fashioned way: | |
955 | for(unsigned i = start - 1; i > k; --i) | |
956 | { | |
957 | result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i)); | |
958 | } | |
959 | } | |
960 | else | |
961 | { | |
962 | T term = result; | |
963 | T start_term = result; | |
964 | for(unsigned i = start - 1; i > k; --i) | |
965 | { | |
966 | term *= ((i + 1) * y) / ((n - i) * x); | |
967 | result += term; | |
968 | } | |
969 | term = start_term; | |
970 | for(unsigned i = start + 1; i <= n; ++i) | |
971 | { | |
972 | term *= (n - i + 1) * x / (i * y); | |
973 | result += term; | |
974 | } | |
975 | } | |
976 | } | |
977 | ||
978 | return result; | |
979 | } | |
980 | ||
981 | ||
982 | // | |
983 | // The incomplete beta function implementation: | |
984 | // This is just a big bunch of spagetti code to divide up the | |
985 | // input range and select the right implementation method for | |
986 | // each domain: | |
987 | // | |
988 | template <class T, class Policy> | |
989 | T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative) | |
990 | { | |
991 | static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; | |
992 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
993 | BOOST_MATH_STD_USING // for ADL of std math functions. | |
994 | ||
995 | BOOST_MATH_INSTRUMENT_VARIABLE(a); | |
996 | BOOST_MATH_INSTRUMENT_VARIABLE(b); | |
997 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
998 | BOOST_MATH_INSTRUMENT_VARIABLE(inv); | |
999 | BOOST_MATH_INSTRUMENT_VARIABLE(normalised); | |
1000 | ||
1001 | bool invert = inv; | |
1002 | T fract; | |
1003 | T y = 1 - x; | |
1004 | ||
1005 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
1006 | ||
1007 | if(p_derivative) | |
1008 | *p_derivative = -1; // value not set. | |
1009 | ||
1010 | if((x < 0) || (x > 1)) | |
1011 | return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); | |
1012 | ||
1013 | if(normalised) | |
1014 | { | |
1015 | if(a < 0) | |
1016 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); | |
1017 | if(b < 0) | |
1018 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); | |
1019 | // extend to a few very special cases: | |
1020 | if(a == 0) | |
1021 | { | |
1022 | if(b == 0) | |
1023 | return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol); | |
1024 | if(b > 0) | |
1025 | return static_cast<T>(inv ? 0 : 1); | |
1026 | } | |
1027 | else if(b == 0) | |
1028 | { | |
1029 | if(a > 0) | |
1030 | return static_cast<T>(inv ? 1 : 0); | |
1031 | } | |
1032 | } | |
1033 | else | |
1034 | { | |
1035 | if(a <= 0) | |
1036 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); | |
1037 | if(b <= 0) | |
1038 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); | |
1039 | } | |
1040 | ||
1041 | if(x == 0) | |
1042 | { | |
1043 | if(p_derivative) | |
1044 | { | |
1045 | *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); | |
1046 | } | |
1047 | return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0)); | |
1048 | } | |
1049 | if(x == 1) | |
1050 | { | |
1051 | if(p_derivative) | |
1052 | { | |
1053 | *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); | |
1054 | } | |
1055 | return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0); | |
1056 | } | |
1057 | if((a == 0.5f) && (b == 0.5f)) | |
1058 | { | |
1059 | // We have an arcsine distribution: | |
1060 | if(p_derivative) | |
1061 | { | |
1062 | *p_derivative = 1 / constants::pi<T>() * sqrt(y * x); | |
1063 | } | |
1064 | T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>(); | |
1065 | if(!normalised) | |
1066 | p *= constants::pi<T>(); | |
1067 | return p; | |
1068 | } | |
1069 | if(a == 1) | |
1070 | { | |
1071 | std::swap(a, b); | |
1072 | std::swap(x, y); | |
1073 | invert = !invert; | |
1074 | } | |
1075 | if(b == 1) | |
1076 | { | |
1077 | // | |
1078 | // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/ | |
1079 | // | |
1080 | if(a == 1) | |
1081 | { | |
1082 | if(p_derivative) | |
1083 | *p_derivative = 1; | |
1084 | return invert ? y : x; | |
1085 | } | |
1086 | ||
1087 | if(p_derivative) | |
1088 | { | |
1089 | *p_derivative = a * pow(x, a - 1); | |
1090 | } | |
1091 | T p; | |
1092 | if(y < 0.5) | |
1093 | p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol))); | |
1094 | else | |
1095 | p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a)); | |
1096 | if(!normalised) | |
1097 | p /= a; | |
1098 | return p; | |
1099 | } | |
1100 | ||
1101 | if((std::min)(a, b) <= 1) | |
1102 | { | |
1103 | if(x > 0.5) | |
1104 | { | |
1105 | std::swap(a, b); | |
1106 | std::swap(x, y); | |
1107 | invert = !invert; | |
1108 | BOOST_MATH_INSTRUMENT_VARIABLE(invert); | |
1109 | } | |
1110 | if((std::max)(a, b) <= 1) | |
1111 | { | |
1112 | // Both a,b < 1: | |
1113 | if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9)) | |
1114 | { | |
1115 | if(!invert) | |
1116 | { | |
1117 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1118 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1119 | } | |
1120 | else | |
1121 | { | |
1122 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1123 | invert = false; | |
1124 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1125 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1126 | } | |
1127 | } | |
1128 | else | |
1129 | { | |
1130 | std::swap(a, b); | |
1131 | std::swap(x, y); | |
1132 | invert = !invert; | |
1133 | if(y >= 0.3) | |
1134 | { | |
1135 | if(!invert) | |
1136 | { | |
1137 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1138 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1139 | } | |
1140 | else | |
1141 | { | |
1142 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1143 | invert = false; | |
1144 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1145 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1146 | } | |
1147 | } | |
1148 | else | |
1149 | { | |
1150 | // Sidestep on a, and then use the series representation: | |
1151 | T prefix; | |
1152 | if(!normalised) | |
1153 | { | |
1154 | prefix = rising_factorial_ratio(T(a+b), a, 20); | |
1155 | } | |
1156 | else | |
1157 | { | |
1158 | prefix = 1; | |
1159 | } | |
1160 | fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); | |
1161 | if(!invert) | |
1162 | { | |
1163 | fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1164 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1165 | } | |
1166 | else | |
1167 | { | |
1168 | fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); | |
1169 | invert = false; | |
1170 | fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1171 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1172 | } | |
1173 | } | |
1174 | } | |
1175 | } | |
1176 | else | |
1177 | { | |
1178 | // One of a, b < 1 only: | |
1179 | if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7))) | |
1180 | { | |
1181 | if(!invert) | |
1182 | { | |
1183 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1184 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1185 | } | |
1186 | else | |
1187 | { | |
1188 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1189 | invert = false; | |
1190 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1191 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1192 | } | |
1193 | } | |
1194 | else | |
1195 | { | |
1196 | std::swap(a, b); | |
1197 | std::swap(x, y); | |
1198 | invert = !invert; | |
1199 | ||
1200 | if(y >= 0.3) | |
1201 | { | |
1202 | if(!invert) | |
1203 | { | |
1204 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1205 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1206 | } | |
1207 | else | |
1208 | { | |
1209 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1210 | invert = false; | |
1211 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1212 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1213 | } | |
1214 | } | |
1215 | else if(a >= 15) | |
1216 | { | |
1217 | if(!invert) | |
1218 | { | |
1219 | fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised); | |
1220 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1221 | } | |
1222 | else | |
1223 | { | |
1224 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1225 | invert = false; | |
1226 | fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised); | |
1227 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1228 | } | |
1229 | } | |
1230 | else | |
1231 | { | |
1232 | // Sidestep to improve errors: | |
1233 | T prefix; | |
1234 | if(!normalised) | |
1235 | { | |
1236 | prefix = rising_factorial_ratio(T(a+b), a, 20); | |
1237 | } | |
1238 | else | |
1239 | { | |
1240 | prefix = 1; | |
1241 | } | |
1242 | fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); | |
1243 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1244 | if(!invert) | |
1245 | { | |
1246 | fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1247 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1248 | } | |
1249 | else | |
1250 | { | |
1251 | fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); | |
1252 | invert = false; | |
1253 | fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1254 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1255 | } | |
1256 | } | |
1257 | } | |
1258 | } | |
1259 | } | |
1260 | else | |
1261 | { | |
1262 | // Both a,b >= 1: | |
1263 | T lambda; | |
1264 | if(a < b) | |
1265 | { | |
1266 | lambda = a - (a + b) * x; | |
1267 | } | |
1268 | else | |
1269 | { | |
1270 | lambda = (a + b) * y - b; | |
1271 | } | |
1272 | if(lambda < 0) | |
1273 | { | |
1274 | std::swap(a, b); | |
1275 | std::swap(x, y); | |
1276 | invert = !invert; | |
1277 | BOOST_MATH_INSTRUMENT_VARIABLE(invert); | |
1278 | } | |
1279 | ||
1280 | if(b < 40) | |
1281 | { | |
1282 | if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100) && (y != 1)) | |
1283 | { | |
1284 | // relate to the binomial distribution and use a finite sum: | |
1285 | T k = a - 1; | |
1286 | T n = b + k; | |
1287 | fract = binomial_ccdf(n, k, x, y); | |
1288 | if(!normalised) | |
1289 | fract *= boost::math::beta(a, b, pol); | |
1290 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1291 | } | |
1292 | else if(b * x <= 0.7) | |
1293 | { | |
1294 | if(!invert) | |
1295 | { | |
1296 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1297 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1298 | } | |
1299 | else | |
1300 | { | |
1301 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1302 | invert = false; | |
1303 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1304 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1305 | } | |
1306 | } | |
1307 | else if(a > 15) | |
1308 | { | |
1309 | // sidestep so we can use the series representation: | |
1310 | int n = itrunc(T(floor(b)), pol); | |
1311 | if(n == b) | |
1312 | --n; | |
1313 | T bbar = b - n; | |
1314 | T prefix; | |
1315 | if(!normalised) | |
1316 | { | |
1317 | prefix = rising_factorial_ratio(T(a+bbar), bbar, n); | |
1318 | } | |
1319 | else | |
1320 | { | |
1321 | prefix = 1; | |
1322 | } | |
1323 | fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); | |
1324 | fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised); | |
1325 | fract /= prefix; | |
1326 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1327 | } | |
1328 | else if(normalised) | |
1329 | { | |
1330 | // The formula here for the non-normalised case is tricky to figure | |
1331 | // out (for me!!), and requires two pochhammer calculations rather | |
1332 | // than one, so leave it for now and only use this in the normalized case.... | |
1333 | int n = itrunc(T(floor(b)), pol); | |
1334 | T bbar = b - n; | |
1335 | if(bbar <= 0) | |
1336 | { | |
1337 | --n; | |
1338 | bbar += 1; | |
1339 | } | |
1340 | fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); | |
1341 | fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0)); | |
1342 | if(invert) | |
1343 | fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case | |
1344 | fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised); | |
1345 | if(invert) | |
1346 | { | |
1347 | fract = -fract; | |
1348 | invert = false; | |
1349 | } | |
1350 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1351 | } | |
1352 | else | |
1353 | { | |
1354 | fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); | |
1355 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1356 | } | |
1357 | } | |
1358 | else | |
1359 | { | |
1360 | fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); | |
1361 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1362 | } | |
1363 | } | |
1364 | if(p_derivative) | |
1365 | { | |
1366 | if(*p_derivative < 0) | |
1367 | { | |
1368 | *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol); | |
1369 | } | |
1370 | T div = y * x; | |
1371 | ||
1372 | if(*p_derivative != 0) | |
1373 | { | |
1374 | if((tools::max_value<T>() * div < *p_derivative)) | |
1375 | { | |
1376 | // overflow, return an arbitarily large value: | |
1377 | *p_derivative = tools::max_value<T>() / 2; | |
1378 | } | |
1379 | else | |
1380 | { | |
1381 | *p_derivative /= div; | |
1382 | } | |
1383 | } | |
1384 | } | |
1385 | return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; | |
1386 | } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised) | |
1387 | ||
1388 | template <class T, class Policy> | |
1389 | inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised) | |
1390 | { | |
1391 | return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(0)); | |
1392 | } | |
1393 | ||
1394 | template <class T, class Policy> | |
1395 | T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) | |
1396 | { | |
1397 | static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)"; | |
1398 | // | |
1399 | // start with the usual error checks: | |
1400 | // | |
1401 | if(a <= 0) | |
1402 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); | |
1403 | if(b <= 0) | |
1404 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); | |
1405 | if((x < 0) || (x > 1)) | |
1406 | return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); | |
1407 | // | |
1408 | // Now the corner cases: | |
1409 | // | |
1410 | if(x == 0) | |
1411 | { | |
1412 | return (a > 1) ? 0 : | |
1413 | (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol); | |
1414 | } | |
1415 | else if(x == 1) | |
1416 | { | |
1417 | return (b > 1) ? 0 : | |
1418 | (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol); | |
1419 | } | |
1420 | // | |
1421 | // Now the regular cases: | |
1422 | // | |
1423 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
1424 | T y = (1 - x) * x; | |
1425 | T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function); | |
1426 | return f1; | |
1427 | } | |
1428 | // | |
1429 | // Some forwarding functions that dis-ambiguate the third argument type: | |
1430 | // | |
1431 | template <class RT1, class RT2, class Policy> | |
1432 | inline typename tools::promote_args<RT1, RT2>::type | |
1433 | beta(RT1 a, RT2 b, const Policy&, const mpl::true_*) | |
1434 | { | |
1435 | BOOST_FPU_EXCEPTION_GUARD | |
1436 | typedef typename tools::promote_args<RT1, RT2>::type result_type; | |
1437 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1438 | typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; | |
1439 | typedef typename policies::normalise< | |
1440 | Policy, | |
1441 | policies::promote_float<false>, | |
1442 | policies::promote_double<false>, | |
1443 | policies::discrete_quantile<>, | |
1444 | policies::assert_undefined<> >::type forwarding_policy; | |
1445 | ||
1446 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::beta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::beta<%1%>(%1%,%1%)"); | |
1447 | } | |
1448 | template <class RT1, class RT2, class RT3> | |
1449 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1450 | beta(RT1 a, RT2 b, RT3 x, const mpl::false_*) | |
1451 | { | |
1452 | return boost::math::beta(a, b, x, policies::policy<>()); | |
1453 | } | |
1454 | } // namespace detail | |
1455 | ||
1456 | // | |
1457 | // The actual function entry-points now follow, these just figure out | |
1458 | // which Lanczos approximation to use | |
1459 | // and forward to the implementation functions: | |
1460 | // | |
1461 | template <class RT1, class RT2, class A> | |
1462 | inline typename tools::promote_args<RT1, RT2, A>::type | |
1463 | beta(RT1 a, RT2 b, A arg) | |
1464 | { | |
1465 | typedef typename policies::is_policy<A>::type tag; | |
1466 | return boost::math::detail::beta(a, b, arg, static_cast<tag*>(0)); | |
1467 | } | |
1468 | ||
1469 | template <class RT1, class RT2> | |
1470 | inline typename tools::promote_args<RT1, RT2>::type | |
1471 | beta(RT1 a, RT2 b) | |
1472 | { | |
1473 | return boost::math::beta(a, b, policies::policy<>()); | |
1474 | } | |
1475 | ||
1476 | template <class RT1, class RT2, class RT3, class Policy> | |
1477 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1478 | beta(RT1 a, RT2 b, RT3 x, const Policy&) | |
1479 | { | |
1480 | BOOST_FPU_EXCEPTION_GUARD | |
1481 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1482 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1483 | typedef typename policies::normalise< | |
1484 | Policy, | |
1485 | policies::promote_float<false>, | |
1486 | policies::promote_double<false>, | |
1487 | policies::discrete_quantile<>, | |
1488 | policies::assert_undefined<> >::type forwarding_policy; | |
1489 | ||
1490 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, false), "boost::math::beta<%1%>(%1%,%1%,%1%)"); | |
1491 | } | |
1492 | ||
1493 | template <class RT1, class RT2, class RT3, class Policy> | |
1494 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1495 | betac(RT1 a, RT2 b, RT3 x, const Policy&) | |
1496 | { | |
1497 | BOOST_FPU_EXCEPTION_GUARD | |
1498 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1499 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1500 | typedef typename policies::normalise< | |
1501 | Policy, | |
1502 | policies::promote_float<false>, | |
1503 | policies::promote_double<false>, | |
1504 | policies::discrete_quantile<>, | |
1505 | policies::assert_undefined<> >::type forwarding_policy; | |
1506 | ||
1507 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, false), "boost::math::betac<%1%>(%1%,%1%,%1%)"); | |
1508 | } | |
1509 | template <class RT1, class RT2, class RT3> | |
1510 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1511 | betac(RT1 a, RT2 b, RT3 x) | |
1512 | { | |
1513 | return boost::math::betac(a, b, x, policies::policy<>()); | |
1514 | } | |
1515 | ||
1516 | template <class RT1, class RT2, class RT3, class Policy> | |
1517 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1518 | ibeta(RT1 a, RT2 b, RT3 x, const Policy&) | |
1519 | { | |
1520 | BOOST_FPU_EXCEPTION_GUARD | |
1521 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1522 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1523 | typedef typename policies::normalise< | |
1524 | Policy, | |
1525 | policies::promote_float<false>, | |
1526 | policies::promote_double<false>, | |
1527 | policies::discrete_quantile<>, | |
1528 | policies::assert_undefined<> >::type forwarding_policy; | |
1529 | ||
1530 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)"); | |
1531 | } | |
1532 | template <class RT1, class RT2, class RT3> | |
1533 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1534 | ibeta(RT1 a, RT2 b, RT3 x) | |
1535 | { | |
1536 | return boost::math::ibeta(a, b, x, policies::policy<>()); | |
1537 | } | |
1538 | ||
1539 | template <class RT1, class RT2, class RT3, class Policy> | |
1540 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1541 | ibetac(RT1 a, RT2 b, RT3 x, const Policy&) | |
1542 | { | |
1543 | BOOST_FPU_EXCEPTION_GUARD | |
1544 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1545 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1546 | typedef typename policies::normalise< | |
1547 | Policy, | |
1548 | policies::promote_float<false>, | |
1549 | policies::promote_double<false>, | |
1550 | policies::discrete_quantile<>, | |
1551 | policies::assert_undefined<> >::type forwarding_policy; | |
1552 | ||
1553 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, true), "boost::math::ibetac<%1%>(%1%,%1%,%1%)"); | |
1554 | } | |
1555 | template <class RT1, class RT2, class RT3> | |
1556 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1557 | ibetac(RT1 a, RT2 b, RT3 x) | |
1558 | { | |
1559 | return boost::math::ibetac(a, b, x, policies::policy<>()); | |
1560 | } | |
1561 | ||
1562 | template <class RT1, class RT2, class RT3, class Policy> | |
1563 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1564 | ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&) | |
1565 | { | |
1566 | BOOST_FPU_EXCEPTION_GUARD | |
1567 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1568 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1569 | typedef typename policies::normalise< | |
1570 | Policy, | |
1571 | policies::promote_float<false>, | |
1572 | policies::promote_double<false>, | |
1573 | policies::discrete_quantile<>, | |
1574 | policies::assert_undefined<> >::type forwarding_policy; | |
1575 | ||
1576 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy()), "boost::math::ibeta_derivative<%1%>(%1%,%1%,%1%)"); | |
1577 | } | |
1578 | template <class RT1, class RT2, class RT3> | |
1579 | inline typename tools::promote_args<RT1, RT2, RT3>::type | |
1580 | ibeta_derivative(RT1 a, RT2 b, RT3 x) | |
1581 | { | |
1582 | return boost::math::ibeta_derivative(a, b, x, policies::policy<>()); | |
1583 | } | |
1584 | ||
1585 | } // namespace math | |
1586 | } // namespace boost | |
1587 | ||
1588 | #include <boost/math/special_functions/detail/ibeta_inverse.hpp> | |
1589 | #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp> | |
1590 | ||
1591 | #endif // BOOST_MATH_SPECIAL_BETA_HPP | |
1592 | ||
1593 | ||
1594 | ||
1595 | ||
1596 |