1 // boost\math\special_functions\negative_binomial.hpp 2 3 // Copyright Paul A. Bristow 2007. 4 // Copyright John Maddock 2007. 5 6 // Use, modification and distribution are subject to the 7 // Boost Software License, Version 1.0. 8 // (See accompanying file LICENSE_1_0.txt 9 // or copy at http://www.boost.org/LICENSE_1_0.txt) 10 11 // http://en.wikipedia.org/wiki/negative_binomial_distribution 12 // http://mathworld.wolfram.com/NegativeBinomialDistribution.html 13 // http://documents.wolfram.com/teachersedition/Teacher/Statistics/DiscreteDistributions.html 14 15 // The negative binomial distribution NegativeBinomialDistribution[n, p] 16 // is the distribution of the number (k) of failures that occur in a sequence of trials before 17 // r successes have occurred, where the probability of success in each trial is p. 18 19 // In a sequence of Bernoulli trials or events 20 // (independent, yes or no, succeed or fail) with success_fraction probability p, 21 // negative_binomial is the probability that k or fewer failures 22 // preceed the r th trial's success. 23 // random variable k is the number of failures (NOT the probability). 24 25 // Negative_binomial distribution is a discrete probability distribution. 26 // But note that the negative binomial distribution 27 // (like others including the binomial, Poisson & Bernoulli) 28 // is strictly defined as a discrete function: only integral values of k are envisaged. 29 // However because of the method of calculation using a continuous gamma function, 30 // it is convenient to treat it as if a continous function, 31 // and permit non-integral values of k. 32 33 // However, by default the policy is to use discrete_quantile_policy. 34 35 // To enforce the strict mathematical model, users should use conversion 36 // on k outside this function to ensure that k is integral. 37 38 // MATHCAD cumulative negative binomial pnbinom(k, n, p) 39 40 // Implementation note: much greater speed, and perhaps greater accuracy, 41 // might be achieved for extreme values by using a normal approximation. 42 // This is NOT been tested or implemented. 43 44 #ifndef BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP 45 #define BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP 46 47 #include <boost/math/distributions/fwd.hpp> 48 #include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b). 49 #include <boost/math/distributions/complement.hpp> // complement. 50 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error. 51 #include <boost/math/special_functions/fpclassify.hpp> // isnan. 52 #include <boost/math/tools/roots.hpp> // for root finding. 53 #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> 54 55 #include <boost/type_traits/is_floating_point.hpp> 56 #include <boost/type_traits/is_integral.hpp> 57 #include <boost/type_traits/is_same.hpp> 58 #include <boost/mpl/if.hpp> 59 60 #include <limits> // using std::numeric_limits; 61 #include <utility> 62 63 #if defined (BOOST_MSVC) 64 # pragma warning(push) 65 // This believed not now necessary, so commented out. 66 //# pragma warning(disable: 4702) // unreachable code. 67 // in domain_error_imp in error_handling. 68 #endif 69 70 namespace boost 71 { 72 namespace math 73 { 74 namespace negative_binomial_detail 75 { 76 // Common error checking routines for negative binomial distribution functions: 77 template <class RealType, class Policy> check_successes(const char * function,const RealType & r,RealType * result,const Policy & pol)78 inline bool check_successes(const char* function, const RealType& r, RealType* result, const Policy& pol) 79 { 80 if( !(boost::math::isfinite)(r) || (r <= 0) ) 81 { 82 *result = policies::raise_domain_error<RealType>( 83 function, 84 "Number of successes argument is %1%, but must be > 0 !", r, pol); 85 return false; 86 } 87 return true; 88 } 89 template <class RealType, class Policy> check_success_fraction(const char * function,const RealType & p,RealType * result,const Policy & pol)90 inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol) 91 { 92 if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) ) 93 { 94 *result = policies::raise_domain_error<RealType>( 95 function, 96 "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol); 97 return false; 98 } 99 return true; 100 } 101 template <class RealType, class Policy> check_dist(const char * function,const RealType & r,const RealType & p,RealType * result,const Policy & pol)102 inline bool check_dist(const char* function, const RealType& r, const RealType& p, RealType* result, const Policy& pol) 103 { 104 return check_success_fraction(function, p, result, pol) 105 && check_successes(function, r, result, pol); 106 } 107 template <class RealType, class Policy> check_dist_and_k(const char * function,const RealType & r,const RealType & p,RealType k,RealType * result,const Policy & pol)108 inline bool check_dist_and_k(const char* function, const RealType& r, const RealType& p, RealType k, RealType* result, const Policy& pol) 109 { 110 if(check_dist(function, r, p, result, pol) == false) 111 { 112 return false; 113 } 114 if( !(boost::math::isfinite)(k) || (k < 0) ) 115 { // Check k failures. 116 *result = policies::raise_domain_error<RealType>( 117 function, 118 "Number of failures argument is %1%, but must be >= 0 !", k, pol); 119 return false; 120 } 121 return true; 122 } // Check_dist_and_k 123 124 template <class RealType, class Policy> check_dist_and_prob(const char * function,const RealType & r,RealType p,RealType prob,RealType * result,const Policy & pol)125 inline bool check_dist_and_prob(const char* function, const RealType& r, RealType p, RealType prob, RealType* result, const Policy& pol) 126 { 127 if(check_dist(function, r, p, result, pol) && detail::check_probability(function, prob, result, pol) == false) 128 { 129 return false; 130 } 131 return true; 132 } // check_dist_and_prob 133 } // namespace negative_binomial_detail 134 135 template <class RealType = double, class Policy = policies::policy<> > 136 class negative_binomial_distribution 137 { 138 public: 139 typedef RealType value_type; 140 typedef Policy policy_type; 141 negative_binomial_distribution(RealType r,RealType p)142 negative_binomial_distribution(RealType r, RealType p) : m_r(r), m_p(p) 143 { // Constructor. 144 RealType result; 145 negative_binomial_detail::check_dist( 146 "negative_binomial_distribution<%1%>::negative_binomial_distribution", 147 m_r, // Check successes r > 0. 148 m_p, // Check success_fraction 0 <= p <= 1. 149 &result, Policy()); 150 } // negative_binomial_distribution constructor. 151 152 // Private data getter class member functions. success_fraction() const153 RealType success_fraction() const 154 { // Probability of success as fraction in range 0 to 1. 155 return m_p; 156 } successes() const157 RealType successes() const 158 { // Total number of successes r. 159 return m_r; 160 } 161 find_lower_bound_on_p(RealType trials,RealType successes,RealType alpha)162 static RealType find_lower_bound_on_p( 163 RealType trials, 164 RealType successes, 165 RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. 166 { 167 static const char* function = "boost::math::negative_binomial<%1%>::find_lower_bound_on_p"; 168 RealType result = 0; // of error checks. 169 RealType failures = trials - successes; 170 if(false == detail::check_probability(function, alpha, &result, Policy()) 171 && negative_binomial_detail::check_dist_and_k( 172 function, successes, RealType(0), failures, &result, Policy())) 173 { 174 return result; 175 } 176 // Use complement ibeta_inv function for lower bound. 177 // This is adapted from the corresponding binomial formula 178 // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm 179 // This is a Clopper-Pearson interval, and may be overly conservative, 180 // see also "A Simple Improved Inferential Method for Some 181 // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY 182 // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf 183 // 184 return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(0), Policy()); 185 } // find_lower_bound_on_p 186 find_upper_bound_on_p(RealType trials,RealType successes,RealType alpha)187 static RealType find_upper_bound_on_p( 188 RealType trials, 189 RealType successes, 190 RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test. 191 { 192 static const char* function = "boost::math::negative_binomial<%1%>::find_upper_bound_on_p"; 193 RealType result = 0; // of error checks. 194 RealType failures = trials - successes; 195 if(false == negative_binomial_detail::check_dist_and_k( 196 function, successes, RealType(0), failures, &result, Policy()) 197 && detail::check_probability(function, alpha, &result, Policy())) 198 { 199 return result; 200 } 201 if(failures == 0) 202 return 1; 203 // Use complement ibetac_inv function for upper bound. 204 // Note adjusted failures value: *not* failures+1 as usual. 205 // This is adapted from the corresponding binomial formula 206 // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm 207 // This is a Clopper-Pearson interval, and may be overly conservative, 208 // see also "A Simple Improved Inferential Method for Some 209 // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY 210 // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf 211 // 212 return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(0), Policy()); 213 } // find_upper_bound_on_p 214 215 // Estimate number of trials : 216 // "How many trials do I need to be P% sure of seeing k or fewer failures?" 217 find_minimum_number_of_trials(RealType k,RealType p,RealType alpha)218 static RealType find_minimum_number_of_trials( 219 RealType k, // number of failures (k >= 0). 220 RealType p, // success fraction 0 <= p <= 1. 221 RealType alpha) // risk level threshold 0 <= alpha <= 1. 222 { 223 static const char* function = "boost::math::negative_binomial<%1%>::find_minimum_number_of_trials"; 224 // Error checks: 225 RealType result = 0; 226 if(false == negative_binomial_detail::check_dist_and_k( 227 function, RealType(1), p, k, &result, Policy()) 228 && detail::check_probability(function, alpha, &result, Policy())) 229 { return result; } 230 231 result = ibeta_inva(k + 1, p, alpha, Policy()); // returns n - k 232 return result + k; 233 } // RealType find_number_of_failures 234 find_maximum_number_of_trials(RealType k,RealType p,RealType alpha)235 static RealType find_maximum_number_of_trials( 236 RealType k, // number of failures (k >= 0). 237 RealType p, // success fraction 0 <= p <= 1. 238 RealType alpha) // risk level threshold 0 <= alpha <= 1. 239 { 240 static const char* function = "boost::math::negative_binomial<%1%>::find_maximum_number_of_trials"; 241 // Error checks: 242 RealType result = 0; 243 if(false == negative_binomial_detail::check_dist_and_k( 244 function, RealType(1), p, k, &result, Policy()) 245 && detail::check_probability(function, alpha, &result, Policy())) 246 { return result; } 247 248 result = ibetac_inva(k + 1, p, alpha, Policy()); // returns n - k 249 return result + k; 250 } // RealType find_number_of_trials complemented 251 252 private: 253 RealType m_r; // successes. 254 RealType m_p; // success_fraction 255 }; // template <class RealType, class Policy> class negative_binomial_distribution 256 257 typedef negative_binomial_distribution<double> negative_binomial; // Reserved name of type double. 258 259 template <class RealType, class Policy> range(const negative_binomial_distribution<RealType,Policy> &)260 inline const std::pair<RealType, RealType> range(const negative_binomial_distribution<RealType, Policy>& /* dist */) 261 { // Range of permissible values for random variable k. 262 using boost::math::tools::max_value; 263 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? 264 } 265 266 template <class RealType, class Policy> support(const negative_binomial_distribution<RealType,Policy> &)267 inline const std::pair<RealType, RealType> support(const negative_binomial_distribution<RealType, Policy>& /* dist */) 268 { // Range of supported values for random variable k. 269 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 270 using boost::math::tools::max_value; 271 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer? 272 } 273 274 template <class RealType, class Policy> mean(const negative_binomial_distribution<RealType,Policy> & dist)275 inline RealType mean(const negative_binomial_distribution<RealType, Policy>& dist) 276 { // Mean of Negative Binomial distribution = r(1-p)/p. 277 return dist.successes() * (1 - dist.success_fraction() ) / dist.success_fraction(); 278 } // mean 279 280 //template <class RealType, class Policy> 281 //inline RealType median(const negative_binomial_distribution<RealType, Policy>& dist) 282 //{ // Median of negative_binomial_distribution is not defined. 283 // return policies::raise_domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); 284 //} // median 285 // Now implemented via quantile(half) in derived accessors. 286 287 template <class RealType, class Policy> mode(const negative_binomial_distribution<RealType,Policy> & dist)288 inline RealType mode(const negative_binomial_distribution<RealType, Policy>& dist) 289 { // Mode of Negative Binomial distribution = floor[(r-1) * (1 - p)/p] 290 BOOST_MATH_STD_USING // ADL of std functions. 291 return floor((dist.successes() -1) * (1 - dist.success_fraction()) / dist.success_fraction()); 292 } // mode 293 294 template <class RealType, class Policy> skewness(const negative_binomial_distribution<RealType,Policy> & dist)295 inline RealType skewness(const negative_binomial_distribution<RealType, Policy>& dist) 296 { // skewness of Negative Binomial distribution = 2-p / (sqrt(r(1-p)) 297 BOOST_MATH_STD_USING // ADL of std functions. 298 RealType p = dist.success_fraction(); 299 RealType r = dist.successes(); 300 301 return (2 - p) / 302 sqrt(r * (1 - p)); 303 } // skewness 304 305 template <class RealType, class Policy> kurtosis(const negative_binomial_distribution<RealType,Policy> & dist)306 inline RealType kurtosis(const negative_binomial_distribution<RealType, Policy>& dist) 307 { // kurtosis of Negative Binomial distribution 308 // http://en.wikipedia.org/wiki/Negative_binomial is kurtosis_excess so add 3 309 RealType p = dist.success_fraction(); 310 RealType r = dist.successes(); 311 return 3 + (6 / r) + ((p * p) / (r * (1 - p))); 312 } // kurtosis 313 314 template <class RealType, class Policy> kurtosis_excess(const negative_binomial_distribution<RealType,Policy> & dist)315 inline RealType kurtosis_excess(const negative_binomial_distribution<RealType, Policy>& dist) 316 { // kurtosis excess of Negative Binomial distribution 317 // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess 318 RealType p = dist.success_fraction(); 319 RealType r = dist.successes(); 320 return (6 - p * (6-p)) / (r * (1-p)); 321 } // kurtosis_excess 322 323 template <class RealType, class Policy> variance(const negative_binomial_distribution<RealType,Policy> & dist)324 inline RealType variance(const negative_binomial_distribution<RealType, Policy>& dist) 325 { // Variance of Binomial distribution = r (1-p) / p^2. 326 return dist.successes() * (1 - dist.success_fraction()) 327 / (dist.success_fraction() * dist.success_fraction()); 328 } // variance 329 330 // RealType standard_deviation(const negative_binomial_distribution<RealType, Policy>& dist) 331 // standard_deviation provided by derived accessors. 332 // RealType hazard(const negative_binomial_distribution<RealType, Policy>& dist) 333 // hazard of Negative Binomial distribution provided by derived accessors. 334 // RealType chf(const negative_binomial_distribution<RealType, Policy>& dist) 335 // chf of Negative Binomial distribution provided by derived accessors. 336 337 template <class RealType, class Policy> pdf(const negative_binomial_distribution<RealType,Policy> & dist,const RealType & k)338 inline RealType pdf(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& k) 339 { // Probability Density/Mass Function. 340 BOOST_FPU_EXCEPTION_GUARD 341 342 static const char* function = "boost::math::pdf(const negative_binomial_distribution<%1%>&, %1%)"; 343 344 RealType r = dist.successes(); 345 RealType p = dist.success_fraction(); 346 RealType result = 0; 347 if(false == negative_binomial_detail::check_dist_and_k( 348 function, 349 r, 350 dist.success_fraction(), 351 k, 352 &result, Policy())) 353 { 354 return result; 355 } 356 357 result = (p/(r + k)) * ibeta_derivative(r, static_cast<RealType>(k+1), p, Policy()); 358 // Equivalent to: 359 // return exp(lgamma(r + k) - lgamma(r) - lgamma(k+1)) * pow(p, r) * pow((1-p), k); 360 return result; 361 } // negative_binomial_pdf 362 363 template <class RealType, class Policy> cdf(const negative_binomial_distribution<RealType,Policy> & dist,const RealType & k)364 inline RealType cdf(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& k) 365 { // Cumulative Distribution Function of Negative Binomial. 366 static const char* function = "boost::math::cdf(const negative_binomial_distribution<%1%>&, %1%)"; 367 using boost::math::ibeta; // Regularized incomplete beta function. 368 // k argument may be integral, signed, or unsigned, or floating point. 369 // If necessary, it has already been promoted from an integral type. 370 RealType p = dist.success_fraction(); 371 RealType r = dist.successes(); 372 // Error check: 373 RealType result = 0; 374 if(false == negative_binomial_detail::check_dist_and_k( 375 function, 376 r, 377 dist.success_fraction(), 378 k, 379 &result, Policy())) 380 { 381 return result; 382 } 383 384 RealType probability = ibeta(r, static_cast<RealType>(k+1), p, Policy()); 385 // Ip(r, k+1) = ibeta(r, k+1, p) 386 return probability; 387 } // cdf Cumulative Distribution Function Negative Binomial. 388 389 template <class RealType, class Policy> cdf(const complemented2_type<negative_binomial_distribution<RealType,Policy>,RealType> & c)390 inline RealType cdf(const complemented2_type<negative_binomial_distribution<RealType, Policy>, RealType>& c) 391 { // Complemented Cumulative Distribution Function Negative Binomial. 392 393 static const char* function = "boost::math::cdf(const negative_binomial_distribution<%1%>&, %1%)"; 394 using boost::math::ibetac; // Regularized incomplete beta function complement. 395 // k argument may be integral, signed, or unsigned, or floating point. 396 // If necessary, it has already been promoted from an integral type. 397 RealType const& k = c.param; 398 negative_binomial_distribution<RealType, Policy> const& dist = c.dist; 399 RealType p = dist.success_fraction(); 400 RealType r = dist.successes(); 401 // Error check: 402 RealType result = 0; 403 if(false == negative_binomial_detail::check_dist_and_k( 404 function, 405 r, 406 p, 407 k, 408 &result, Policy())) 409 { 410 return result; 411 } 412 // Calculate cdf negative binomial using the incomplete beta function. 413 // Use of ibeta here prevents cancellation errors in calculating 414 // 1-p if p is very small, perhaps smaller than machine epsilon. 415 // Ip(k+1, r) = ibetac(r, k+1, p) 416 // constrain_probability here? 417 RealType probability = ibetac(r, static_cast<RealType>(k+1), p, Policy()); 418 // Numerical errors might cause probability to be slightly outside the range < 0 or > 1. 419 // This might cause trouble downstream, so warn, possibly throw exception, but constrain to the limits. 420 return probability; 421 } // cdf Cumulative Distribution Function Negative Binomial. 422 423 template <class RealType, class Policy> quantile(const negative_binomial_distribution<RealType,Policy> & dist,const RealType & P)424 inline RealType quantile(const negative_binomial_distribution<RealType, Policy>& dist, const RealType& P) 425 { // Quantile, percentile/100 or Percent Point Negative Binomial function. 426 // Return the number of expected failures k for a given probability p. 427 428 // Inverse cumulative Distribution Function or Quantile (percentile / 100) of negative_binomial Probability. 429 // MAthCAD pnbinom return smallest k such that negative_binomial(k, n, p) >= probability. 430 // k argument may be integral, signed, or unsigned, or floating point. 431 // BUT Cephes/CodeCogs says: finds argument p (0 to 1) such that cdf(k, n, p) = y 432 static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)"; 433 BOOST_MATH_STD_USING // ADL of std functions. 434 435 RealType p = dist.success_fraction(); 436 RealType r = dist.successes(); 437 // Check dist and P. 438 RealType result = 0; 439 if(false == negative_binomial_detail::check_dist_and_prob 440 (function, r, p, P, &result, Policy())) 441 { 442 return result; 443 } 444 445 // Special cases. 446 if (P == 1) 447 { // Would need +infinity failures for total confidence. 448 result = policies::raise_overflow_error<RealType>( 449 function, 450 "Probability argument is 1, which implies infinite failures !", Policy()); 451 return result; 452 // usually means return +std::numeric_limits<RealType>::infinity(); 453 // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR 454 } 455 if (P == 0) 456 { // No failures are expected if P = 0. 457 return 0; // Total trials will be just dist.successes. 458 } 459 if (P <= pow(dist.success_fraction(), dist.successes())) 460 { // p <= pdf(dist, 0) == cdf(dist, 0) 461 return 0; 462 } 463 /* 464 // Calculate quantile of negative_binomial using the inverse incomplete beta function. 465 using boost::math::ibeta_invb; 466 return ibeta_invb(r, p, P, Policy()) - 1; // 467 */ 468 RealType guess = 0; 469 RealType factor = 5; 470 if(r * r * r * P * p > 0.005) 471 guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), P, RealType(1-P), Policy()); 472 473 if(guess < 10) 474 { 475 // 476 // Cornish-Fisher Negative binomial approximation not accurate in this area: 477 // 478 guess = (std::min)(RealType(r * 2), RealType(10)); 479 } 480 else 481 factor = (1-P < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); 482 BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); 483 // 484 // Max iterations permitted: 485 // 486 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 487 typedef typename Policy::discrete_quantile_type discrete_type; 488 return detail::inverse_discrete_quantile( 489 dist, 490 P, 491 1-P, 492 guess, 493 factor, 494 RealType(1), 495 discrete_type(), 496 max_iter); 497 } // RealType quantile(const negative_binomial_distribution dist, p) 498 499 template <class RealType, class Policy> quantile(const complemented2_type<negative_binomial_distribution<RealType,Policy>,RealType> & c)500 inline RealType quantile(const complemented2_type<negative_binomial_distribution<RealType, Policy>, RealType>& c) 501 { // Quantile or Percent Point Binomial function. 502 // Return the number of expected failures k for a given 503 // complement of the probability Q = 1 - P. 504 static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)"; 505 BOOST_MATH_STD_USING 506 507 // Error checks: 508 RealType Q = c.param; 509 const negative_binomial_distribution<RealType, Policy>& dist = c.dist; 510 RealType p = dist.success_fraction(); 511 RealType r = dist.successes(); 512 RealType result = 0; 513 if(false == negative_binomial_detail::check_dist_and_prob( 514 function, 515 r, 516 p, 517 Q, 518 &result, Policy())) 519 { 520 return result; 521 } 522 523 // Special cases: 524 // 525 if(Q == 1) 526 { // There may actually be no answer to this question, 527 // since the probability of zero failures may be non-zero, 528 return 0; // but zero is the best we can do: 529 } 530 if (-Q <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy())) 531 { // q <= cdf(complement(dist, 0)) == pdf(dist, 0) 532 return 0; // 533 } 534 if(Q == 0) 535 { // Probability 1 - Q == 1 so infinite failures to achieve certainty. 536 // Would need +infinity failures for total confidence. 537 result = policies::raise_overflow_error<RealType>( 538 function, 539 "Probability argument complement is 0, which implies infinite failures !", Policy()); 540 return result; 541 // usually means return +std::numeric_limits<RealType>::infinity(); 542 // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR 543 } 544 //return ibetac_invb(r, p, Q, Policy()) -1; 545 RealType guess = 0; 546 RealType factor = 5; 547 if(r * r * r * (1-Q) * p > 0.005) 548 guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), RealType(1-Q), Q, Policy()); 549 550 if(guess < 10) 551 { 552 // 553 // Cornish-Fisher Negative binomial approximation not accurate in this area: 554 // 555 guess = (std::min)(RealType(r * 2), RealType(10)); 556 } 557 else 558 factor = (Q < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); 559 BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); 560 // 561 // Max iterations permitted: 562 // 563 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 564 typedef typename Policy::discrete_quantile_type discrete_type; 565 return detail::inverse_discrete_quantile( 566 dist, 567 1-Q, 568 Q, 569 guess, 570 factor, 571 RealType(1), 572 discrete_type(), 573 max_iter); 574 } // quantile complement 575 576 } // namespace math 577 } // namespace boost 578 579 // This include must be at the end, *after* the accessors 580 // for this distribution have been defined, in order to 581 // keep compilers that support two-phase lookup happy. 582 #include <boost/math/distributions/detail/derived_accessors.hpp> 583 584 #if defined (BOOST_MSVC) 585 # pragma warning(pop) 586 #endif 587 588 #endif // BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP 589