math/distributions/non_central_beta.hpp

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