ceph/src/boost/libs/math/include/boost/math/distributions/non_central_beta.hpp

   1 // boost\math\distributions\non_central_beta.hpp
   2
   3 // Copyright John Maddock 2008.
   4
   5 // Use, modification and distribution are subject to the
   6 // Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt
   8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
  11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
  12
  13 #include <boost/math/distributions/fwd.hpp>
  14 #include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q
  15 #include <boost/math/distributions/complement.hpp> // complements
  16 #include <boost/math/distributions/beta.hpp> // central distribution
  17 #include <boost/math/distributions/detail/generic_mode.hpp>
  18 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
  19 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
  20 #include <boost/math/tools/roots.hpp> // for root finding.
  21 #include <boost/math/tools/series.hpp>
  22
  23 namespace boost
  24 {
  25    namespace math
  26    {
  27
  28       template <class RealType, class Policy>
  29       class non_central_beta_distribution;
  30
  31       namespace detail{
  32
  33          template <class T, class Policy>
  34          T non_central_beta_p(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
  35          {
  36             BOOST_MATH_STD_USING
  37                using namespace boost::math;
  38             //
  39             // Variables come first:
  40             //
  41             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  42             T errtol = boost::math::policies::get_epsilon<T, Policy>();
  43             T l2 = lam / 2;
  44             //
  45             // k is the starting point for iteration, and is the
  46             // maximum of the poisson weighting term,
  47             // note that unlike other similar code, we do not set
  48             // k to zero, when l2 is small, as forward iteration
  49             // is unstable:
  50             //
  51             int k = itrunc(l2);
  52             if(k == 0)
  53                k = 1;
  54                // Starting Poisson weight:
  55             T pois = gamma_p_derivative(T(k+1), l2, pol);
  56             if(pois == 0)
  57                return init_val;
  58             // recurance term:
  59             T xterm;
  60             // Starting beta term:
  61             T beta = x < y
  62                ? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
  63                : detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
  64
  65             xterm *= y / (a + b + k - 1);
  66             T poisf(pois), betaf(beta), xtermf(xterm);
  67             T sum = init_val;
  68
  69             if((beta == 0) && (xterm == 0))
  70                return init_val;
  71
  72             //
  73             // Backwards recursion first, this is the stable
  74             // direction for recursion:
  75             //
  76             T last_term = 0;
  77             boost::uintmax_t count = k;
  78             for(int i = k; i >= 0; --i)
  79             {
  80                T term = beta * pois;
  81                sum += term;
  82                if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
  83                {
  84                   count = k - i;
  85                   break;
  86                }
  87                pois *= i / l2;
  88                beta += xterm;
  89                xterm *= (a + i - 1) / (x * (a + b + i - 2));
  90                last_term = term;
  91             }
  92             for(int i = k + 1; ; ++i)
  93             {
  94                poisf *= l2 / i;
  95                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
  96                betaf -= xtermf;
  97
  98                T term = poisf * betaf;
  99                sum += term;
 100                if((fabs(term/sum) < errtol) || (term == 0))
 101                {
 102                   break;
 103                }
 104                if(static_cast<boost::uintmax_t>(count + i - k) > max_iter)
 105                {
 106                   return policies::raise_evaluation_error(
 107                      "cdf(non_central_beta_distribution<%1%>, %1%)",
 108                      "Series did not converge, closest value was %1%", sum, pol);
 109                }
 110             }
 111             return sum;
 112          }
 113
 114          template <class T, class Policy>
 115          T non_central_beta_q(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
 116          {
 117             BOOST_MATH_STD_USING
 118                using namespace boost::math;
 119             //
 120             // Variables come first:
 121             //
 122             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 123             T errtol = boost::math::policies::get_epsilon<T, Policy>();
 124             T l2 = lam / 2;
 125             //
 126             // k is the starting point for iteration, and is the
 127             // maximum of the poisson weighting term:
 128             //
 129             int k = itrunc(l2);
 130             T pois;
 131             if(k <= 30)
 132             {
 133                //
 134                // Might as well start at 0 since we'll likely have this number of terms anyway:
 135                //
 136                if(a + b > 1)
 137                   k = 0;
 138                else if(k == 0)
 139                   k = 1;
 140             }
 141             if(k == 0)
 142             {
 143                // Starting Poisson weight:
 144                pois = exp(-l2);
 145             }
 146             else
 147             {
 148                // Starting Poisson weight:
 149                pois = gamma_p_derivative(T(k+1), l2, pol);
 150             }
 151             if(pois == 0)
 152                return init_val;
 153             // recurance term:
 154             T xterm;
 155             // Starting beta term:
 156             T beta = x < y
 157                ? detail::ibeta_imp(T(a + k), b, x, pol, true, true, &xterm)
 158                : detail::ibeta_imp(b, T(a + k), y, pol, false, true, &xterm);
 159
 160             xterm *= y / (a + b + k - 1);
 161             T poisf(pois), betaf(beta), xtermf(xterm);
 162             T sum = init_val;
 163             if((beta == 0) && (xterm == 0))
 164                return init_val;
 165             //
 166             // Forwards recursion first, this is the stable
 167             // direction for recursion, and the location
 168             // of the bulk of the sum:
 169             //
 170             T last_term = 0;
 171             boost::uintmax_t count = 0;
 172             for(int i = k + 1; ; ++i)
 173             {
 174                poisf *= l2 / i;
 175                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
 176                betaf += xtermf;
 177
 178                T term = poisf * betaf;
 179                sum += term;
 180                if((fabs(term/sum) < errtol) && (last_term >= term))
 181                {
 182                   count = i - k;
 183                   break;
 184                }
 185                if(static_cast<boost::uintmax_t>(i - k) > max_iter)
 186                {
 187                   return policies::raise_evaluation_error(
 188                      "cdf(non_central_beta_distribution<%1%>, %1%)",
 189                      "Series did not converge, closest value was %1%", sum, pol);
 190                }
 191                last_term = term;
 192             }
 193             for(int i = k; i >= 0; --i)
 194             {
 195                T term = beta * pois;
 196                sum += term;
 197                if(fabs(term/sum) < errtol)
 198                {
 199                   break;
 200                }
 201                if(static_cast<boost::uintmax_t>(count + k - i) > max_iter)
 202                {
 203                   return policies::raise_evaluation_error(
 204                      "cdf(non_central_beta_distribution<%1%>, %1%)",
 205                      "Series did not converge, closest value was %1%", sum, pol);
 206                }
 207                pois *= i / l2;
 208                beta -= xterm;
 209                xterm *= (a + i - 1) / (x * (a + b + i - 2));
 210             }
 211             return sum;
 212          }
 213
 214          template <class RealType, class Policy>
 215          inline RealType non_central_beta_cdf(RealType x, RealType y, RealType a, RealType b, RealType l, bool invert, const Policy&)
 216          {
 217             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 218             typedef typename policies::normalise<
 219                Policy,
 220                policies::promote_float<false>,
 221                policies::promote_double<false>,
 222                policies::discrete_quantile<>,
 223                policies::assert_undefined<> >::type forwarding_policy;
 224
 225             BOOST_MATH_STD_USING
 226
 227             if(x == 0)
 228                return invert ? 1.0f : 0.0f;
 229             if(y == 0)
 230                return invert ? 0.0f : 1.0f;
 231             value_type result;
 232             value_type c = a + b + l / 2;
 233             value_type cross = 1 - (b / c) * (1 + l / (2 * c * c));
 234             if(l == 0)
 235                result = cdf(boost::math::beta_distribution<RealType, Policy>(a, b), x);
 236             else if(x > cross)
 237             {
 238                // Complement is the smaller of the two:
 239                result = detail::non_central_beta_q(
 240                   static_cast<value_type>(a),
 241                   static_cast<value_type>(b),
 242                   static_cast<value_type>(l),
 243                   static_cast<value_type>(x),
 244                   static_cast<value_type>(y),
 245                   forwarding_policy(),
 246                   static_cast<value_type>(invert ? 0 : -1));
 247                invert = !invert;
 248             }
 249             else
 250             {
 251                result = detail::non_central_beta_p(
 252                   static_cast<value_type>(a),
 253                   static_cast<value_type>(b),
 254                   static_cast<value_type>(l),
 255                   static_cast<value_type>(x),
 256                   static_cast<value_type>(y),
 257                   forwarding_policy(),
 258                   static_cast<value_type>(invert ? -1 : 0));
 259             }
 260             if(invert)
 261                result = -result;
 262             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 263                result,
 264                "boost::math::non_central_beta_cdf<%1%>(%1%, %1%, %1%)");
 265          }
 266
 267          template <class T, class Policy>
 268          struct nc_beta_quantile_functor
 269          {
 270             nc_beta_quantile_functor(const non_central_beta_distribution<T,Policy>& d, T t, bool c)
 271                : dist(d), target(t), comp(c) {}
 272
 273             T operator()(const T& x)
 274             {
 275                return comp ?
 276                   T(target - cdf(complement(dist, x)))
 277                   : T(cdf(dist, x) - target);
 278             }
 279
 280          private:
 281             non_central_beta_distribution<T,Policy> dist;
 282             T target;
 283             bool comp;
 284          };
 285
 286          //
 287          // This is more or less a copy of bracket_and_solve_root, but
 288          // modified to search only the interval [0,1] using similar
 289          // heuristics.
 290          //
 291          template <class F, class T, class Tol, class Policy>
 292          std::pair<T, T> bracket_and_solve_root_01(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
 293          {
 294             BOOST_MATH_STD_USING
 295                static const char* function = "boost::math::tools::bracket_and_solve_root_01<%1%>";
 296             //
 297             // Set up inital brackets:
 298             //
 299             T a = guess;
 300             T b = a;
 301             T fa = f(a);
 302             T fb = fa;
 303             //
 304             // Set up invocation count:
 305             //
 306             boost::uintmax_t count = max_iter - 1;
 307
 308             if((fa < 0) == (guess < 0 ? !rising : rising))
 309             {
 310                //
 311                // Zero is to the right of b, so walk upwards
 312                // until we find it:
 313                //
 314                while((boost::math::sign)(fb) == (boost::math::sign)(fa))
 315                {
 316                   if(count == 0)
 317                   {
 318                      b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol);
 319                      return std::make_pair(a, b);
 320                   }
 321                   //
 322                   // Heuristic: every 20 iterations we double the growth factor in case the
 323                   // initial guess was *really* bad !
 324                   //
 325                   if((max_iter - count) % 20 == 0)
 326                      factor *= 2;
 327                   //
 328                   // Now go ahead and move are guess by "factor",
 329                   // we do this by reducing 1-guess by factor:
 330                   //
 331                   a = b;
 332                   fa = fb;
 333                   b = 1 - ((1 - b) / factor);
 334                   fb = f(b);
 335                   --count;
 336                   BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 337                }
 338             }
 339             else
 340             {
 341                //
 342                // Zero is to the left of a, so walk downwards
 343                // until we find it:
 344                //
 345                while((boost::math::sign)(fb) == (boost::math::sign)(fa))
 346                {
 347                   if(fabs(a) < tools::min_value<T>())
 348                   {
 349                      // Escape route just in case the answer is zero!
 350                      max_iter -= count;
 351                      max_iter += 1;
 352                      return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
 353                   }
 354                   if(count == 0)
 355                   {
 356                      a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol);
 357                      return std::make_pair(a, b);
 358                   }
 359                   //
 360                   // Heuristic: every 20 iterations we double the growth factor in case the
 361                   // initial guess was *really* bad !
 362                   //
 363                   if((max_iter - count) % 20 == 0)
 364                      factor *= 2;
 365                   //
 366                   // Now go ahead and move are guess by "factor":
 367                   //
 368                   b = a;
 369                   fb = fa;
 370                   a /= factor;
 371                   fa = f(a);
 372                   --count;
 373                   BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 374                }
 375             }
 376             max_iter -= count;
 377             max_iter += 1;
 378             std::pair<T, T> r = toms748_solve(
 379                f,
 380                (a < 0 ? b : a),
 381                (a < 0 ? a : b),
 382                (a < 0 ? fb : fa),
 383                (a < 0 ? fa : fb),
 384                tol,
 385                count,
 386                pol);
 387             max_iter += count;
 388             BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
 389             return r;
 390          }
 391
 392          template <class RealType, class Policy>
 393          RealType nc_beta_quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
 394          {
 395             static const char* function = "quantile(non_central_beta_distribution<%1%>, %1%)";
 396             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 397             typedef typename policies::normalise<
 398                Policy,
 399                policies::promote_float<false>,
 400                policies::promote_double<false>,
 401                policies::discrete_quantile<>,
 402                policies::assert_undefined<> >::type forwarding_policy;
 403
 404             value_type a = dist.alpha();
 405             value_type b = dist.beta();
 406             value_type l = dist.non_centrality();
 407             value_type r;
 408             if(!beta_detail::check_alpha(
 409                function,
 410                a, &r, Policy())
 411                ||
 412             !beta_detail::check_beta(
 413                function,
 414                b, &r, Policy())
 415                ||
 416             !detail::check_non_centrality(
 417                function,
 418                l,
 419                &r,
 420                Policy())
 421                ||
 422             !detail::check_probability(
 423                function,
 424                static_cast<value_type>(p),
 425                &r,
 426                Policy()))
 427                   return (RealType)r;
 428             //
 429             // Special cases first:
 430             //
 431             if(p == 0)
 432                return comp
 433                ? 1.0f
 434                : 0.0f;
 435             if(p == 1)
 436                return !comp
 437                ? 1.0f
 438                : 0.0f;
 439
 440             value_type c = a + b + l / 2;
 441             value_type mean = 1 - (b / c) * (1 + l / (2 * c * c));
 442             /*
 443             //
 444             // Calculate a normal approximation to the quantile,
 445             // uses mean and variance approximations from:
 446             // Algorithm AS 310:
 447             // Computing the Non-Central Beta Distribution Function
 448             // R. Chattamvelli; R. Shanmugam
 449             // Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156.
 450             //
 451             // Unfortunately, when this is wrong it tends to be *very*
 452             // wrong, so it's disabled for now, even though it often
 453             // gets the initial guess quite close.  Probably we could
 454             // do much better by factoring in the skewness if only
 455             // we could calculate it....
 456             //
 457             value_type delta = l / 2;
 458             value_type delta2 = delta * delta;
 459             value_type delta3 = delta * delta2;
 460             value_type delta4 = delta2 * delta2;
 461             value_type G = c * (c + 1) + delta;
 462             value_type alpha = a + b;
 463             value_type alpha2 = alpha * alpha;
 464             value_type eta = (2 * alpha + 1) * (2 * alpha + 1) + 1;
 465             value_type H = 3 * alpha2 + 5 * alpha + 2;
 466             value_type F = alpha2 * (alpha + 1) + H * delta
 467                + (2 * alpha + 4) * delta2 + delta3;
 468             value_type P = (3 * alpha + 1) * (9 * alpha + 17)
 469                + 2 * alpha * (3 * alpha + 2) * (3 * alpha + 4) + 15;
 470             value_type Q = 54 * alpha2 + 162 * alpha + 130;
 471             value_type R = 6 * (6 * alpha + 11);
 472             value_type D = delta
 473                * (H * H + 2 * P * delta + Q * delta2 + R * delta3 + 9 * delta4);
 474             value_type variance = (b / G)
 475                * (1 + delta * (l * l + 3 * l + eta) / (G * G))
 476                - (b * b / F) * (1 + D / (F * F));
 477             value_type sd = sqrt(variance);
 478
 479             value_type guess = comp
 480                ? quantile(complement(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p))
 481                : quantile(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p);
 482
 483             if(guess >= 1)
 484                guess = mean;
 485             if(guess <= tools::min_value<value_type>())
 486                guess = mean;
 487             */
 488             value_type guess = mean;
 489             detail::nc_beta_quantile_functor<value_type, Policy>
 490                f(non_central_beta_distribution<value_type, Policy>(a, b, l), p, comp);
 491             tools::eps_tolerance<value_type> tol(policies::digits<RealType, Policy>());
 492             boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 493
 494             std::pair<value_type, value_type> ir
 495                = bracket_and_solve_root_01(
 496                   f, guess, value_type(2.5), true, tol,
 497                   max_iter, Policy());
 498             value_type result = ir.first + (ir.second - ir.first) / 2;
 499
 500             if(max_iter >= policies::get_max_root_iterations<Policy>())
 501             {
 502                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
 503                   " either there is no answer to quantile of the non central beta distribution"
 504                   " or the answer is infinite.  Current best guess is %1%",
 505                   policies::checked_narrowing_cast<RealType, forwarding_policy>(
 506                      result,
 507                      function), Policy());
 508             }
 509             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 510                result,
 511                function);
 512          }
 513
 514          template <class T, class Policy>
 515          T non_central_beta_pdf(T a, T b, T lam, T x, T y, const Policy& pol)
 516          {
 517             BOOST_MATH_STD_USING
 518             //
 519             // Special cases:
 520             //
 521             if((x == 0) || (y == 0))
 522                return 0;
 523             //
 524             // Variables come first:
 525             //
 526             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 527             T errtol = boost::math::policies::get_epsilon<T, Policy>();
 528             T l2 = lam / 2;
 529             //
 530             // k is the starting point for iteration, and is the
 531             // maximum of the poisson weighting term:
 532             //
 533             int k = itrunc(l2);
 534             // Starting Poisson weight:
 535             T pois = gamma_p_derivative(T(k+1), l2, pol);
 536             // Starting beta term:
 537             T beta = x < y ?
 538                ibeta_derivative(a + k, b, x, pol)
 539                : ibeta_derivative(b, a + k, y, pol);
 540             T sum = 0;
 541             T poisf(pois);
 542             T betaf(beta);
 543
 544             //
 545             // Stable backwards recursion first:
 546             //
 547             boost::uintmax_t count = k;
 548             for(int i = k; i >= 0; --i)
 549             {
 550                T term = beta * pois;
 551                sum += term;
 552                if((fabs(term/sum) < errtol) || (term == 0))
 553                {
 554                   count = k - i;
 555                   break;
 556                }
 557                pois *= i / l2;
 558                beta *= (a + i - 1) / (x * (a + i + b - 1));
 559             }
 560             for(int i = k + 1; ; ++i)
 561             {
 562                poisf *= l2 / i;
 563                betaf *= x * (a + b + i - 1) / (a + i - 1);
 564
 565                T term = poisf * betaf;
 566                sum += term;
 567                if((fabs(term/sum) < errtol) || (term == 0))
 568                {
 569                   break;
 570                }
 571                if(static_cast<boost::uintmax_t>(count + i - k) > max_iter)
 572                {
 573                   return policies::raise_evaluation_error(
 574                      "pdf(non_central_beta_distribution<%1%>, %1%)",
 575                      "Series did not converge, closest value was %1%", sum, pol);
 576                }
 577             }
 578             return sum;
 579          }
 580
 581          template <class RealType, class Policy>
 582          RealType nc_beta_pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
 583          {
 584             BOOST_MATH_STD_USING
 585             static const char* function = "pdf(non_central_beta_distribution<%1%>, %1%)";
 586             typedef typename policies::evaluation<RealType, Policy>::type value_type;
 587             typedef typename policies::normalise<
 588                Policy,
 589                policies::promote_float<false>,
 590                policies::promote_double<false>,
 591                policies::discrete_quantile<>,
 592                policies::assert_undefined<> >::type forwarding_policy;
 593
 594             value_type a = dist.alpha();
 595             value_type b = dist.beta();
 596             value_type l = dist.non_centrality();
 597             value_type r;
 598             if(!beta_detail::check_alpha(
 599                function,
 600                a, &r, Policy())
 601                ||
 602             !beta_detail::check_beta(
 603                function,
 604                b, &r, Policy())
 605                ||
 606             !detail::check_non_centrality(
 607                function,
 608                l,
 609                &r,
 610                Policy())
 611                ||
 612             !beta_detail::check_x(
 613                function,
 614                static_cast<value_type>(x),
 615                &r,
 616                Policy()))
 617                   return (RealType)r;
 618
 619             if(l == 0)
 620                return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
 621             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 622                non_central_beta_pdf(a, b, l, static_cast<value_type>(x), value_type(1 - static_cast<value_type>(x)), forwarding_policy()),
 623                "function");
 624          }
 625
 626          template <class T>
 627          struct hypergeometric_2F2_sum
 628          {
 629             typedef T result_type;
 630             hypergeometric_2F2_sum(T a1_, T a2_, T b1_, T b2_, T z_) : a1(a1_), a2(a2_), b1(b1_), b2(b2_), z(z_), term(1), k(0) {}
 631             T operator()()
 632             {
 633                T result = term;
 634                term *= a1 * a2 / (b1 * b2);
 635                a1 += 1;
 636                a2 += 1;
 637                b1 += 1;
 638                b2 += 1;
 639                k += 1;
 640                term /= k;
 641                term *= z;
 642                return result;
 643             }
 644             T a1, a2, b1, b2, z, term, k;
 645          };
 646
 647          template <class T, class Policy>
 648          T hypergeometric_2F2(T a1, T a2, T b1, T b2, T z, const Policy& pol)
 649          {
 650             typedef typename policies::evaluation<T, Policy>::type value_type;
 651
 652             const char* function = "boost::math::detail::hypergeometric_2F2<%1%>(%1%,%1%,%1%,%1%,%1%)";
 653
 654             hypergeometric_2F2_sum<value_type> s(a1, a2, b1, b2, z);
 655             boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 656 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 657             value_type zero = 0;
 658             value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter, zero);
 659 #else
 660             value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter);
 661 #endif
 662             policies::check_series_iterations<T>(function, max_iter, pol);
 663             return policies::checked_narrowing_cast<T, Policy>(result, function);
 664          }
 665
 666       } // namespace detail
 667
 668       template <class RealType = double, class Policy = policies::policy<> >
 669       class non_central_beta_distribution
 670       {
 671       public:
 672          typedef RealType value_type;
 673          typedef Policy policy_type;
 674
 675          non_central_beta_distribution(RealType a_, RealType b_, RealType lambda) : a(a_), b(b_), ncp(lambda)
 676          {
 677             const char* function = "boost::math::non_central_beta_distribution<%1%>::non_central_beta_distribution(%1%,%1%)";
 678             RealType r;
 679             beta_detail::check_alpha(
 680                function,
 681                a, &r, Policy());
 682             beta_detail::check_beta(
 683                function,
 684                b, &r, Policy());
 685             detail::check_non_centrality(
 686                function,
 687                lambda,
 688                &r,
 689                Policy());
 690          } // non_central_beta_distribution constructor.
 691
 692          RealType alpha() const
 693          { // Private data getter function.
 694             return a;
 695          }
 696          RealType beta() const
 697          { // Private data getter function.
 698             return b;
 699          }
 700          RealType non_centrality() const
 701          { // Private data getter function.
 702             return ncp;
 703          }
 704       private:
 705          // Data member, initialized by constructor.
 706          RealType a;   // alpha.
 707          RealType b;   // beta.
 708          RealType ncp; // non-centrality parameter
 709       }; // template <class RealType, class Policy> class non_central_beta_distribution
 710
 711       typedef non_central_beta_distribution<double> non_central_beta; // Reserved name of type double.
 712
 713       // Non-member functions to give properties of the distribution.
 714
 715       template <class RealType, class Policy>
 716       inline const std::pair<RealType, RealType> range(const non_central_beta_distribution<RealType, Policy>& /* dist */)
 717       { // Range of permissible values for random variable k.
 718          using boost::math::tools::max_value;
 719          return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
 720       }
 721
 722       template <class RealType, class Policy>
 723       inline const std::pair<RealType, RealType> support(const non_central_beta_distribution<RealType, Policy>& /* dist */)
 724       { // Range of supported values for random variable k.
 725          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 726          using boost::math::tools::max_value;
 727          return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
 728       }
 729
 730       template <class RealType, class Policy>
 731       inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist)
 732       { // mode.
 733          static const char* function = "mode(non_central_beta_distribution<%1%> const&)";
 734
 735          RealType a = dist.alpha();
 736          RealType b = dist.beta();
 737          RealType l = dist.non_centrality();
 738          RealType r;
 739          if(!beta_detail::check_alpha(
 740                function,
 741                a, &r, Policy())
 742                ||
 743             !beta_detail::check_beta(
 744                function,
 745                b, &r, Policy())
 746                ||
 747             !detail::check_non_centrality(
 748                function,
 749                l,
 750                &r,
 751                Policy()))
 752                   return (RealType)r;
 753          RealType c = a + b + l / 2;
 754          RealType mean = 1 - (b / c) * (1 + l / (2 * c * c));
 755          return detail::generic_find_mode_01(
 756             dist,
 757             mean,
 758             function);
 759       }
 760
 761       //
 762       // We don't have the necessary information to implement
 763       // these at present.  These are just disabled for now,
 764       // prototypes retained so we can fill in the blanks
 765       // later:
 766       //
 767       template <class RealType, class Policy>
 768       inline RealType mean(const non_central_beta_distribution<RealType, Policy>& dist)
 769       {
 770          BOOST_MATH_STD_USING
 771          RealType a = dist.alpha();
 772          RealType b = dist.beta();
 773          RealType d = dist.non_centrality();
 774          RealType apb = a + b;
 775          return exp(-d / 2) * a * detail::hypergeometric_2F2<RealType, Policy>(1 + a, apb, a, 1 + apb, d / 2, Policy()) / apb;
 776       } // mean
 777
 778       template <class RealType, class Policy>
 779       inline RealType variance(const non_central_beta_distribution<RealType, Policy>& dist)
 780       {
 781          //
 782          // Relative error of this function may be arbitarily large... absolute
 783          // error will be small however... that's the best we can do for now.
 784          //
 785          BOOST_MATH_STD_USING
 786          RealType a = dist.alpha();
 787          RealType b = dist.beta();
 788          RealType d = dist.non_centrality();
 789          RealType apb = a + b;
 790          RealType result = detail::hypergeometric_2F2(RealType(1 + a), apb, a, RealType(1 + apb), RealType(d / 2), Policy());
 791          result *= result * -exp(-d) * a * a / (apb * apb);
 792          result += exp(-d / 2) * a * (1 + a) * detail::hypergeometric_2F2(RealType(2 + a), apb, a, RealType(2 + apb), RealType(d / 2), Policy()) / (apb * (1 + apb));
 793          return result;
 794       }
 795
 796       // RealType standard_deviation(const non_central_beta_distribution<RealType, Policy>& dist)
 797       // standard_deviation provided by derived accessors.
 798       template <class RealType, class Policy>
 799       inline RealType skewness(const non_central_beta_distribution<RealType, Policy>& /*dist*/)
 800       { // skewness = sqrt(l).
 801          const char* function = "boost::math::non_central_beta_distribution<%1%>::skewness()";
 802          typedef typename Policy::assert_undefined_type assert_type;
 803          BOOST_STATIC_ASSERT(assert_type::value == 0);
 804
 805          return policies::raise_evaluation_error<RealType>(
 806             function,
 807             "This function is not yet implemented, the only sensible result is %1%.",
 808             std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
 809       }
 810
 811       template <class RealType, class Policy>
 812       inline RealType kurtosis_excess(const non_central_beta_distribution<RealType, Policy>& /*dist*/)
 813       {
 814          const char* function = "boost::math::non_central_beta_distribution<%1%>::kurtosis_excess()";
 815          typedef typename Policy::assert_undefined_type assert_type;
 816          BOOST_STATIC_ASSERT(assert_type::value == 0);
 817
 818          return policies::raise_evaluation_error<RealType>(
 819             function,
 820             "This function is not yet implemented, the only sensible result is %1%.",
 821             std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
 822       } // kurtosis_excess
 823
 824       template <class RealType, class Policy>
 825       inline RealType kurtosis(const non_central_beta_distribution<RealType, Policy>& dist)
 826       {
 827          return kurtosis_excess(dist) + 3;
 828       }
 829
 830       template <class RealType, class Policy>
 831       inline RealType pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
 832       { // Probability Density/Mass Function.
 833          return detail::nc_beta_pdf(dist, x);
 834       } // pdf
 835
 836       template <class RealType, class Policy>
 837       RealType cdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
 838       {
 839          const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)";
 840             RealType a = dist.alpha();
 841             RealType b = dist.beta();
 842             RealType l = dist.non_centrality();
 843             RealType r;
 844             if(!beta_detail::check_alpha(
 845                function,
 846                a, &r, Policy())
 847                ||
 848             !beta_detail::check_beta(
 849                function,
 850                b, &r, Policy())
 851                ||
 852             !detail::check_non_centrality(
 853                function,
 854                l,
 855                &r,
 856                Policy())
 857                ||
 858             !beta_detail::check_x(
 859                function,
 860                x,
 861                &r,
 862                Policy()))
 863                   return (RealType)r;
 864
 865          if(l == 0)
 866             return cdf(beta_distribution<RealType, Policy>(a, b), x);
 867
 868          return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, false, Policy());
 869       } // cdf
 870
 871       template <class RealType, class Policy>
 872       RealType cdf(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c)
 873       { // Complemented Cumulative Distribution Function
 874          const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)";
 875          non_central_beta_distribution<RealType, Policy> const& dist = c.dist;
 876             RealType a = dist.alpha();
 877             RealType b = dist.beta();
 878             RealType l = dist.non_centrality();
 879             RealType x = c.param;
 880             RealType r;
 881             if(!beta_detail::check_alpha(
 882                function,
 883                a, &r, Policy())
 884                ||
 885             !beta_detail::check_beta(
 886                function,
 887                b, &r, Policy())
 888                ||
 889             !detail::check_non_centrality(
 890                function,
 891                l,
 892                &r,
 893                Policy())
 894                ||
 895             !beta_detail::check_x(
 896                function,
 897                x,
 898                &r,
 899                Policy()))
 900                   return (RealType)r;
 901
 902          if(l == 0)
 903             return cdf(complement(beta_distribution<RealType, Policy>(a, b), x));
 904
 905          return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, true, Policy());
 906       } // ccdf
 907
 908       template <class RealType, class Policy>
 909       inline RealType quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p)
 910       { // Quantile (or Percent Point) function.
 911          return detail::nc_beta_quantile(dist, p, false);
 912       } // quantile
 913
 914       template <class RealType, class Policy>
 915       inline RealType quantile(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c)
 916       { // Quantile (or Percent Point) function.
 917          return detail::nc_beta_quantile(c.dist, c.param, true);
 918       } // quantile complement.
 919
 920    } // namespace math
 921 } // namespace boost
 922
 923 // This include must be at the end, *after* the accessors
 924 // for this distribution have been defined, in order to
 925 // keep compilers that support two-phase lookup happy.
 926 #include <boost/math/distributions/detail/derived_accessors.hpp>
 927
 928 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
 929