]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
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> | |
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> | |
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> | |
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> | |
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 | ||
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. | |
153 | RealType success_fraction() const | |
154 | { // Probability of success as fraction in range 0 to 1. | |
155 | return m_p; | |
156 | } | |
157 | RealType successes() const | |
158 | { // Total number of successes r. | |
159 | return m_r; | |
160 | } | |
161 | ||
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 | ||
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 | ||
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 | ||
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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 | if(p == 0) | |
464 | { // Would need +infinity failures for total confidence. | |
465 | result = policies::raise_overflow_error<RealType>( | |
466 | function, | |
467 | "Success fraction is 0, which implies infinite failures !", Policy()); | |
468 | return result; | |
469 | // usually means return +std::numeric_limits<RealType>::infinity(); | |
470 | // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR | |
471 | } | |
472 | /* | |
473 | // Calculate quantile of negative_binomial using the inverse incomplete beta function. | |
474 | using boost::math::ibeta_invb; | |
475 | return ibeta_invb(r, p, P, Policy()) - 1; // | |
476 | */ | |
477 | RealType guess = 0; | |
478 | RealType factor = 5; | |
479 | if(r * r * r * P * p > 0.005) | |
480 | guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), P, RealType(1-P), Policy()); | |
481 | ||
482 | if(guess < 10) | |
483 | { | |
484 | // | |
485 | // Cornish-Fisher Negative binomial approximation not accurate in this area: | |
486 | // | |
487 | guess = (std::min)(RealType(r * 2), RealType(10)); | |
488 | } | |
489 | else | |
490 | factor = (1-P < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); | |
491 | BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); | |
492 | // | |
493 | // Max iterations permitted: | |
494 | // | |
495 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
496 | typedef typename Policy::discrete_quantile_type discrete_type; | |
497 | return detail::inverse_discrete_quantile( | |
498 | dist, | |
499 | P, | |
500 | false, | |
501 | guess, | |
502 | factor, | |
503 | RealType(1), | |
504 | discrete_type(), | |
505 | max_iter); | |
506 | } // RealType quantile(const negative_binomial_distribution dist, p) | |
507 | ||
508 | template <class RealType, class Policy> | |
509 | inline RealType quantile(const complemented2_type<negative_binomial_distribution<RealType, Policy>, RealType>& c) | |
510 | { // Quantile or Percent Point Binomial function. | |
511 | // Return the number of expected failures k for a given | |
512 | // complement of the probability Q = 1 - P. | |
513 | static const char* function = "boost::math::quantile(const negative_binomial_distribution<%1%>&, %1%)"; | |
514 | BOOST_MATH_STD_USING | |
515 | ||
516 | // Error checks: | |
517 | RealType Q = c.param; | |
518 | const negative_binomial_distribution<RealType, Policy>& dist = c.dist; | |
519 | RealType p = dist.success_fraction(); | |
520 | RealType r = dist.successes(); | |
521 | RealType result = 0; | |
522 | if(false == negative_binomial_detail::check_dist_and_prob( | |
523 | function, | |
524 | r, | |
525 | p, | |
526 | Q, | |
527 | &result, Policy())) | |
528 | { | |
529 | return result; | |
530 | } | |
531 | ||
532 | // Special cases: | |
533 | // | |
534 | if(Q == 1) | |
535 | { // There may actually be no answer to this question, | |
536 | // since the probability of zero failures may be non-zero, | |
537 | return 0; // but zero is the best we can do: | |
538 | } | |
539 | if(Q == 0) | |
540 | { // Probability 1 - Q == 1 so infinite failures to achieve certainty. | |
541 | // Would need +infinity failures for total confidence. | |
542 | result = policies::raise_overflow_error<RealType>( | |
543 | function, | |
544 | "Probability argument complement is 0, which implies infinite failures !", Policy()); | |
545 | return result; | |
546 | // usually means return +std::numeric_limits<RealType>::infinity(); | |
547 | // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR | |
548 | } | |
549 | if (-Q <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy())) | |
550 | { // q <= cdf(complement(dist, 0)) == pdf(dist, 0) | |
551 | return 0; // | |
552 | } | |
553 | if(p == 0) | |
554 | { // Success fraction is 0 so infinite failures to achieve certainty. | |
555 | // Would need +infinity failures for total confidence. | |
556 | result = policies::raise_overflow_error<RealType>( | |
557 | function, | |
558 | "Success fraction is 0, which implies infinite failures !", Policy()); | |
559 | return result; | |
560 | // usually means return +std::numeric_limits<RealType>::infinity(); | |
561 | // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR | |
562 | } | |
563 | //return ibetac_invb(r, p, Q, Policy()) -1; | |
564 | RealType guess = 0; | |
565 | RealType factor = 5; | |
566 | if(r * r * r * (1-Q) * p > 0.005) | |
567 | guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), RealType(1-Q), Q, Policy()); | |
568 | ||
569 | if(guess < 10) | |
570 | { | |
571 | // | |
572 | // Cornish-Fisher Negative binomial approximation not accurate in this area: | |
573 | // | |
574 | guess = (std::min)(RealType(r * 2), RealType(10)); | |
575 | } | |
576 | else | |
577 | factor = (Q < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f); | |
578 | BOOST_MATH_INSTRUMENT_CODE("guess = " << guess); | |
579 | // | |
580 | // Max iterations permitted: | |
581 | // | |
582 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
583 | typedef typename Policy::discrete_quantile_type discrete_type; | |
584 | return detail::inverse_discrete_quantile( | |
585 | dist, | |
586 | Q, | |
587 | true, | |
588 | guess, | |
589 | factor, | |
590 | RealType(1), | |
591 | discrete_type(), | |
592 | max_iter); | |
593 | } // quantile complement | |
594 | ||
595 | } // namespace math | |
596 | } // namespace boost | |
597 | ||
598 | // This include must be at the end, *after* the accessors | |
599 | // for this distribution have been defined, in order to | |
600 | // keep compilers that support two-phase lookup happy. | |
601 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
602 | ||
603 | #if defined (BOOST_MSVC) | |
604 | # pragma warning(pop) | |
605 | #endif | |
606 | ||
607 | #endif // BOOST_MATH_SPECIAL_NEGATIVE_BINOMIAL_HPP |