Ticket #9183: math-discrete-distro.patch

File math-discrete-distro.patch, 15.7 KB (added by John Maddock, 9 years ago)

  • math/distributions/binomial.hpp

     
    196196         }
    197197
    198198      template <class RealType, class Policy>
    199       RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q)
     199      RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp)
    200200      { // Quantile or Percent Point Binomial function.
    201201        // Return the number of expected successes k,
    202202        // for a given probability p.
    … …  
    264264        boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
    265265        return detail::inverse_discrete_quantile(
    266266            dist,
    267             p,
    268             q,
     267            comp ? q : p,
     268            comp,
    269269            guess,
    270270            factor,
    271271            RealType(1),
    … …  
    653653      template <class RealType, class Policy>
    654654      inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
    655655      {
    656          return binomial_detail::quantile_imp(dist, p, RealType(1-p));
     656         return binomial_detail::quantile_imp(dist, p, RealType(1-p), false);
    657657      } // quantile
    658658
    659659      template <class RealType, class Policy>
    660660      RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
    661661      {
    662          return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param);
     662         return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true);
    663663      } // quantile
    664664
    665665      template <class RealType, class Policy>

  • math/distributions/detail/inv_discrete_quantile.hpp

     
    1919   typedef typename Dist::value_type value_type;
    2020   typedef typename Dist::policy_type policy_type;
    2121
    22    distribution_quantile_finder(const Dist d, value_type p, value_type q)
    23       : dist(d), target(p < q ? p : q), comp(p < q ? false : true) {}
     22   distribution_quantile_finder(const Dist d, value_type p, bool c)
     23      : dist(d), target(p), comp(c) {}
    2424
    2525   value_type operator()(value_type const& x)
    2626   {
    … …  
    7373   do_inverse_discrete_quantile(
    7474      const Dist& dist,
    7575      const typename Dist::value_type& p,
    76       const typename Dist::value_type& q,
     76      bool comp,
    7777      typename Dist::value_type guess,
    7878      const typename Dist::value_type& multiplier,
    7979      typename Dist::value_type adder,
    … …  
    8787
    8888   BOOST_MATH_STD_USING
    8989
    90    distribution_quantile_finder<Dist> f(dist, p, q);
     90   distribution_quantile_finder<Dist> f(dist, p, comp);
    9191   //
    9292   // Max bounds of the distribution:
    9393   //
    … …  
    280280   return (r.first + r.second) / 2;
    281281}
    282282//
     283// Some special routine for rounding up and down:
     284// We want to check and see if we are very close to an integer, and if so test to see if
     285// that integer is an exact root of the cdf.  We do this because our root finder only
     286// guarantees to find *a root*, and there can sometimes be many consecutive floating
     287// point values which are all roots.  This is especially true if the target probability
     288// is very close 1.
     289//
     290template <class Dist>
     291inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
     292{
     293   BOOST_MATH_STD_USING
     294   typename Dist::value_type cc = ceil(result);
     295   typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;
     296   if(pp == p)
     297      result = cc;
     298   else
     299      result = floor(result);
     300   //
     301   // Now find the smallest integer <= result for which we get an exact root:
     302   //
     303   while(result != 0)
     304   {
     305      cc = result - 1;
     306      if(cc < support(d).first)
     307         break;
     308      typename Dist::value_type pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
     309      if(pp == p)
     310         result = cc;
     311      else if(c ? pp > p : pp < p)
     312         break;
     313      result -= 1;
     314   }
     315
     316   return result;
     317}
     318template <class Dist>
     319inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)
     320{
     321   BOOST_MATH_STD_USING
     322   typename Dist::value_type cc = floor(result);
     323   typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;
     324   if(pp == p)
     325      result = cc;
     326   else
     327      result = ceil(result);
     328   //
     329   // Now find the largest integer >= result for which we get an exact root:
     330   //
     331   while(true)
     332   {
     333      cc = result + 1;
     334      if(cc > support(d).second)
     335         break;
     336      typename Dist::value_type pp = c ? cdf(complement(d, cc)) : cdf(d, cc);
     337      if(pp == p)
     338         result = cc;
     339      else if(c ? pp < p : pp > p)
     340         break;
     341      result += 1;
     342   }
     343
     344   return result;
     345}
     346//
    283347// Now finally are the public API functions.
    284348// There is one overload for each policy,
    285349// each one is responsible for selecting the correct
    … …  
    290354inline typename Dist::value_type
    291355   inverse_discrete_quantile(
    292356      const Dist& dist,
    293       const typename Dist::value_type& p,
    294       const typename Dist::value_type& q,
     357      typename Dist::value_type p,
     358      bool c,
    295359      const typename Dist::value_type& guess,
    296360      const typename Dist::value_type& multiplier,
    297361      const typename Dist::value_type& adder,
    298362      const policies::discrete_quantile<policies::real>&,
    299363      boost::uintmax_t& max_iter)
    300364{
    301    if(p <= pdf(dist, 0))
     365   if(p > 0.5)
     366   {
     367      p = 1 - p;
     368      c = !c;
     369   }
     370   typename Dist::value_type pp = c ? 1 - p : p;
     371   if(pp <= pdf(dist, 0))
    302372      return 0;
    303373   return do_inverse_discrete_quantile(
    304374      dist,
    305375      p,
    306       q,
     376      c,
    307377      guess,
    308378      multiplier,
    309379      adder,
    … …  
    316386   inverse_discrete_quantile(
    317387      const Dist& dist,
    318388      const typename Dist::value_type& p,
    319       const typename Dist::value_type& q,
     389      bool c,
    320390      const typename Dist::value_type& guess,
    321391      const typename Dist::value_type& multiplier,
    322392      const typename Dist::value_type& adder,
    … …  
    325395{
    326396   typedef typename Dist::value_type value_type;
    327397   BOOST_MATH_STD_USING
    328    if(p <= pdf(dist, 0))
     398   typename Dist::value_type pp = c ? 1 - p : p;
     399   if(pp <= pdf(dist, 0))
    329400      return 0;
    330401   //
    331402   // What happens next depends on whether we're looking for an
    332403   // upper or lower quantile:
    333404   //
    334    if(p < 0.5f)
    335       return floor(do_inverse_discrete_quantile(
     405   if(pp < 0.5f)
     406      return round_to_floor(dist, do_inverse_discrete_quantile(
    336407         dist,
    337408         p,
    338          q,
     409         c,
    339410         (guess < 1 ? value_type(1) : (value_type)floor(guess)),
    340411         multiplier,
    341412         adder,
    342413         tools::equal_floor(),
    343          max_iter));
     414         max_iter), p, c);
    344415   // else:
    345    return ceil(do_inverse_discrete_quantile(
     416   return round_to_ceil(dist, do_inverse_discrete_quantile(
    346417      dist,
    347418      p,
    348       q,
     419      c,
    349420      (value_type)ceil(guess),
    350421      multiplier,
    351422      adder,
    352423      tools::equal_ceil(),
    353       max_iter));
     424      max_iter), p, c);
    354425}
    355426
    356427template <class Dist>
    … …  
    358429   inverse_discrete_quantile(
    359430      const Dist& dist,
    360431      const typename Dist::value_type& p,
    361       const typename Dist::value_type& q,
     432      bool c,
    362433      const typename Dist::value_type& guess,
    363434      const typename Dist::value_type& multiplier,
    364435      const typename Dist::value_type& adder,
    … …  
    367438{
    368439   typedef typename Dist::value_type value_type;
    369440   BOOST_MATH_STD_USING
    370    if(p <= pdf(dist, 0))
     441   typename Dist::value_type pp = c ? 1 - p : p;
     442   if(pp <= pdf(dist, 0))
    371443      return 0;
    372444   //
    373445   // What happens next depends on whether we're looking for an
    374446   // upper or lower quantile:
    375447   //
    376    if(p < 0.5f)
    377       return ceil(do_inverse_discrete_quantile(
     448   if(pp < 0.5f)
     449      return round_to_ceil(dist, do_inverse_discrete_quantile(
    378450         dist,
    379451         p,
    380          q,
     452         c,
    381453         ceil(guess),
    382454         multiplier,
    383455         adder,
    384456         tools::equal_ceil(),
    385          max_iter));
     457         max_iter), p, c);
    386458   // else:
    387    return floor(do_inverse_discrete_quantile(
     459   return round_to_floor(dist, do_inverse_discrete_quantile(
    388460      dist,
    389461      p,
    390       q,
     462      c,
    391463      (guess < 1 ? value_type(1) : floor(guess)),
    392464      multiplier,
    393465      adder,
    394466      tools::equal_floor(),
    395       max_iter));
     467      max_iter), p, c);
    396468}
    397469
    398470template <class Dist>
    … …  
    400472   inverse_discrete_quantile(
    401473      const Dist& dist,
    402474      const typename Dist::value_type& p,
    403       const typename Dist::value_type& q,
     475      bool c,
    404476      const typename Dist::value_type& guess,
    405477      const typename Dist::value_type& multiplier,
    406478      const typename Dist::value_type& adder,
    … …  
    409481{
    410482   typedef typename Dist::value_type value_type;
    411483   BOOST_MATH_STD_USING
    412    if(p <= pdf(dist, 0))
     484   typename Dist::value_type pp = c ? 1 - p : p;
     485   if(pp <= pdf(dist, 0))
    413486      return 0;
    414    return floor(do_inverse_discrete_quantile(
     487   return round_to_floor(dist, do_inverse_discrete_quantile(
    415488      dist,
    416489      p,
    417       q,
     490      c,
    418491      (guess < 1 ? value_type(1) : floor(guess)),
    419492      multiplier,
    420493      adder,
    421494      tools::equal_floor(),
    422       max_iter));
     495      max_iter), p, c);
    423496}
    424497
    425498template <class Dist>
    … …  
    427500   inverse_discrete_quantile(
    428501      const Dist& dist,
    429502      const typename Dist::value_type& p,
    430       const typename Dist::value_type& q,
     503      bool c,
    431504      const typename Dist::value_type& guess,
    432505      const typename Dist::value_type& multiplier,
    433506      const typename Dist::value_type& adder,
    … …  
    435508      boost::uintmax_t& max_iter)
    436509{
    437510   BOOST_MATH_STD_USING
    438    if(p <= pdf(dist, 0))
     511   typename Dist::value_type pp = c ? 1 - p : p;
     512   if(pp <= pdf(dist, 0))
    439513      return 0;
    440    return ceil(do_inverse_discrete_quantile(
     514   return round_to_ceil(dist, do_inverse_discrete_quantile(
    441515      dist,
    442516      p,
    443       q,
     517      c,
    444518      ceil(guess),
    445519      multiplier,
    446520      adder,
    447521      tools::equal_ceil(),
    448       max_iter));
     522      max_iter), p, c);
    449523}
    450524
    451525template <class Dist>
    … …  
    453527   inverse_discrete_quantile(
    454528      const Dist& dist,
    455529      const typename Dist::value_type& p,
    456       const typename Dist::value_type& q,
     530      bool c,
    457531      const typename Dist::value_type& guess,
    458532      const typename Dist::value_type& multiplier,
    459533      const typename Dist::value_type& adder,
    … …  
    462536{
    463537   typedef typename Dist::value_type value_type;
    464538   BOOST_MATH_STD_USING
    465    if(p <= pdf(dist, 0))
     539   typename Dist::value_type pp = c ? 1 - p : p;
     540   if(pp <= pdf(dist, 0))
    466541      return 0;
    467542   //
    468543   // Note that we adjust the guess to the nearest half-integer:
    469544   // this increase the chances that we will bracket the root
    470545   // with two results that both round to the same integer quickly.
    471546   //
    472    return floor(do_inverse_discrete_quantile(
     547   return round_to_floor(dist, do_inverse_discrete_quantile(
    473548      dist,
    474549      p,
    475       q,
     550      c,
    476551      (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f),
    477552      multiplier,
    478553      adder,
    479554      tools::equal_nearest_integer(),
    480       max_iter) + 0.5f);
     555      max_iter) + 0.5f, p, c);
    481556}
    482557
    483558}}} // namespaces
    484559
    485560#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE
    486561
    487 

  • math/distributions/negative_binomial.hpp

     
    488488      return detail::inverse_discrete_quantile(
    489489         dist,
    490490         P,
    491          1-P,
     491         false,
    492492         guess,
    493493         factor,
    494494         RealType(1),
    … …  
    564564       typedef typename Policy::discrete_quantile_type discrete_type;
    565565       return detail::inverse_discrete_quantile(
    566566          dist,
    567           1-Q,
    568567          Q,
     568          true,
    569569          guess,
    570570          factor,
    571571          RealType(1),

  • math/distributions/poisson.hpp

     
    5252{
    5353  namespace math
    5454  {
    55      namespace detail{
    56       template <class Dist>
    57       inline typename Dist::value_type
    58          inverse_discrete_quantile(
    59             const Dist& dist,
    60             const typename Dist::value_type& p,
    61             const typename Dist::value_type& guess,
    62             const typename Dist::value_type& multiplier,
    63             const typename Dist::value_type& adder,
    64             const policies::discrete_quantile<policies::integer_round_nearest>&,
    65             boost::uintmax_t& max_iter);
    66       template <class Dist>
    67       inline typename Dist::value_type
    68          inverse_discrete_quantile(
    69             const Dist& dist,
    70             const typename Dist::value_type& p,
    71             const typename Dist::value_type& guess,
    72             const typename Dist::value_type& multiplier,
    73             const typename Dist::value_type& adder,
    74             const policies::discrete_quantile<policies::integer_round_up>&,
    75             boost::uintmax_t& max_iter);
    76       template <class Dist>
    77       inline typename Dist::value_type
    78          inverse_discrete_quantile(
    79             const Dist& dist,
    80             const typename Dist::value_type& p,
    81             const typename Dist::value_type& guess,
    82             const typename Dist::value_type& multiplier,
    83             const typename Dist::value_type& adder,
    84             const policies::discrete_quantile<policies::integer_round_down>&,
    85             boost::uintmax_t& max_iter);
    86       template <class Dist>
    87       inline typename Dist::value_type
    88          inverse_discrete_quantile(
    89             const Dist& dist,
    90             const typename Dist::value_type& p,
    91             const typename Dist::value_type& guess,
    92             const typename Dist::value_type& multiplier,
    93             const typename Dist::value_type& adder,
    94             const policies::discrete_quantile<policies::integer_round_outwards>&,
    95             boost::uintmax_t& max_iter);
    96       template <class Dist>
    97       inline typename Dist::value_type
    98          inverse_discrete_quantile(
    99             const Dist& dist,
    100             const typename Dist::value_type& p,
    101             const typename Dist::value_type& guess,
    102             const typename Dist::value_type& multiplier,
    103             const typename Dist::value_type& adder,
    104             const policies::discrete_quantile<policies::integer_round_inwards>&,
    105             boost::uintmax_t& max_iter);
    106       template <class Dist>
    107       inline typename Dist::value_type
    108          inverse_discrete_quantile(
    109             const Dist& dist,
    110             const typename Dist::value_type& p,
    111             const typename Dist::value_type& guess,
    112             const typename Dist::value_type& multiplier,
    113             const typename Dist::value_type& adder,
    114             const policies::discrete_quantile<policies::real>&,
    115             boost::uintmax_t& max_iter);
    116      }
    11755    namespace poisson_detail
    11856    {
    11957      // Common error checking routines for Poisson distribution functions.
    … …  
    496434      return detail::inverse_discrete_quantile(
    497435         dist,
    498436         p,
    499          1-p,
     437         false,
    500438         guess,
    501439         factor,
    502440         RealType(1),
    … …  
    565503
    566504      return detail::inverse_discrete_quantile(
    567505         dist,
    568          1-q,
    569506         q,
     507         true,
    570508         guess,
    571509         factor,
    572510         RealType(1),