1// Copyright John Maddock 2007.
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_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
7#define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
8
9#include <algorithm>
10
11namespace boost{ namespace math{ namespace detail{
12
13//
14// Functor for root finding algorithm:
15//
16template <class Dist>
17struct distribution_quantile_finder
18{
19 typedef typename Dist::value_type value_type;
20 typedef typename Dist::policy_type policy_type;
21
22 distribution_quantile_finder(const Dist d, value_type p, bool c)
23 : dist(d), target(p), comp(c) {}
24
25 value_type operator()(value_type const& x)
26 {
27 return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target);
28 }
29
30private:
31 Dist dist;
32 value_type target;
33 bool comp;
34};
35//
36// The purpose of adjust_bounds, is to toggle the last bit of the
37// range so that both ends round to the same integer, if possible.
38// If they do both round the same then we terminate the search
39// for the root *very* quickly when finding an integer result.
40// At the point that this function is called we know that "a" is
41// below the root and "b" above it, so this change can not result
42// in the root no longer being bracketed.
43//
44template <class Real, class Tol>
45void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){}
46
47template <class Real>
48void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */)
49{
50 BOOST_MATH_STD_USING
51 b -= tools::epsilon<Real>() * b;
52}
53
54template <class Real>
55void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */)
56{
57 BOOST_MATH_STD_USING
58 a += tools::epsilon<Real>() * a;
59}
60
61template <class Real>
62void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */)
63{
64 BOOST_MATH_STD_USING
65 a += tools::epsilon<Real>() * a;
66 b -= tools::epsilon<Real>() * b;
67}
68//
69// This is where all the work is done:
70//
71template <class Dist, class Tolerance>
72typename Dist::value_type
73 do_inverse_discrete_quantile(
74 const Dist& dist,
75 const typename Dist::value_type& p,
76 bool comp,
77 typename Dist::value_type guess,
78 const typename Dist::value_type& multiplier,
79 typename Dist::value_type adder,
80 const Tolerance& tol,
81 std::uintmax_t& max_iter)
82{
83 typedef typename Dist::value_type value_type;
84 typedef typename Dist::policy_type policy_type;
85
86 static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>";
87
88 BOOST_MATH_STD_USING
89
90 distribution_quantile_finder<Dist> f(dist, p, comp);
91 //
92 // Max bounds of the distribution:
93 //
94 value_type min_bound, max_bound;
95 boost::math::tie(min_bound, max_bound) = support(dist);
96
97 if(guess > max_bound)
98 guess = max_bound;
99 if(guess < min_bound)
100 guess = min_bound;
101
102 value_type fa = f(guess);
103 std::uintmax_t count = max_iter - 1;
104 value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used
105
106 if(fa == 0)
107 return guess;
108
109 //
110 // For small expected results, just use a linear search:
111 //
112 if(guess < 10)
113 {
114 b = a;
115 while((a < 10) && (fa * fb >= 0))
116 {
117 if(fb <= 0)
118 {
119 a = b;
120 b = a + 1;
121 if(b > max_bound)
122 b = max_bound;
123 fb = f(b);
124 --count;
125 if(fb == 0)
126 return b;
127 if(a == b)
128 return b; // can't go any higher!
129 }
130 else
131 {
132 b = a;
133 a = (std::max)(value_type(b - 1), value_type(0));
134 if(a < min_bound)
135 a = min_bound;
136 fa = f(a);
137 --count;
138 if(fa == 0)
139 return a;
140 if(a == b)
141 return a; // We can't go any lower than this!
142 }
143 }
144 }
145 //
146 // Try and bracket using a couple of additions first,
147 // we're assuming that "guess" is likely to be accurate
148 // to the nearest int or so:
149 //
150 else if((adder != 0) && (a + adder != a))
151 {
152 //
153 // If we're looking for a large result, then bump "adder" up
154 // by a bit to increase our chances of bracketing the root:
155 //
156 //adder = (std::max)(adder, 0.001f * guess);
157 if(fa < 0)
158 {
159 b = a + adder;
160 if(b > max_bound)
161 b = max_bound;
162 }
163 else
164 {
165 b = (std::max)(value_type(a - adder), value_type(0));
166 if(b < min_bound)
167 b = min_bound;
168 }
169 fb = f(b);
170 --count;
171 if(fb == 0)
172 return b;
173 if(count && (fa * fb >= 0))
174 {
175 //
176 // We didn't bracket the root, try
177 // once more:
178 //
179 a = b;
180 fa = fb;
181 if(fa < 0)
182 {
183 b = a + adder;
184 if(b > max_bound)
185 b = max_bound;
186 }
187 else
188 {
189 b = (std::max)(value_type(a - adder), value_type(0));
190 if(b < min_bound)
191 b = min_bound;
192 }
193 fb = f(b);
194 --count;
195 }
196 if(a > b)
197 {
198 using std::swap;
199 swap(a, b);
200 swap(fa, fb);
201 }
202 }
203 //
204 // If the root hasn't been bracketed yet, try again
205 // using the multiplier this time:
206 //
207 if((boost::math::sign)(fb) == (boost::math::sign)(fa))
208 {
209 if(fa < 0)
210 {
211 //
212 // Zero is to the right of x2, so walk upwards
213 // until we find it:
214 //
215 while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
216 {
217 if(count == 0)
218 return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); // LCOV_EXCL_LINE
219 a = b;
220 fa = fb;
221 b *= multiplier;
222 if(b > max_bound)
223 b = max_bound;
224 fb = f(b);
225 --count;
226 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
227 }
228 }
229 else
230 {
231 //
232 // Zero is to the left of a, so walk downwards
233 // until we find it:
234 //
235 while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))
236 {
237 if(fabs(a) < tools::min_value<value_type>())
238 {
239 // Escape route just in case the answer is zero!
240 max_iter -= count;
241 max_iter += 1;
242 return 0;
243 }
244 if(count == 0)
245 return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); // LCOV_EXCL_LINE
246 b = a;
247 fb = fa;
248 a /= multiplier;
249 if(a < min_bound)
250 a = min_bound;
251 fa = f(a);
252 --count;
253 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
254 }
255 }
256 }
257 max_iter -= count;
258 if(fa == 0)
259 return a;
260 if(fb == 0)
261 return b;
262 if(a == b)
263 return b; // Ran out of bounds trying to bracket - there is no answer!
264 //
265 // Adjust bounds so that if we're looking for an integer
266 // result, then both ends round the same way:
267 //
268 adjust_bounds(a, b, tol);
269 //
270 // We don't want zero or denorm lower bounds:
271 //
272 if(a < tools::min_value<value_type>())
273 a = tools::min_value<value_type>();
274 //
275 // Go ahead and find the root:
276 //
277 std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
278 max_iter += count;
279 if (max_iter >= policies::get_max_root_iterations<policy_type>())
280 {
281 return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
282 " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", r.first, policy_type()); // LCOV_EXCL_LINE
283 }
284 BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
285 return (r.first + r.second) / 2;
286}
287//
288// Some special routine for rounding up and down:
289// We want to check and see if we are very close to an integer, and if so test to see if
290// that integer is an exact root of the cdf. We do this because our root finder only
291// guarantees to find *a root*, and there can sometimes be many consecutive floating
292// point values which are all roots. This is especially true if the target probability
293// is very close 1.
294//
295template <class Dist>
296inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
297{
298 BOOST_MATH_STD_USING
299 typename Dist::value_type cc = ceil(result);
300 typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;
301 if(pp == p)
302 result = cc;
303 else
304 result = floor(result);
305 //
306 // Now find the smallest integer <= result for which we get an exact root:
307 //
308 while(result != 0)
309 {
310 cc = floor(float_prior(result));
311 if(cc < support(d).first)
312 break;
313 pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
314 if(c ? pp > p : pp < p)
315 break;
316 result = cc;
317 }
318
319 return result;
320}
321
322#ifdef _MSC_VER
323#pragma warning(push)
324#pragma warning(disable:4127)
325#endif
326
327template <class Dist>
328inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
329{
330 BOOST_MATH_STD_USING
331 typename Dist::value_type cc = floor(result);
332 typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;
333 if(pp == p)
334 result = cc;
335 else
336 result = ceil(result);
337 //
338 // Now find the largest integer >= result for which we get an exact root:
339 //
340 while(true)
341 {
342 cc = ceil(float_next(result));
343 if(cc > support(d).second)
344 break;
345 pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
346 if(c ? pp < p : pp > p)
347 break;
348 result = cc;
349 }
350
351 return result;
352}
353
354#ifdef _MSC_VER
355#pragma warning(pop)
356#endif
357//
358// Now finally are the public API functions.
359// There is one overload for each policy,
360// each one is responsible for selecting the correct
361// termination condition, and rounding the result
362// to an int where required.
363//
364template <class Dist>
365inline typename Dist::value_type
366 inverse_discrete_quantile(
367 const Dist& dist,
368 typename Dist::value_type p,
369 bool c,
370 const typename Dist::value_type& guess,
371 const typename Dist::value_type& multiplier,
372 const typename Dist::value_type& adder,
373 const policies::discrete_quantile<policies::real>&,
374 std::uintmax_t& max_iter)
375{
376 if(p > 0.5)
377 {
378 p = 1 - p;
379 c = !c;
380 }
381 typename Dist::value_type pp = c ? 1 - p : p;
382 if(pp <= pdf(dist, 0))
383 return 0;
384 return do_inverse_discrete_quantile(
385 dist,
386 p,
387 c,
388 guess,
389 multiplier,
390 adder,
391 tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),
392 max_iter);
393}
394
395template <class Dist>
396inline typename Dist::value_type
397 inverse_discrete_quantile(
398 const Dist& dist,
399 const typename Dist::value_type& p,
400 bool c,
401 const typename Dist::value_type& guess,
402 const typename Dist::value_type& multiplier,
403 const typename Dist::value_type& adder,
404 const policies::discrete_quantile<policies::integer_round_outwards>&,
405 std::uintmax_t& max_iter)
406{
407 typedef typename Dist::value_type value_type;
408 BOOST_MATH_STD_USING
409 typename Dist::value_type pp = c ? 1 - p : p;
410 if(pp <= pdf(dist, 0))
411 return 0;
412 //
413 // What happens next depends on whether we're looking for an
414 // upper or lower quantile:
415 //
416 if(pp < 0.5f)
417 return round_to_floor(dist, do_inverse_discrete_quantile(
418 dist,
419 p,
420 c,
421 (guess < 1 ? value_type(1) : (value_type)floor(guess)),
422 multiplier,
423 adder,
424 tools::equal_floor(),
425 max_iter), p, c);
426 // else:
427 return round_to_ceil(dist, do_inverse_discrete_quantile(
428 dist,
429 p,
430 c,
431 (value_type)ceil(guess),
432 multiplier,
433 adder,
434 tools::equal_ceil(),
435 max_iter), p, c);
436}
437
438template <class Dist>
439inline typename Dist::value_type
440 inverse_discrete_quantile(
441 const Dist& dist,
442 const typename Dist::value_type& p,
443 bool c,
444 const typename Dist::value_type& guess,
445 const typename Dist::value_type& multiplier,
446 const typename Dist::value_type& adder,
447 const policies::discrete_quantile<policies::integer_round_inwards>&,
448 std::uintmax_t& max_iter)
449{
450 typedef typename Dist::value_type value_type;
451 BOOST_MATH_STD_USING
452 typename Dist::value_type pp = c ? 1 - p : p;
453 if(pp <= pdf(dist, 0))
454 return 0;
455 //
456 // What happens next depends on whether we're looking for an
457 // upper or lower quantile:
458 //
459 if(pp < 0.5f)
460 return round_to_ceil(dist, do_inverse_discrete_quantile(
461 dist,
462 p,
463 c,
464 ceil(guess),
465 multiplier,
466 adder,
467 tools::equal_ceil(),
468 max_iter), p, c);
469 // else:
470 return round_to_floor(dist, do_inverse_discrete_quantile(
471 dist,
472 p,
473 c,
474 (guess < 1 ? value_type(1) : floor(guess)),
475 multiplier,
476 adder,
477 tools::equal_floor(),
478 max_iter), p, c);
479}
480
481template <class Dist>
482inline typename Dist::value_type
483 inverse_discrete_quantile(
484 const Dist& dist,
485 const typename Dist::value_type& p,
486 bool c,
487 const typename Dist::value_type& guess,
488 const typename Dist::value_type& multiplier,
489 const typename Dist::value_type& adder,
490 const policies::discrete_quantile<policies::integer_round_down>&,
491 std::uintmax_t& max_iter)
492{
493 typedef typename Dist::value_type value_type;
494 BOOST_MATH_STD_USING
495 typename Dist::value_type pp = c ? 1 - p : p;
496 if(pp <= pdf(dist, 0))
497 return 0;
498 return round_to_floor(dist, do_inverse_discrete_quantile(
499 dist,
500 p,
501 c,
502 (guess < 1 ? value_type(1) : floor(guess)),
503 multiplier,
504 adder,
505 tools::equal_floor(),
506 max_iter), p, c);
507}
508
509template <class Dist>
510inline typename Dist::value_type
511 inverse_discrete_quantile(
512 const Dist& dist,
513 const typename Dist::value_type& p,
514 bool c,
515 const typename Dist::value_type& guess,
516 const typename Dist::value_type& multiplier,
517 const typename Dist::value_type& adder,
518 const policies::discrete_quantile<policies::integer_round_up>&,
519 std::uintmax_t& max_iter)
520{
521 BOOST_MATH_STD_USING
522 typename Dist::value_type pp = c ? 1 - p : p;
523 if(pp <= pdf(dist, 0))
524 return 0;
525 return round_to_ceil(dist, do_inverse_discrete_quantile(
526 dist,
527 p,
528 c,
529 ceil(guess),
530 multiplier,
531 adder,
532 tools::equal_ceil(),
533 max_iter), p, c);
534}
535
536template <class Dist>
537inline typename Dist::value_type
538 inverse_discrete_quantile(
539 const Dist& dist,
540 const typename Dist::value_type& p,
541 bool c,
542 const typename Dist::value_type& guess,
543 const typename Dist::value_type& multiplier,
544 const typename Dist::value_type& adder,
545 const policies::discrete_quantile<policies::integer_round_nearest>&,
546 std::uintmax_t& max_iter)
547{
548 typedef typename Dist::value_type value_type;
549 BOOST_MATH_STD_USING
550 typename Dist::value_type pp = c ? 1 - p : p;
551 if(pp <= pdf(dist, 0))
552 return 0;
553 //
554 // Note that we adjust the guess to the nearest half-integer:
555 // this increase the chances that we will bracket the root
556 // with two results that both round to the same integer quickly.
557 //
558 return round_to_floor(dist, do_inverse_discrete_quantile(
559 dist,
560 p,
561 c,
562 (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
563 multiplier,
564 adder,
565 tools::equal_nearest_integer(),
566 max_iter) + 0.5f, p, c);
567}
568
569}}} // namespaces
570
571#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
572
573

source code of boost/libs/math/include/boost/math/distributions/detail/inv_discrete_quantile.hpp