math/distributions/non_central_beta.hpp

0001 // boost\math\distributions\non_central_beta.hpp
0002
0003 // Copyright John Maddock 2008.
0004
0005 // Use, modification and distribution are subject to the
0006 // Boost Software License, Version 1.0.
0007 // (See accompanying file LICENSE_1_0.txt
0008 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
0011 #define BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
0012
0013 #include <boost/math/distributions/fwd.hpp>
0014 #include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q
0015 #include <boost/math/distributions/complement.hpp> // complements
0016 #include <boost/math/distributions/beta.hpp> // central distribution
0017 #include <boost/math/distributions/detail/generic_mode.hpp>
0018 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0019 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
0020 #include <boost/math/special_functions/trunc.hpp>
0021 #include <boost/math/tools/roots.hpp> // for root finding.
0022 #include <boost/math/tools/series.hpp>
0023
0024 namespace boost
0025 {
0026    namespace math
0027    {
0028
0029       template <class RealType, class Policy>
0030       class non_central_beta_distribution;
0031
0032       namespace detail{
0033
0034          template <class T, class Policy>
0035          T non_central_beta_p(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
0036          {
0037             BOOST_MATH_STD_USING
0038                using namespace boost::math;
0039             //
0040             // Variables come first:
0041             //
0042             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0043             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0044             T l2 = lam / 2;
0045             //
0046             // k is the starting point for iteration, and is the
0047             // maximum of the poisson weighting term,
0048             // note that unlike other similar code, we do not set
0049             // k to zero, when l2 is small, as forward iteration
0050             // is unstable:
0051             //
0052             long long k = lltrunc(l2);
0053             if(k == 0)
0054                k = 1;
0055                // Starting Poisson weight:
0056             T pois = gamma_p_derivative(T(k+1), l2, pol);
0057             if(pois == 0)
0058                return init_val;
0059             // recurance term:
0060             T xterm;
0061             // Starting beta term:
0062             T beta = x < y
0063                ? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
0064                : detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
0065
0066             xterm *= y / (a + b + k - 1);
0067             T poisf(pois), betaf(beta), xtermf(xterm);
0068             T sum = init_val;
0069
0070             if((beta == 0) && (xterm == 0))
0071                return init_val;
0072
0073             //
0074             // Backwards recursion first, this is the stable
0075             // direction for recursion:
0076             //
0077             T last_term = 0;
0078             std::uintmax_t count = k;
0079             for(auto i = k; i >= 0; --i)
0080             {
0081                T term = beta * pois;
0082                sum += term;
0083                if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
0084                {
0085                   count = k - i;
0086                   break;
0087                }
0088                pois *= i / l2;
0089                beta += xterm;
0090
0091                if (a + b + i != 2)
0092                {
0093                   xterm *= (a + i - 1) / (x * (a + b + i - 2));
0094                }
0095
0096                last_term = term;
0097             }
0098             for(auto i = k + 1; ; ++i)
0099             {
0100                poisf *= l2 / i;
0101                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
0102                betaf -= xtermf;
0103
0104                T term = poisf * betaf;
0105                sum += term;
0106                if((fabs(term/sum) < errtol) || (term == 0))
0107                {
0108                   break;
0109                }
0110                if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
0111                {
0112                   return policies::raise_evaluation_error(
0113                      "cdf(non_central_beta_distribution<%1%>, %1%)",
0114                      "Series did not converge, closest value was %1%", sum, pol);
0115                }
0116             }
0117             return sum;
0118          }
0119
0120          template <class T, class Policy>
0121          T non_central_beta_q(T a, T b, T lam, T x, T y, const Policy& pol, T init_val = 0)
0122          {
0123             BOOST_MATH_STD_USING
0124                using namespace boost::math;
0125             //
0126             // Variables come first:
0127             //
0128             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0129             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0130             T l2 = lam / 2;
0131             //
0132             // k is the starting point for iteration, and is the
0133             // maximum of the poisson weighting term:
0134             //
0135             long long k = lltrunc(l2);
0136             T pois;
0137             if(k <= 30)
0138             {
0139                //
0140                // Might as well start at 0 since we'll likely have this number of terms anyway:
0141                //
0142                if(a + b > 1)
0143                   k = 0;
0144                else if(k == 0)
0145                   k = 1;
0146             }
0147             if(k == 0)
0148             {
0149                // Starting Poisson weight:
0150                pois = exp(-l2);
0151             }
0152             else
0153             {
0154                // Starting Poisson weight:
0155                pois = gamma_p_derivative(T(k+1), l2, pol);
0156             }
0157             if(pois == 0)
0158                return init_val;
0159             // recurance term:
0160             T xterm;
0161             // Starting beta term:
0162             T beta = x < y
0163                ? detail::ibeta_imp(T(a + k), b, x, pol, true, true, &xterm)
0164                : detail::ibeta_imp(b, T(a + k), y, pol, false, true, &xterm);
0165
0166             xterm *= y / (a + b + k - 1);
0167             T poisf(pois), betaf(beta), xtermf(xterm);
0168             T sum = init_val;
0169             if((beta == 0) && (xterm == 0))
0170                return init_val;
0171             //
0172             // Forwards recursion first, this is the stable
0173             // direction for recursion, and the location
0174             // of the bulk of the sum:
0175             //
0176             T last_term = 0;
0177             std::uintmax_t count = 0;
0178             for(auto i = k + 1; ; ++i)
0179             {
0180                poisf *= l2 / i;
0181                xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
0182                betaf += xtermf;
0183
0184                T term = poisf * betaf;
0185                sum += term;
0186                if((fabs(term/sum) < errtol) && (last_term >= term))
0187                {
0188                   count = i - k;
0189                   break;
0190                }
0191                if(static_cast<std::uintmax_t>(i - k) > max_iter)
0192                {
0193                   return policies::raise_evaluation_error(
0194                      "cdf(non_central_beta_distribution<%1%>, %1%)",
0195                      "Series did not converge, closest value was %1%", sum, pol);
0196                }
0197                last_term = term;
0198             }
0199             for(auto i = k; i >= 0; --i)
0200             {
0201                T term = beta * pois;
0202                sum += term;
0203                if(fabs(term/sum) < errtol)
0204                {
0205                   break;
0206                }
0207                if(static_cast<std::uintmax_t>(count + k - i) > max_iter)
0208                {
0209                   return policies::raise_evaluation_error(
0210                      "cdf(non_central_beta_distribution<%1%>, %1%)",
0211                      "Series did not converge, closest value was %1%", sum, pol);
0212                }
0213                pois *= i / l2;
0214                beta -= xterm;
0215                xterm *= (a + i - 1) / (x * (a + b + i - 2));
0216             }
0217             return sum;
0218          }
0219
0220          template <class RealType, class Policy>
0221          inline RealType non_central_beta_cdf(RealType x, RealType y, RealType a, RealType b, RealType l, bool invert, const Policy&)
0222          {
0223             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0224             typedef typename policies::normalise<
0225                Policy,
0226                policies::promote_float<false>,
0227                policies::promote_double<false>,
0228                policies::discrete_quantile<>,
0229                policies::assert_undefined<> >::type forwarding_policy;
0230
0231             BOOST_MATH_STD_USING
0232
0233             if(x == 0)
0234                return invert ? 1.0f : 0.0f;
0235             if(y == 0)
0236                return invert ? 0.0f : 1.0f;
0237             value_type result;
0238             value_type c = a + b + l / 2;
0239             value_type cross = 1 - (b / c) * (1 + l / (2 * c * c));
0240             if(l == 0)
0241                result = cdf(boost::math::beta_distribution<RealType, Policy>(a, b), x);
0242             else if(x > cross)
0243             {
0244                // Complement is the smaller of the two:
0245                result = detail::non_central_beta_q(
0246                   static_cast<value_type>(a),
0247                   static_cast<value_type>(b),
0248                   static_cast<value_type>(l),
0249                   static_cast<value_type>(x),
0250                   static_cast<value_type>(y),
0251                   forwarding_policy(),
0252                   static_cast<value_type>(invert ? 0 : -1));
0253                invert = !invert;
0254             }
0255             else
0256             {
0257                result = detail::non_central_beta_p(
0258                   static_cast<value_type>(a),
0259                   static_cast<value_type>(b),
0260                   static_cast<value_type>(l),
0261                   static_cast<value_type>(x),
0262                   static_cast<value_type>(y),
0263                   forwarding_policy(),
0264                   static_cast<value_type>(invert ? -1 : 0));
0265             }
0266             if(invert)
0267                result = -result;
0268             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0269                result,
0270                "boost::math::non_central_beta_cdf<%1%>(%1%, %1%, %1%)");
0271          }
0272
0273          template <class T, class Policy>
0274          struct nc_beta_quantile_functor
0275          {
0276             nc_beta_quantile_functor(const non_central_beta_distribution<T,Policy>& d, T t, bool c)
0277                : dist(d), target(t), comp(c) {}
0278
0279             T operator()(const T& x)
0280             {
0281                return comp ?
0282                   T(target - cdf(complement(dist, x)))
0283                   : T(cdf(dist, x) - target);
0284             }
0285
0286          private:
0287             non_central_beta_distribution<T,Policy> dist;
0288             T target;
0289             bool comp;
0290          };
0291
0292          //
0293          // This is more or less a copy of bracket_and_solve_root, but
0294          // modified to search only the interval [0,1] using similar
0295          // heuristics.
0296          //
0297          template <class F, class T, class Tol, class Policy>
0298          std::pair<T, T> bracket_and_solve_root_01(F f, const T& guess, T factor, bool rising, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
0299          {
0300             BOOST_MATH_STD_USING
0301                static const char* function = "boost::math::tools::bracket_and_solve_root_01<%1%>";
0302             //
0303             // Set up initial brackets:
0304             //
0305             T a = guess;
0306             T b = a;
0307             T fa = f(a);
0308             T fb = fa;
0309             //
0310             // Set up invocation count:
0311             //
0312             std::uintmax_t count = max_iter - 1;
0313
0314             if((fa < 0) == (guess < 0 ? !rising : rising))
0315             {
0316                //
0317                // Zero is to the right of b, so walk upwards
0318                // until we find it:
0319                //
0320                while((boost::math::sign)(fb) == (boost::math::sign)(fa))
0321                {
0322                   if(count == 0)
0323                   {
0324                      b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol);
0325                      return std::make_pair(a, b);
0326                   }
0327                   //
0328                   // Heuristic: every 20 iterations we double the growth factor in case the
0329                   // initial guess was *really* bad !
0330                   //
0331                   if((max_iter - count) % 20 == 0)
0332                      factor *= 2;
0333                   //
0334                   // Now go ahead and move are guess by "factor",
0335                   // we do this by reducing 1-guess by factor:
0336                   //
0337                   a = b;
0338                   fa = fb;
0339                   b = 1 - ((1 - b) / factor);
0340                   fb = f(b);
0341                   --count;
0342                   BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0343                }
0344             }
0345             else
0346             {
0347                //
0348                // Zero is to the left of a, so walk downwards
0349                // until we find it:
0350                //
0351                while((boost::math::sign)(fb) == (boost::math::sign)(fa))
0352                {
0353                   if(fabs(a) < tools::min_value<T>())
0354                   {
0355                      // Escape route just in case the answer is zero!
0356                      max_iter -= count;
0357                      max_iter += 1;
0358                      return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
0359                   }
0360                   if(count == 0)
0361                   {
0362                      a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol);
0363                      return std::make_pair(a, b);
0364                   }
0365                   //
0366                   // Heuristic: every 20 iterations we double the growth factor in case the
0367                   // initial guess was *really* bad !
0368                   //
0369                   if((max_iter - count) % 20 == 0)
0370                      factor *= 2;
0371                   //
0372                   // Now go ahead and move are guess by "factor":
0373                   //
0374                   b = a;
0375                   fb = fa;
0376                   a /= factor;
0377                   fa = f(a);
0378                   --count;
0379                   BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0380                }
0381             }
0382             max_iter -= count;
0383             max_iter += 1;
0384             std::pair<T, T> r = toms748_solve(
0385                f,
0386                (a < 0 ? b : a),
0387                (a < 0 ? a : b),
0388                (a < 0 ? fb : fa),
0389                (a < 0 ? fa : fb),
0390                tol,
0391                count,
0392                pol);
0393             max_iter += count;
0394             BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
0395             return r;
0396          }
0397
0398          template <class RealType, class Policy>
0399          RealType nc_beta_quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
0400          {
0401             static const char* function = "quantile(non_central_beta_distribution<%1%>, %1%)";
0402             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0403             typedef typename policies::normalise<
0404                Policy,
0405                policies::promote_float<false>,
0406                policies::promote_double<false>,
0407                policies::discrete_quantile<>,
0408                policies::assert_undefined<> >::type forwarding_policy;
0409
0410             value_type a = dist.alpha();
0411             value_type b = dist.beta();
0412             value_type l = dist.non_centrality();
0413             value_type r;
0414             if(!beta_detail::check_alpha(
0415                function,
0416                a, &r, Policy())
0417                ||
0418             !beta_detail::check_beta(
0419                function,
0420                b, &r, Policy())
0421                ||
0422             !detail::check_non_centrality(
0423                function,
0424                l,
0425                &r,
0426                Policy())
0427                ||
0428             !detail::check_probability(
0429                function,
0430                static_cast<value_type>(p),
0431                &r,
0432                Policy()))
0433                   return static_cast<RealType>(r);
0434             //
0435             // Special cases first:
0436             //
0437             if(p == 0)
0438                return comp
0439                ? 1.0f
0440                : 0.0f;
0441             if(p == 1)
0442                return !comp
0443                ? 1.0f
0444                : 0.0f;
0445
0446             value_type c = a + b + l / 2;
0447             value_type mean = 1 - (b / c) * (1 + l / (2 * c * c));
0448             /*
0449             //
0450             // Calculate a normal approximation to the quantile,
0451             // uses mean and variance approximations from:
0452             // Algorithm AS 310:
0453             // Computing the Non-Central Beta Distribution Function
0454             // R. Chattamvelli; R. Shanmugam
0455             // Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156.
0456             //
0457             // Unfortunately, when this is wrong it tends to be *very*
0458             // wrong, so it's disabled for now, even though it often
0459             // gets the initial guess quite close.  Probably we could
0460             // do much better by factoring in the skewness if only
0461             // we could calculate it....
0462             //
0463             value_type delta = l / 2;
0464             value_type delta2 = delta * delta;
0465             value_type delta3 = delta * delta2;
0466             value_type delta4 = delta2 * delta2;
0467             value_type G = c * (c + 1) + delta;
0468             value_type alpha = a + b;
0469             value_type alpha2 = alpha * alpha;
0470             value_type eta = (2 * alpha + 1) * (2 * alpha + 1) + 1;
0471             value_type H = 3 * alpha2 + 5 * alpha + 2;
0472             value_type F = alpha2 * (alpha + 1) + H * delta
0473                + (2 * alpha + 4) * delta2 + delta3;
0474             value_type P = (3 * alpha + 1) * (9 * alpha + 17)
0475                + 2 * alpha * (3 * alpha + 2) * (3 * alpha + 4) + 15;
0476             value_type Q = 54 * alpha2 + 162 * alpha + 130;
0477             value_type R = 6 * (6 * alpha + 11);
0478             value_type D = delta
0479                * (H * H + 2 * P * delta + Q * delta2 + R * delta3 + 9 * delta4);
0480             value_type variance = (b / G)
0481                * (1 + delta * (l * l + 3 * l + eta) / (G * G))
0482                - (b * b / F) * (1 + D / (F * F));
0483             value_type sd = sqrt(variance);
0484
0485             value_type guess = comp
0486                ? quantile(complement(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p))
0487                : quantile(normal_distribution<RealType, Policy>(static_cast<RealType>(mean), static_cast<RealType>(sd)), p);
0488
0489             if(guess >= 1)
0490                guess = mean;
0491             if(guess <= tools::min_value<value_type>())
0492                guess = mean;
0493             */
0494             value_type guess = mean;
0495             detail::nc_beta_quantile_functor<value_type, Policy>
0496                f(non_central_beta_distribution<value_type, Policy>(a, b, l), p, comp);
0497             tools::eps_tolerance<value_type> tol(policies::digits<RealType, Policy>());
0498             std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0499
0500             std::pair<value_type, value_type> ir
0501                = bracket_and_solve_root_01(
0502                   f, guess, value_type(2.5), true, tol,
0503                   max_iter, Policy());
0504             value_type result = ir.first + (ir.second - ir.first) / 2;
0505
0506             if(max_iter >= policies::get_max_root_iterations<Policy>())
0507             {
0508                return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
0509                   " either there is no answer to quantile of the non central beta distribution"
0510                   " or the answer is infinite.  Current best guess is %1%",
0511                   policies::checked_narrowing_cast<RealType, forwarding_policy>(
0512                      result,
0513                      function), Policy());
0514             }
0515             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0516                result,
0517                function);
0518          }
0519
0520          template <class T, class Policy>
0521          T non_central_beta_pdf(T a, T b, T lam, T x, T y, const Policy& pol)
0522          {
0523             BOOST_MATH_STD_USING
0524             //
0525             // Special cases:
0526             //
0527             if((x == 0) || (y == 0))
0528                return 0;
0529             //
0530             // Variables come first:
0531             //
0532             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0533             T errtol = boost::math::policies::get_epsilon<T, Policy>();
0534             T l2 = lam / 2;
0535             //
0536             // k is the starting point for iteration, and is the
0537             // maximum of the poisson weighting term:
0538             //
0539             long long k = lltrunc(l2);
0540             // Starting Poisson weight:
0541             T pois = gamma_p_derivative(T(k+1), l2, pol);
0542             // Starting beta term:
0543             T beta = x < y ?
0544                ibeta_derivative(a + k, b, x, pol)
0545                : ibeta_derivative(b, a + k, y, pol);
0546             T sum = 0;
0547             T poisf(pois);
0548             T betaf(beta);
0549
0550             //
0551             // Stable backwards recursion first:
0552             //
0553             std::uintmax_t count = k;
0554             for(auto i = k; i >= 0; --i)
0555             {
0556                T term = beta * pois;
0557                sum += term;
0558                if((fabs(term/sum) < errtol) || (term == 0))
0559                {
0560                   count = k - i;
0561                   break;
0562                }
0563                pois *= i / l2;
0564
0565                if (a + b + i != 1)
0566                {
0567                   beta *= (a + i - 1) / (x * (a + i + b - 1));
0568                }
0569             }
0570             for(auto i = k + 1; ; ++i)
0571             {
0572                poisf *= l2 / i;
0573                betaf *= x * (a + b + i - 1) / (a + i - 1);
0574
0575                T term = poisf * betaf;
0576                sum += term;
0577                if((fabs(term/sum) < errtol) || (term == 0))
0578                {
0579                   break;
0580                }
0581                if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
0582                {
0583                   return policies::raise_evaluation_error(
0584                      "pdf(non_central_beta_distribution<%1%>, %1%)",
0585                      "Series did not converge, closest value was %1%", sum, pol);
0586                }
0587             }
0588             return sum;
0589          }
0590
0591          template <class RealType, class Policy>
0592          RealType nc_beta_pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
0593          {
0594             BOOST_MATH_STD_USING
0595             static const char* function = "pdf(non_central_beta_distribution<%1%>, %1%)";
0596             typedef typename policies::evaluation<RealType, Policy>::type value_type;
0597             typedef typename policies::normalise<
0598                Policy,
0599                policies::promote_float<false>,
0600                policies::promote_double<false>,
0601                policies::discrete_quantile<>,
0602                policies::assert_undefined<> >::type forwarding_policy;
0603
0604             value_type a = dist.alpha();
0605             value_type b = dist.beta();
0606             value_type l = dist.non_centrality();
0607             value_type r;
0608             if(!beta_detail::check_alpha(
0609                function,
0610                a, &r, Policy())
0611                ||
0612             !beta_detail::check_beta(
0613                function,
0614                b, &r, Policy())
0615                ||
0616             !detail::check_non_centrality(
0617                function,
0618                l,
0619                &r,
0620                Policy())
0621                ||
0622             !beta_detail::check_x(
0623                function,
0624                static_cast<value_type>(x),
0625                &r,
0626                Policy()))
0627                   return static_cast<RealType>(r);
0628
0629             if(l == 0)
0630                return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
0631             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0632                non_central_beta_pdf(a, b, l, static_cast<value_type>(x), value_type(1 - static_cast<value_type>(x)), forwarding_policy()),
0633                "function");
0634          }
0635
0636          template <class T>
0637          struct hypergeometric_2F2_sum
0638          {
0639             typedef T result_type;
0640             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) {}
0641             T operator()()
0642             {
0643                T result = term;
0644                term *= a1 * a2 / (b1 * b2);
0645                a1 += 1;
0646                a2 += 1;
0647                b1 += 1;
0648                b2 += 1;
0649                k += 1;
0650                term /= k;
0651                term *= z;
0652                return result;
0653             }
0654             T a1, a2, b1, b2, z, term, k;
0655          };
0656
0657          template <class T, class Policy>
0658          T hypergeometric_2F2(T a1, T a2, T b1, T b2, T z, const Policy& pol)
0659          {
0660             typedef typename policies::evaluation<T, Policy>::type value_type;
0661
0662             const char* function = "boost::math::detail::hypergeometric_2F2<%1%>(%1%,%1%,%1%,%1%,%1%)";
0663
0664             hypergeometric_2F2_sum<value_type> s(a1, a2, b1, b2, z);
0665             std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0666
0667             value_type result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<value_type, Policy>(), max_iter);
0668
0669             policies::check_series_iterations<T>(function, max_iter, pol);
0670             return policies::checked_narrowing_cast<T, Policy>(result, function);
0671          }
0672
0673       } // namespace detail
0674
0675       template <class RealType = double, class Policy = policies::policy<> >
0676       class non_central_beta_distribution
0677       {
0678       public:
0679          typedef RealType value_type;
0680          typedef Policy policy_type;
0681
0682          non_central_beta_distribution(RealType a_, RealType b_, RealType lambda) : a(a_), b(b_), ncp(lambda)
0683          {
0684             const char* function = "boost::math::non_central_beta_distribution<%1%>::non_central_beta_distribution(%1%,%1%)";
0685             RealType r;
0686             beta_detail::check_alpha(
0687                function,
0688                a, &r, Policy());
0689             beta_detail::check_beta(
0690                function,
0691                b, &r, Policy());
0692             detail::check_non_centrality(
0693                function,
0694                lambda,
0695                &r,
0696                Policy());
0697          } // non_central_beta_distribution constructor.
0698
0699          RealType alpha() const
0700          { // Private data getter function.
0701             return a;
0702          }
0703          RealType beta() const
0704          { // Private data getter function.
0705             return b;
0706          }
0707          RealType non_centrality() const
0708          { // Private data getter function.
0709             return ncp;
0710          }
0711       private:
0712          // Data member, initialized by constructor.
0713          RealType a;   // alpha.
0714          RealType b;   // beta.
0715          RealType ncp; // non-centrality parameter
0716       }; // template <class RealType, class Policy> class non_central_beta_distribution
0717
0718       typedef non_central_beta_distribution<double> non_central_beta; // Reserved name of type double.
0719
0720       #ifdef __cpp_deduction_guides
0721       template <class RealType>
0722       non_central_beta_distribution(RealType,RealType,RealType)->non_central_beta_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0723       #endif
0724
0725       // Non-member functions to give properties of the distribution.
0726
0727       template <class RealType, class Policy>
0728       inline const std::pair<RealType, RealType> range(const non_central_beta_distribution<RealType, Policy>& /* dist */)
0729       { // Range of permissible values for random variable k.
0730          using boost::math::tools::max_value;
0731          return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0732       }
0733
0734       template <class RealType, class Policy>
0735       inline const std::pair<RealType, RealType> support(const non_central_beta_distribution<RealType, Policy>& /* dist */)
0736       { // Range of supported values for random variable k.
0737          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0738          using boost::math::tools::max_value;
0739          return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0740       }
0741
0742       template <class RealType, class Policy>
0743       inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist)
0744       { // mode.
0745          static const char* function = "mode(non_central_beta_distribution<%1%> const&)";
0746
0747          RealType a = dist.alpha();
0748          RealType b = dist.beta();
0749          RealType l = dist.non_centrality();
0750          RealType r;
0751          if(!beta_detail::check_alpha(
0752                function,
0753                a, &r, Policy())
0754                ||
0755             !beta_detail::check_beta(
0756                function,
0757                b, &r, Policy())
0758                ||
0759             !detail::check_non_centrality(
0760                function,
0761                l,
0762                &r,
0763                Policy()))
0764                   return static_cast<RealType>(r);
0765          RealType c = a + b + l / 2;
0766          RealType mean = 1 - (b / c) * (1 + l / (2 * c * c));
0767          return detail::generic_find_mode_01(
0768             dist,
0769             mean,
0770             function);
0771       }
0772
0773       //
0774       // We don't have the necessary information to implement
0775       // these at present.  These are just disabled for now,
0776       // prototypes retained so we can fill in the blanks
0777       // later:
0778       //
0779       template <class RealType, class Policy>
0780       inline RealType mean(const non_central_beta_distribution<RealType, Policy>& dist)
0781       {
0782          BOOST_MATH_STD_USING
0783          RealType a = dist.alpha();
0784          RealType b = dist.beta();
0785          RealType d = dist.non_centrality();
0786          RealType apb = a + b;
0787          return exp(-d / 2) * a * detail::hypergeometric_2F2<RealType, Policy>(1 + a, apb, a, 1 + apb, d / 2, Policy()) / apb;
0788       } // mean
0789
0790       template <class RealType, class Policy>
0791       inline RealType variance(const non_central_beta_distribution<RealType, Policy>& dist)
0792       {
0793          //
0794          // Relative error of this function may be arbitrarily large... absolute
0795          // error will be small however... that's the best we can do for now.
0796          //
0797          BOOST_MATH_STD_USING
0798          RealType a = dist.alpha();
0799          RealType b = dist.beta();
0800          RealType d = dist.non_centrality();
0801          RealType apb = a + b;
0802          RealType result = detail::hypergeometric_2F2(RealType(1 + a), apb, a, RealType(1 + apb), RealType(d / 2), Policy());
0803          result *= result * -exp(-d) * a * a / (apb * apb);
0804          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));
0805          return result;
0806       }
0807
0808       // RealType standard_deviation(const non_central_beta_distribution<RealType, Policy>& dist)
0809       // standard_deviation provided by derived accessors.
0810       template <class RealType, class Policy>
0811       inline RealType skewness(const non_central_beta_distribution<RealType, Policy>& /*dist*/)
0812       { // skewness = sqrt(l).
0813          const char* function = "boost::math::non_central_beta_distribution<%1%>::skewness()";
0814          typedef typename Policy::assert_undefined_type assert_type;
0815          static_assert(assert_type::value == 0, "Assert type is undefined.");
0816
0817          return policies::raise_evaluation_error<RealType>(
0818             function,
0819             "This function is not yet implemented, the only sensible result is %1%.",
0820             std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
0821       }
0822
0823       template <class RealType, class Policy>
0824       inline RealType kurtosis_excess(const non_central_beta_distribution<RealType, Policy>& /*dist*/)
0825       {
0826          const char* function = "boost::math::non_central_beta_distribution<%1%>::kurtosis_excess()";
0827          typedef typename Policy::assert_undefined_type assert_type;
0828          static_assert(assert_type::value == 0, "Assert type is undefined.");
0829
0830          return policies::raise_evaluation_error<RealType>(
0831             function,
0832             "This function is not yet implemented, the only sensible result is %1%.",
0833             std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
0834       } // kurtosis_excess
0835
0836       template <class RealType, class Policy>
0837       inline RealType kurtosis(const non_central_beta_distribution<RealType, Policy>& dist)
0838       {
0839          return kurtosis_excess(dist) + 3;
0840       }
0841
0842       template <class RealType, class Policy>
0843       inline RealType pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
0844       { // Probability Density/Mass Function.
0845          return detail::nc_beta_pdf(dist, x);
0846       } // pdf
0847
0848       template <class RealType, class Policy>
0849       RealType cdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
0850       {
0851          const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)";
0852             RealType a = dist.alpha();
0853             RealType b = dist.beta();
0854             RealType l = dist.non_centrality();
0855             RealType r;
0856             if(!beta_detail::check_alpha(
0857                function,
0858                a, &r, Policy())
0859                ||
0860             !beta_detail::check_beta(
0861                function,
0862                b, &r, Policy())
0863                ||
0864             !detail::check_non_centrality(
0865                function,
0866                l,
0867                &r,
0868                Policy())
0869                ||
0870             !beta_detail::check_x(
0871                function,
0872                x,
0873                &r,
0874                Policy()))
0875                   return static_cast<RealType>(r);
0876
0877          if(l == 0)
0878             return cdf(beta_distribution<RealType, Policy>(a, b), x);
0879
0880          return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, false, Policy());
0881       } // cdf
0882
0883       template <class RealType, class Policy>
0884       RealType cdf(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c)
0885       { // Complemented Cumulative Distribution Function
0886          const char* function = "boost::math::non_central_beta_distribution<%1%>::cdf(%1%)";
0887          non_central_beta_distribution<RealType, Policy> const& dist = c.dist;
0888             RealType a = dist.alpha();
0889             RealType b = dist.beta();
0890             RealType l = dist.non_centrality();
0891             RealType x = c.param;
0892             RealType r;
0893             if(!beta_detail::check_alpha(
0894                function,
0895                a, &r, Policy())
0896                ||
0897             !beta_detail::check_beta(
0898                function,
0899                b, &r, Policy())
0900                ||
0901             !detail::check_non_centrality(
0902                function,
0903                l,
0904                &r,
0905                Policy())
0906                ||
0907             !beta_detail::check_x(
0908                function,
0909                x,
0910                &r,
0911                Policy()))
0912                   return static_cast<RealType>(r);
0913
0914          if(l == 0)
0915             return cdf(complement(beta_distribution<RealType, Policy>(a, b), x));
0916
0917          return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, true, Policy());
0918       } // ccdf
0919
0920       template <class RealType, class Policy>
0921       inline RealType quantile(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& p)
0922       { // Quantile (or Percent Point) function.
0923          return detail::nc_beta_quantile(dist, p, false);
0924       } // quantile
0925
0926       template <class RealType, class Policy>
0927       inline RealType quantile(const complemented2_type<non_central_beta_distribution<RealType, Policy>, RealType>& c)
0928       { // Quantile (or Percent Point) function.
0929          return detail::nc_beta_quantile(c.dist, c.param, true);
0930       } // quantile complement.
0931
0932    } // namespace math
0933 } // namespace boost
0934
0935 // This include must be at the end, *after* the accessors
0936 // for this distribution have been defined, in order to
0937 // keep compilers that support two-phase lookup happy.
0938 #include <boost/math/distributions/detail/derived_accessors.hpp>
0939
0940 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_BETA_HPP
0941