boost/math/distributions/detail/generic_mode.hpp

   1 // Copyright John Maddock 2008.
   2
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
   9 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
  10
  11 #include <boost/math/tools/minima.hpp> // function minimization for mode
  12 #include <boost/math/policies/error_handling.hpp>
  13 #include <boost/math/distributions/fwd.hpp>
  14
  15 namespace boost{ namespace math{ namespace detail{
  16
  17 template <class Dist>
  18 struct pdf_minimizer
  19 {
  20    pdf_minimizer(const Dist& d)
  21       : dist(d) {}
  22
  23    typename Dist::value_type operator()(const typename Dist::value_type& x)
  24    {
  25       return -pdf(dist, x);
  26    }
  27 private:
  28    Dist dist;
  29 };
  30
  31 template <class Dist>
  32 typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0)
  33 {
  34    BOOST_MATH_STD_USING
  35    typedef typename Dist::value_type value_type;
  36    typedef typename Dist::policy_type policy_type;
  37    //
  38    // Need to begin by bracketing the maxima of the PDF:
  39    //
  40    value_type maxval;
  41    value_type upper_bound = guess;
  42    value_type lower_bound;
  43    value_type v = pdf(dist, guess);
  44    if(v == 0)
  45    {
  46       //
  47       // Oops we don't know how to handle this, or even in which
  48       // direction we should move in, treat as an evaluation error:
  49       //
  50       return policies::raise_evaluation_error(
  51          function,
  52          "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
  53    }
  54    do
  55    {
  56       maxval = v;
  57       if(step != 0)
  58          upper_bound += step;
  59       else
  60          upper_bound *= 2;
  61       v = pdf(dist, upper_bound);
  62    }while(maxval < v);
  63
  64    lower_bound = upper_bound;
  65    do
  66    {
  67       maxval = v;
  68       if(step != 0)
  69          lower_bound -= step;
  70       else
  71          lower_bound /= 2;
  72       v = pdf(dist, lower_bound);
  73    }while(maxval < v);
  74
  75    boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
  76
  77    value_type result = tools::brent_find_minima(
  78       pdf_minimizer<Dist>(dist),
  79       lower_bound,
  80       upper_bound,
  81       policies::digits<value_type, policy_type>(),
  82       max_iter).first;
  83    if(max_iter >= policies::get_max_root_iterations<policy_type>())
  84    {
  85       return policies::raise_evaluation_error<value_type>(
  86          function,
  87          "Unable to locate solution in a reasonable time:"
  88          " either there is no answer to the mode of the distribution"
  89          " or the answer is infinite.  Current best guess is %1%", result, policy_type());
  90    }
  91    return result;
  92 }
  93 //
  94 // As above,but confined to the interval [0,1]:
  95 //
  96 template <class Dist>
  97 typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
  98 {
  99    BOOST_MATH_STD_USING
 100    typedef typename Dist::value_type value_type;
 101    typedef typename Dist::policy_type policy_type;
 102    //
 103    // Need to begin by bracketing the maxima of the PDF:
 104    //
 105    value_type maxval;
 106    value_type upper_bound = guess;
 107    value_type lower_bound;
 108    value_type v = pdf(dist, guess);
 109    do
 110    {
 111       maxval = v;
 112       upper_bound = 1 - (1 - upper_bound) / 2;
 113       if(upper_bound == 1)
 114          return 1;
 115       v = pdf(dist, upper_bound);
 116    }while(maxval < v);
 117
 118    lower_bound = upper_bound;
 119    do
 120    {
 121       maxval = v;
 122       lower_bound /= 2;
 123       if(lower_bound < tools::min_value<value_type>())
 124          return 0;
 125       v = pdf(dist, lower_bound);
 126    }while(maxval < v);
 127
 128    boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
 129
 130    value_type result = tools::brent_find_minima(
 131       pdf_minimizer<Dist>(dist),
 132       lower_bound,
 133       upper_bound,
 134       policies::digits<value_type, policy_type>(),
 135       max_iter).first;
 136    if(max_iter >= policies::get_max_root_iterations<policy_type>())
 137    {
 138       return policies::raise_evaluation_error<value_type>(
 139          function,
 140          "Unable to locate solution in a reasonable time:"
 141          " either there is no answer to the mode of the distribution"
 142          " or the answer is infinite.  Current best guess is %1%", result, policy_type());
 143    }
 144    return result;
 145 }
 146
 147 }}} // namespaces
 148
 149 #endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP