non_central_chi_squared.hpp 40 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993
  1. // boost\math\distributions\non_central_chi_squared.hpp
  2. // Copyright John Maddock 2008.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0.
  5. // (See accompanying file LICENSE_1_0.txt
  6. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  7. #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  8. #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  9. #include <boost/math/distributions/fwd.hpp>
  10. #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
  11. #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
  12. #include <boost/math/special_functions/round.hpp> // for llround
  13. #include <boost/math/distributions/complement.hpp> // complements
  14. #include <boost/math/distributions/chi_squared.hpp> // central distribution
  15. #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
  16. #include <boost/math/special_functions/fpclassify.hpp> // isnan.
  17. #include <boost/math/tools/roots.hpp> // for root finding.
  18. #include <boost/math/distributions/detail/generic_mode.hpp>
  19. #include <boost/math/distributions/detail/generic_quantile.hpp>
  20. namespace boost
  21. {
  22. namespace math
  23. {
  24. template <class RealType, class Policy>
  25. class non_central_chi_squared_distribution;
  26. namespace detail{
  27. template <class T, class Policy>
  28. T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
  29. {
  30. //
  31. // Computes the complement of the Non-Central Chi-Square
  32. // Distribution CDF by summing a weighted sum of complements
  33. // of the central-distributions. The weighting factor is
  34. // a Poisson Distribution.
  35. //
  36. // This is an application of the technique described in:
  37. //
  38. // Computing discrete mixtures of continuous
  39. // distributions: noncentral chisquare, noncentral t
  40. // and the distribution of the square of the sample
  41. // multiple correlation coefficient.
  42. // D. Benton, K. Krishnamoorthy.
  43. // Computational Statistics & Data Analysis 43 (2003) 249 - 267
  44. //
  45. BOOST_MATH_STD_USING
  46. // Special case:
  47. if(x == 0)
  48. return 1;
  49. //
  50. // Initialize the variables we'll be using:
  51. //
  52. T lambda = theta / 2;
  53. T del = f / 2;
  54. T y = x / 2;
  55. std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  56. T errtol = boost::math::policies::get_epsilon<T, Policy>();
  57. T sum = init_sum;
  58. //
  59. // k is the starting location for iteration, we'll
  60. // move both forwards and backwards from this point.
  61. // k is chosen as the peek of the Poisson weights, which
  62. // will occur *before* the largest term.
  63. //
  64. long long k = llround(lambda, pol);
  65. // Forwards and backwards Poisson weights:
  66. T poisf = boost::math::gamma_p_derivative(static_cast<T>(1 + k), lambda, pol);
  67. T poisb = poisf * k / lambda;
  68. // Initial forwards central chi squared term:
  69. T gamf = boost::math::gamma_q(del + k, y, pol);
  70. // Forwards and backwards recursion terms on the central chi squared:
  71. T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
  72. T xtermb = xtermf * (del + k) / y;
  73. // Initial backwards central chi squared term:
  74. T gamb = gamf - xtermb;
  75. //
  76. // Forwards iteration first, this is the
  77. // stable direction for the gamma function
  78. // recurrences:
  79. //
  80. long long i;
  81. for(i = k; static_cast<std::uintmax_t>(i-k) < max_iter; ++i)
  82. {
  83. T term = poisf * gamf;
  84. sum += term;
  85. poisf *= lambda / (i + 1);
  86. gamf += xtermf;
  87. xtermf *= y / (del + i + 1);
  88. if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
  89. break;
  90. }
  91. //Error check:
  92. if(static_cast<std::uintmax_t>(i-k) >= max_iter)
  93. return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
  94. //
  95. // Now backwards iteration: the gamma
  96. // function recurrences are unstable in this
  97. // direction, we rely on the terms diminishing in size
  98. // faster than we introduce cancellation errors.
  99. // For this reason it's very important that we start
  100. // *before* the largest term so that backwards iteration
  101. // is strictly converging.
  102. //
  103. for(i = k - 1; i >= 0; --i)
  104. {
  105. T term = poisb * gamb;
  106. sum += term;
  107. poisb *= i / lambda;
  108. xtermb *= (del + i) / y;
  109. gamb -= xtermb;
  110. if((sum == 0) || (fabs(term / sum) < errtol))
  111. break;
  112. }
  113. return sum;
  114. }
  115. template <class T, class Policy>
  116. T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
  117. {
  118. //
  119. // This is an implementation of:
  120. //
  121. // Algorithm AS 275:
  122. // Computing the Non-Central #2 Distribution Function
  123. // Cherng G. Ding
  124. // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
  125. //
  126. // This uses a stable forward iteration to sum the
  127. // CDF, unfortunately this can not be used for large
  128. // values of the non-centrality parameter because:
  129. // * The first term may underflow to zero.
  130. // * We may need an extra-ordinary number of terms
  131. // before we reach the first *significant* term.
  132. //
  133. BOOST_MATH_STD_USING
  134. // Special case:
  135. if(x == 0)
  136. return 0;
  137. T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
  138. T lambda = theta / 2;
  139. T vk = exp(-lambda);
  140. T uk = vk;
  141. T sum = init_sum + tk * vk;
  142. if(sum == 0)
  143. return sum;
  144. std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  145. T errtol = boost::math::policies::get_epsilon<T, Policy>();
  146. int i;
  147. T lterm(0), term(0);
  148. for(i = 1; static_cast<std::uintmax_t>(i) < max_iter; ++i)
  149. {
  150. tk = tk * x / (f + 2 * i);
  151. uk = uk * lambda / i;
  152. vk = vk + uk;
  153. lterm = term;
  154. term = vk * tk;
  155. sum += term;
  156. if((fabs(term / sum) < errtol) && (term <= lterm))
  157. break;
  158. }
  159. //Error check:
  160. if(static_cast<std::uintmax_t>(i) >= max_iter)
  161. return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
  162. return sum;
  163. }
  164. template <class T, class Policy>
  165. T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
  166. {
  167. //
  168. // This is taken more or less directly from:
  169. //
  170. // Computing discrete mixtures of continuous
  171. // distributions: noncentral chisquare, noncentral t
  172. // and the distribution of the square of the sample
  173. // multiple correlation coefficient.
  174. // D. Benton, K. Krishnamoorthy.
  175. // Computational Statistics & Data Analysis 43 (2003) 249 - 267
  176. //
  177. // We're summing a Poisson weighting term multiplied by
  178. // a central chi squared distribution.
  179. //
  180. BOOST_MATH_STD_USING
  181. // Special case:
  182. if(y == 0)
  183. return 0;
  184. std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  185. T errtol = boost::math::policies::get_epsilon<T, Policy>();
  186. T errorf(0), errorb(0);
  187. T x = y / 2;
  188. T del = lambda / 2;
  189. //
  190. // Starting location for the iteration, we'll iterate
  191. // both forwards and backwards from this point. The
  192. // location chosen is the maximum of the Poisson weight
  193. // function, which ocurrs *after* the largest term in the
  194. // sum.
  195. //
  196. long long k = llround(del, pol);
  197. T a = n / 2 + k;
  198. // Central chi squared term for forward iteration:
  199. T gamkf = boost::math::gamma_p(a, x, pol);
  200. if(lambda == 0)
  201. return gamkf;
  202. // Central chi squared term for backward iteration:
  203. T gamkb = gamkf;
  204. // Forwards Poisson weight:
  205. T poiskf = gamma_p_derivative(static_cast<T>(k+1), del, pol);
  206. // Backwards Poisson weight:
  207. T poiskb = poiskf;
  208. // Forwards gamma function recursion term:
  209. T xtermf = boost::math::gamma_p_derivative(a, x, pol);
  210. // Backwards gamma function recursion term:
  211. T xtermb = xtermf * x / a;
  212. T sum = init_sum + poiskf * gamkf;
  213. if(sum == 0)
  214. return sum;
  215. int i = 1;
  216. //
  217. // Backwards recursion first, this is the stable
  218. // direction for gamma function recurrences:
  219. //
  220. while(i <= k)
  221. {
  222. xtermb *= (a - i + 1) / x;
  223. gamkb += xtermb;
  224. poiskb = poiskb * (k - i + 1) / del;
  225. errorf = errorb;
  226. errorb = gamkb * poiskb;
  227. sum += errorb;
  228. if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
  229. break;
  230. ++i;
  231. }
  232. i = 1;
  233. //
  234. // Now forwards recursion, the gamma function
  235. // recurrence relation is unstable in this direction,
  236. // so we rely on the magnitude of successive terms
  237. // decreasing faster than we introduce cancellation error.
  238. // For this reason it's vital that k is chosen to be *after*
  239. // the largest term, so that successive forward iterations
  240. // are strictly (and rapidly) converging.
  241. //
  242. do
  243. {
  244. xtermf = xtermf * x / (a + i - 1);
  245. gamkf = gamkf - xtermf;
  246. poiskf = poiskf * del / (k + i);
  247. errorf = poiskf * gamkf;
  248. sum += errorf;
  249. ++i;
  250. }while((fabs(errorf / sum) > errtol) && (static_cast<std::uintmax_t>(i) < max_iter));
  251. //Error check:
  252. if(static_cast<std::uintmax_t>(i) >= max_iter)
  253. return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
  254. return sum;
  255. }
  256. template <class T, class Policy>
  257. T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
  258. {
  259. //
  260. // As above but for the PDF:
  261. //
  262. BOOST_MATH_STD_USING
  263. std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  264. T errtol = boost::math::policies::get_epsilon<T, Policy>();
  265. T x2 = x / 2;
  266. T n2 = n / 2;
  267. T l2 = lambda / 2;
  268. T sum = 0;
  269. long long k = lltrunc(l2);
  270. T pois = gamma_p_derivative(static_cast<T>(k + 1), l2, pol) * gamma_p_derivative(static_cast<T>(n2 + k), x2);
  271. if(pois == 0)
  272. return 0;
  273. T poisb = pois;
  274. for(long long i = k; ; ++i)
  275. {
  276. sum += pois;
  277. if(pois / sum < errtol)
  278. break;
  279. if(static_cast<std::uintmax_t>(i - k) >= max_iter)
  280. return policies::raise_evaluation_error("pdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
  281. pois *= l2 * x2 / ((i + 1) * (n2 + i));
  282. }
  283. for(long long i = k - 1; i >= 0; --i)
  284. {
  285. poisb *= (i + 1) * (n2 + i) / (l2 * x2);
  286. sum += poisb;
  287. if(poisb / sum < errtol)
  288. break;
  289. }
  290. return sum / 2;
  291. }
  292. template <class RealType, class Policy>
  293. inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
  294. {
  295. typedef typename policies::evaluation<RealType, Policy>::type value_type;
  296. typedef typename policies::normalise<
  297. Policy,
  298. policies::promote_float<false>,
  299. policies::promote_double<false>,
  300. policies::discrete_quantile<>,
  301. policies::assert_undefined<> >::type forwarding_policy;
  302. BOOST_MATH_STD_USING
  303. value_type result;
  304. if(l == 0)
  305. return invert == false ? cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x) : cdf(complement(boost::math::chi_squared_distribution<RealType, Policy>(k), x));
  306. else if(x > k + l)
  307. {
  308. // Complement is the smaller of the two:
  309. result = detail::non_central_chi_square_q(
  310. static_cast<value_type>(x),
  311. static_cast<value_type>(k),
  312. static_cast<value_type>(l),
  313. forwarding_policy(),
  314. static_cast<value_type>(invert ? 0 : -1));
  315. invert = !invert;
  316. }
  317. else if(l < 200)
  318. {
  319. // For small values of the non-centrality parameter
  320. // we can use Ding's method:
  321. result = detail::non_central_chi_square_p_ding(
  322. static_cast<value_type>(x),
  323. static_cast<value_type>(k),
  324. static_cast<value_type>(l),
  325. forwarding_policy(),
  326. static_cast<value_type>(invert ? -1 : 0));
  327. }
  328. else
  329. {
  330. // For largers values of the non-centrality
  331. // parameter Ding's method will consume an
  332. // extra-ordinary number of terms, and worse
  333. // may return zero when the result is in fact
  334. // finite, use Krishnamoorthy's method instead:
  335. result = detail::non_central_chi_square_p(
  336. static_cast<value_type>(x),
  337. static_cast<value_type>(k),
  338. static_cast<value_type>(l),
  339. forwarding_policy(),
  340. static_cast<value_type>(invert ? -1 : 0));
  341. }
  342. if(invert)
  343. result = -result;
  344. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  345. result,
  346. "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
  347. }
  348. template <class T, class Policy>
  349. struct nccs_quantile_functor
  350. {
  351. nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
  352. : dist(d), target(t), comp(c) {}
  353. T operator()(const T& x)
  354. {
  355. return comp ?
  356. target - cdf(complement(dist, x))
  357. : cdf(dist, x) - target;
  358. }
  359. private:
  360. non_central_chi_squared_distribution<T,Policy> dist;
  361. T target;
  362. bool comp;
  363. };
  364. template <class RealType, class Policy>
  365. RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
  366. {
  367. BOOST_MATH_STD_USING
  368. static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
  369. typedef typename policies::evaluation<RealType, Policy>::type value_type;
  370. typedef typename policies::normalise<
  371. Policy,
  372. policies::promote_float<false>,
  373. policies::promote_double<false>,
  374. policies::discrete_quantile<>,
  375. policies::assert_undefined<> >::type forwarding_policy;
  376. value_type k = dist.degrees_of_freedom();
  377. value_type l = dist.non_centrality();
  378. value_type r;
  379. if(!detail::check_df(
  380. function,
  381. k, &r, Policy())
  382. ||
  383. !detail::check_non_centrality(
  384. function,
  385. l,
  386. &r,
  387. Policy())
  388. ||
  389. !detail::check_probability(
  390. function,
  391. static_cast<value_type>(p),
  392. &r,
  393. Policy()))
  394. return static_cast<RealType>(r);
  395. //
  396. // Special cases get short-circuited first:
  397. //
  398. if(p == 0)
  399. return comp ? policies::raise_overflow_error<RealType>(function, 0, Policy()) : 0;
  400. if(p == 1)
  401. return comp ? 0 : policies::raise_overflow_error<RealType>(function, 0, Policy());
  402. //
  403. // This is Pearson's approximation to the quantile, see
  404. // Pearson, E. S. (1959) "Note on an approximation to the distribution of
  405. // noncentral chi squared", Biometrika 46: 364.
  406. // See also:
  407. // "A comparison of approximations to percentiles of the noncentral chi2-distribution",
  408. // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1-2) : 57-76.
  409. // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile.
  410. //
  411. value_type b = -(l * l) / (k + 3 * l);
  412. value_type c = (k + 3 * l) / (k + 2 * l);
  413. value_type ff = (k + 2 * l) / (c * c);
  414. value_type guess;
  415. if(comp)
  416. {
  417. guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
  418. }
  419. else
  420. {
  421. guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
  422. }
  423. //
  424. // Sometimes guess goes very small or negative, in that case we have
  425. // to do something else for the initial guess, this approximation
  426. // was provided in a private communication from Thomas Luu, PhD candidate,
  427. // University College London. It's an asymptotic expansion for the
  428. // quantile which usually gets us within an order of magnitude of the
  429. // correct answer.
  430. // Fast and accurate parallel computation of quantile functions for random number generation,
  431. // Thomas LuuDoctorial Thesis 2016
  432. // http://discovery.ucl.ac.uk/1482128/
  433. //
  434. if(guess < 0.005)
  435. {
  436. value_type pp = comp ? 1 - p : p;
  437. //guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k);
  438. guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k * boost::math::tgamma(k / 2, forwarding_policy()), (2 / k));
  439. if(guess == 0)
  440. guess = tools::min_value<value_type>();
  441. }
  442. value_type result = detail::generic_quantile(
  443. non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
  444. p,
  445. guess,
  446. comp,
  447. function);
  448. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  449. result,
  450. function);
  451. }
  452. template <class RealType, class Policy>
  453. RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
  454. {
  455. BOOST_MATH_STD_USING
  456. static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
  457. typedef typename policies::evaluation<RealType, Policy>::type value_type;
  458. typedef typename policies::normalise<
  459. Policy,
  460. policies::promote_float<false>,
  461. policies::promote_double<false>,
  462. policies::discrete_quantile<>,
  463. policies::assert_undefined<> >::type forwarding_policy;
  464. value_type k = dist.degrees_of_freedom();
  465. value_type l = dist.non_centrality();
  466. value_type r;
  467. if(!detail::check_df(
  468. function,
  469. k, &r, Policy())
  470. ||
  471. !detail::check_non_centrality(
  472. function,
  473. l,
  474. &r,
  475. Policy())
  476. ||
  477. !detail::check_positive_x(
  478. function,
  479. (value_type)x,
  480. &r,
  481. Policy()))
  482. return static_cast<RealType>(r);
  483. if(l == 0)
  484. return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
  485. // Special case:
  486. if(x == 0)
  487. return 0;
  488. if(l > 50)
  489. {
  490. r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
  491. }
  492. else
  493. {
  494. r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
  495. if(fabs(r) >= tools::log_max_value<RealType>() / 4)
  496. {
  497. r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
  498. }
  499. else
  500. {
  501. r = exp(r);
  502. r = 0.5f * r
  503. * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
  504. }
  505. }
  506. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  507. r,
  508. function);
  509. }
  510. template <class RealType, class Policy>
  511. struct degrees_of_freedom_finder
  512. {
  513. degrees_of_freedom_finder(
  514. RealType lam_, RealType x_, RealType p_, bool c)
  515. : lam(lam_), x(x_), p(p_), comp(c) {}
  516. RealType operator()(const RealType& v)
  517. {
  518. non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
  519. return comp ?
  520. RealType(p - cdf(complement(d, x)))
  521. : RealType(cdf(d, x) - p);
  522. }
  523. private:
  524. RealType lam;
  525. RealType x;
  526. RealType p;
  527. bool comp;
  528. };
  529. template <class RealType, class Policy>
  530. inline RealType find_degrees_of_freedom(
  531. RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
  532. {
  533. const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
  534. if((p == 0) || (q == 0))
  535. {
  536. //
  537. // Can't a thing if one of p and q is zero:
  538. //
  539. return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE
  540. RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
  541. }
  542. degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
  543. tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
  544. std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
  545. //
  546. // Pick an initial guess that we know will give us a probability
  547. // right around 0.5.
  548. //
  549. RealType guess = x - lam;
  550. if(guess < 1)
  551. guess = 1;
  552. std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
  553. f, guess, RealType(2), false, tol, max_iter, pol);
  554. RealType result = ir.first + (ir.second - ir.first) / 2;
  555. if(max_iter >= policies::get_max_root_iterations<Policy>())
  556. {
  557. return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
  558. " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
  559. }
  560. return result;
  561. }
  562. template <class RealType, class Policy>
  563. struct non_centrality_finder
  564. {
  565. non_centrality_finder(
  566. RealType v_, RealType x_, RealType p_, bool c)
  567. : v(v_), x(x_), p(p_), comp(c) {}
  568. RealType operator()(const RealType& lam)
  569. {
  570. non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
  571. return comp ?
  572. RealType(p - cdf(complement(d, x)))
  573. : RealType(cdf(d, x) - p);
  574. }
  575. private:
  576. RealType v;
  577. RealType x;
  578. RealType p;
  579. bool comp;
  580. };
  581. template <class RealType, class Policy>
  582. inline RealType find_non_centrality(
  583. RealType v, RealType x, RealType p, RealType q, const Policy& pol)
  584. {
  585. const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
  586. if((p == 0) || (q == 0))
  587. {
  588. //
  589. // Can't do a thing if one of p and q is zero:
  590. //
  591. return policies::raise_evaluation_error<RealType>(function, "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE
  592. RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
  593. }
  594. non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
  595. tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
  596. std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
  597. //
  598. // Pick an initial guess that we know will give us a probability
  599. // right around 0.5.
  600. //
  601. RealType guess = x - v;
  602. if(guess < 1)
  603. guess = 1;
  604. std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
  605. f, guess, RealType(2), false, tol, max_iter, pol);
  606. RealType result = ir.first + (ir.second - ir.first) / 2;
  607. if(max_iter >= policies::get_max_root_iterations<Policy>())
  608. {
  609. return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
  610. " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
  611. }
  612. return result;
  613. }
  614. }
  615. template <class RealType = double, class Policy = policies::policy<> >
  616. class non_central_chi_squared_distribution
  617. {
  618. public:
  619. typedef RealType value_type;
  620. typedef Policy policy_type;
  621. non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
  622. {
  623. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
  624. RealType r;
  625. detail::check_df(
  626. function,
  627. df, &r, Policy());
  628. detail::check_non_centrality(
  629. function,
  630. ncp,
  631. &r,
  632. Policy());
  633. } // non_central_chi_squared_distribution constructor.
  634. RealType degrees_of_freedom() const
  635. { // Private data getter function.
  636. return df;
  637. }
  638. RealType non_centrality() const
  639. { // Private data getter function.
  640. return ncp;
  641. }
  642. static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
  643. {
  644. const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
  645. typedef typename policies::evaluation<RealType, Policy>::type eval_type;
  646. typedef typename policies::normalise<
  647. Policy,
  648. policies::promote_float<false>,
  649. policies::promote_double<false>,
  650. policies::discrete_quantile<>,
  651. policies::assert_undefined<> >::type forwarding_policy;
  652. eval_type result = detail::find_degrees_of_freedom(
  653. static_cast<eval_type>(lam),
  654. static_cast<eval_type>(x),
  655. static_cast<eval_type>(p),
  656. static_cast<eval_type>(1-p),
  657. forwarding_policy());
  658. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  659. result,
  660. function);
  661. }
  662. template <class A, class B, class C>
  663. static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
  664. {
  665. const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
  666. typedef typename policies::evaluation<RealType, Policy>::type eval_type;
  667. typedef typename policies::normalise<
  668. Policy,
  669. policies::promote_float<false>,
  670. policies::promote_double<false>,
  671. policies::discrete_quantile<>,
  672. policies::assert_undefined<> >::type forwarding_policy;
  673. eval_type result = detail::find_degrees_of_freedom(
  674. static_cast<eval_type>(c.dist),
  675. static_cast<eval_type>(c.param1),
  676. static_cast<eval_type>(1-c.param2),
  677. static_cast<eval_type>(c.param2),
  678. forwarding_policy());
  679. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  680. result,
  681. function);
  682. }
  683. static RealType find_non_centrality(RealType v, RealType x, RealType p)
  684. {
  685. const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
  686. typedef typename policies::evaluation<RealType, Policy>::type eval_type;
  687. typedef typename policies::normalise<
  688. Policy,
  689. policies::promote_float<false>,
  690. policies::promote_double<false>,
  691. policies::discrete_quantile<>,
  692. policies::assert_undefined<> >::type forwarding_policy;
  693. eval_type result = detail::find_non_centrality(
  694. static_cast<eval_type>(v),
  695. static_cast<eval_type>(x),
  696. static_cast<eval_type>(p),
  697. static_cast<eval_type>(1-p),
  698. forwarding_policy());
  699. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  700. result,
  701. function);
  702. }
  703. template <class A, class B, class C>
  704. static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
  705. {
  706. const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
  707. typedef typename policies::evaluation<RealType, Policy>::type eval_type;
  708. typedef typename policies::normalise<
  709. Policy,
  710. policies::promote_float<false>,
  711. policies::promote_double<false>,
  712. policies::discrete_quantile<>,
  713. policies::assert_undefined<> >::type forwarding_policy;
  714. eval_type result = detail::find_non_centrality(
  715. static_cast<eval_type>(c.dist),
  716. static_cast<eval_type>(c.param1),
  717. static_cast<eval_type>(1-c.param2),
  718. static_cast<eval_type>(c.param2),
  719. forwarding_policy());
  720. return policies::checked_narrowing_cast<RealType, forwarding_policy>(
  721. result,
  722. function);
  723. }
  724. private:
  725. // Data member, initialized by constructor.
  726. RealType df; // degrees of freedom.
  727. RealType ncp; // non-centrality parameter
  728. }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
  729. typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
  730. #ifdef __cpp_deduction_guides
  731. template <class RealType>
  732. non_central_chi_squared_distribution(RealType,RealType)->non_central_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  733. #endif
  734. // Non-member functions to give properties of the distribution.
  735. template <class RealType, class Policy>
  736. inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
  737. { // Range of permissible values for random variable k.
  738. using boost::math::tools::max_value;
  739. return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer?
  740. }
  741. template <class RealType, class Policy>
  742. inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
  743. { // Range of supported values for random variable k.
  744. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  745. using boost::math::tools::max_value;
  746. return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  747. }
  748. template <class RealType, class Policy>
  749. inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  750. { // Mean of poisson distribution = lambda.
  751. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
  752. RealType k = dist.degrees_of_freedom();
  753. RealType l = dist.non_centrality();
  754. RealType r;
  755. if(!detail::check_df(
  756. function,
  757. k, &r, Policy())
  758. ||
  759. !detail::check_non_centrality(
  760. function,
  761. l,
  762. &r,
  763. Policy()))
  764. return static_cast<RealType>(r);
  765. return k + l;
  766. } // mean
  767. template <class RealType, class Policy>
  768. inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  769. { // mode.
  770. static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
  771. RealType k = dist.degrees_of_freedom();
  772. RealType l = dist.non_centrality();
  773. RealType r;
  774. if(!detail::check_df(
  775. function,
  776. k, &r, Policy())
  777. ||
  778. !detail::check_non_centrality(
  779. function,
  780. l,
  781. &r,
  782. Policy()))
  783. return static_cast<RealType>(r);
  784. bool asymptotic_mode = k < l/4;
  785. RealType starting_point = asymptotic_mode ? k + l - RealType(3) : RealType(1) + k;
  786. return detail::generic_find_mode(dist, starting_point, function);
  787. }
  788. template <class RealType, class Policy>
  789. inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  790. { // variance.
  791. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
  792. RealType k = dist.degrees_of_freedom();
  793. RealType l = dist.non_centrality();
  794. RealType r;
  795. if(!detail::check_df(
  796. function,
  797. k, &r, Policy())
  798. ||
  799. !detail::check_non_centrality(
  800. function,
  801. l,
  802. &r,
  803. Policy()))
  804. return static_cast<RealType>(r);
  805. return 2 * (2 * l + k);
  806. }
  807. // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  808. // standard_deviation provided by derived accessors.
  809. template <class RealType, class Policy>
  810. inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  811. { // skewness = sqrt(l).
  812. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
  813. RealType k = dist.degrees_of_freedom();
  814. RealType l = dist.non_centrality();
  815. RealType r;
  816. if(!detail::check_df(
  817. function,
  818. k, &r, Policy())
  819. ||
  820. !detail::check_non_centrality(
  821. function,
  822. l,
  823. &r,
  824. Policy()))
  825. return static_cast<RealType>(r);
  826. BOOST_MATH_STD_USING
  827. return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
  828. }
  829. template <class RealType, class Policy>
  830. inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  831. {
  832. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
  833. RealType k = dist.degrees_of_freedom();
  834. RealType l = dist.non_centrality();
  835. RealType r;
  836. if(!detail::check_df(
  837. function,
  838. k, &r, Policy())
  839. ||
  840. !detail::check_non_centrality(
  841. function,
  842. l,
  843. &r,
  844. Policy()))
  845. return static_cast<RealType>(r);
  846. return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
  847. } // kurtosis_excess
  848. template <class RealType, class Policy>
  849. inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
  850. {
  851. return kurtosis_excess(dist) + 3;
  852. }
  853. template <class RealType, class Policy>
  854. inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
  855. { // Probability Density/Mass Function.
  856. return detail::nccs_pdf(dist, x);
  857. } // pdf
  858. template <class RealType, class Policy>
  859. RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
  860. {
  861. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
  862. RealType k = dist.degrees_of_freedom();
  863. RealType l = dist.non_centrality();
  864. RealType r;
  865. if(!detail::check_df(
  866. function,
  867. k, &r, Policy())
  868. ||
  869. !detail::check_non_centrality(
  870. function,
  871. l,
  872. &r,
  873. Policy())
  874. ||
  875. !detail::check_positive_x(
  876. function,
  877. x,
  878. &r,
  879. Policy()))
  880. return static_cast<RealType>(r);
  881. return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
  882. } // cdf
  883. template <class RealType, class Policy>
  884. RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
  885. { // Complemented Cumulative Distribution Function
  886. const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
  887. non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
  888. RealType x = c.param;
  889. RealType k = dist.degrees_of_freedom();
  890. RealType l = dist.non_centrality();
  891. RealType r;
  892. if(!detail::check_df(
  893. function,
  894. k, &r, Policy())
  895. ||
  896. !detail::check_non_centrality(
  897. function,
  898. l,
  899. &r,
  900. Policy())
  901. ||
  902. !detail::check_positive_x(
  903. function,
  904. x,
  905. &r,
  906. Policy()))
  907. return static_cast<RealType>(r);
  908. return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
  909. } // ccdf
  910. template <class RealType, class Policy>
  911. inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
  912. { // Quantile (or Percent Point) function.
  913. return detail::nccs_quantile(dist, p, false);
  914. } // quantile
  915. template <class RealType, class Policy>
  916. inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
  917. { // Quantile (or Percent Point) function.
  918. return detail::nccs_quantile(c.dist, c.param, true);
  919. } // quantile complement.
  920. } // namespace math
  921. } // namespace boost
  922. // This include must be at the end, *after* the accessors
  923. // for this distribution have been defined, in order to
  924. // keep compilers that support two-phase lookup happy.
  925. #include <boost/math/distributions/detail/derived_accessors.hpp>
  926. #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP