1[section:nc_chi_squared_dist Noncentral Chi-Squared Distribution] 2 3``#include <boost/math/distributions/non_central_chi_squared.hpp>`` 4 5 namespace boost{ namespace math{ 6 7 template <class RealType = double, 8 class ``__Policy`` = ``__policy_class`` > 9 class non_central_chi_squared_distribution; 10 11 typedef non_central_chi_squared_distribution<> non_central_chi_squared; 12 13 template <class RealType, class ``__Policy``> 14 class non_central_chi_squared_distribution 15 { 16 public: 17 typedef RealType value_type; 18 typedef Policy policy_type; 19 20 // Constructor: 21 non_central_chi_squared_distribution(RealType v, RealType lambda); 22 23 // Accessor to degrees of freedom parameter v: 24 RealType degrees_of_freedom()const; 25 26 // Accessor to non centrality parameter lambda: 27 RealType non_centrality()const; 28 29 // Parameter finders: 30 static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p); 31 template <class A, class B, class C> 32 static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c); 33 34 static RealType find_non_centrality(RealType v, RealType x, RealType p); 35 template <class A, class B, class C> 36 static RealType find_non_centrality(const complemented3_type<A,B,C>& c); 37 }; 38 39 }} // namespaces 40 41The noncentral chi-squared distribution is a generalization of the 42__chi_squared_distrib. If ['X[sub i]] are /[nu]/ independent, normally 43distributed random variables with means /[mu][sub i]/ and variances 44['[sigma][sub i][super 2]], then the random variable 45 46[equation nc_chi_squ_ref1] 47 48is distributed according to the noncentral chi-squared distribution. 49 50The noncentral chi-squared distribution has two parameters: 51/[nu]/ which specifies the number of degrees of freedom 52(i.e. the number of ['X[sub i])], and [lambda] which is related to the 53mean of the random variables ['X[sub i]] by: 54 55[equation nc_chi_squ_ref2] 56 57(Note that some references define [lambda] as one half of the above sum). 58 59This leads to a PDF of: 60 61[equation nc_chi_squ_ref3] 62 63where ['f(x;k)] is the central chi-squared distribution PDF, and 64['I[sub v](x)] is a modified Bessel function of the first kind. 65 66The following graph illustrates how the distribution changes 67for different values of [lambda]: 68 69[graph nccs_pdf] 70 71[h4 Member Functions] 72 73 non_central_chi_squared_distribution(RealType v, RealType lambda); 74 75Constructs a Chi-Squared distribution with /v/ degrees of freedom 76and non-centrality parameter /lambda/. 77 78Requires v > 0 and lambda >= 0, otherwise calls __domain_error. 79 80 RealType degrees_of_freedom()const; 81 82Returns the parameter /v/ from which this object was constructed. 83 84 RealType non_centrality()const; 85 86Returns the parameter /lambda/ from which this object was constructed. 87 88 89 static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p); 90 91This function returns the number of degrees of freedom /v/ such that: 92`cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p` 93 94 template <class A, class B, class C> 95 static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c); 96 97When called with argument `boost::math::complement(lambda, x, q)` 98this function returns the number of degrees of freedom /v/ such that: 99 100`cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`. 101 102 static RealType find_non_centrality(RealType v, RealType x, RealType p); 103 104This function returns the non centrality parameter /lambda/ such that: 105 106`cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p` 107 108 template <class A, class B, class C> 109 static RealType find_non_centrality(const complemented3_type<A,B,C>& c); 110 111When called with argument `boost::math::complement(v, x, q)` 112this function returns the non centrality parameter /lambda/ such that: 113 114`cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`. 115 116[h4 Non-member Accessors] 117 118All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] 119that are generic to all distributions are supported: __usual_accessors. 120 121The domain of the random variable is \[0, +[infin]\]. 122 123[h4 Examples] 124 125There is a 126[link math_toolkit.stat_tut.weg.nccs_eg worked example] 127for the noncentral chi-squared distribution. 128 129[h4 Accuracy] 130 131The following table shows the peak errors 132(in units of [@http://en.wikipedia.org/wiki/Machine_epsilon epsilon]) 133found on various platforms with various floating point types. 134The failures in the comparison to the [@http://www.r-project.org/ R Math library], 135seem to be mostly in the corner cases when the probability would be very small. 136Unless otherwise specified any floating-point type that is narrower 137than the one shown will have __zero_error. 138 139[table_non_central_chi_squared_CDF] 140 141[table_non_central_chi_squared_CDF_complement] 142 143Error rates for the quantile 144functions are broadly similar. Special mention should go to 145the `mode` function: there is no closed form for this function, 146so it is evaluated numerically by finding the maxima of the PDF: 147in principal this can not produce an accuracy greater than the 148square root of the machine epsilon. 149 150[h4 Tests] 151 152There are two sets of test data used to verify this implementation: 153firstly we can compare with published data, for example with 154Table 6 of "Self-Validating Computations of Probabilities for 155Selected Central and Noncentral Univariate Probability Functions", 156Morgan C. Wang and William J. Kennedy, 157Journal of the American Statistical Association, 158Vol. 89, No. 427. (Sep., 1994), pp. 878-887. 159Secondly, we have tables of test data, computed with this 160implementation and using interval arithmetic - this data should 161be accurate to at least 50 decimal digits - and is the used for 162our accuracy tests. 163 164[h4 Implementation] 165 166The CDF and its complement are evaluated as follows: 167 168First we determine which of the two values (the CDF or its 169complement) is likely to be the smaller: for this we can use the 170relation due to Temme (see "Asymptotic and Numerical Aspects of the 171Noncentral Chi-Square Distribution", N. M. Temme, Computers Math. Applic. 172Vol 25, No. 5, 55-63, 1993) that: 173 174F([nu],[lambda];[nu]+[lambda]) [asymp] 0.5 175 176and so compute the CDF when the random variable is less than 177[nu]+[lambda], and its complement when the random variable is 178greater than [nu]+[lambda]. If necessary the computed result 179is then subtracted from 1 to give the desired result (the CDF or its 180complement). 181 182For small values of the non centrality parameter, the CDF is computed 183using the method of Ding (see "Algorithm AS 275: Computing the Non-Central 184#2 Distribution Function", Cherng G. Ding, Applied Statistics, Vol. 41, 185No. 2. (1992), pp. 478-482). This uses the following series representation: 186 187[equation nc_chi_squ_ref4] 188 189which requires just one call to __gamma_p_derivative with the subsequent 190terms being computed by recursion as shown above. 191 192For larger values of the non-centrality parameter, Ding's method can take 193an unreasonable number of terms before convergence is achieved. Furthermore, 194the largest term is not the first term, so in extreme cases the first term may 195be zero, leading to a zero result, even though the true value may be non-zero. 196 197Therefore, when the non-centrality parameter is greater than 200, the method due 198to Krishnamoorthy (see "Computing discrete mixtures of continuous distributions: 199noncentral chisquare, noncentral t and the distribution of the 200square of the sample multiple correlation coefficient", 201Denise Benton and K. Krishnamoorthy, Computational Statistics & 202Data Analysis, 43, (2003), 249-267) is used. 203 204This method uses the well known sum: 205 206[equation nc_chi_squ_ref5] 207 208Where ['P[sub a](x)] is the incomplete gamma function. 209 210The method starts at the [lambda]th term, which is where the Poisson weighting 211function achieves its maximum value, although this is not necessarily 212the largest overall term. Subsequent terms are calculated via the normal 213recurrence relations for the incomplete gamma function, and iteration proceeds 214both forwards and backwards until sufficient precision has been achieved. It 215should be noted that recurrence in the forwards direction of P[sub a](x) is 216numerically unstable. However, since we always start /after/ the largest 217term in the series, numeric instability is introduced more slowly than the 218series converges. 219 220Computation of the complement of the CDF uses an extension of Krishnamoorthy's 221method, given that: 222 223[equation nc_chi_squ_ref6] 224 225we can again start at the [lambda]'th term and proceed in both directions from 226there until the required precision is achieved. This time it is backwards 227recursion on the incomplete gamma function Q[sub a](x) which is unstable. 228However, as long as we start well /before/ the largest term, this is not an 229issue in practice. 230 231The PDF is computed directly using the relation: 232 233[equation nc_chi_squ_ref3] 234 235Where ['f(x; v)] is the PDF of the central __chi_squared_distrib and 236['I[sub v](x)] is a modified Bessel function, see __cyl_bessel_i. 237For small values of the 238non-centrality parameter the relation in terms of __cyl_bessel_i 239is used. However, this method fails for large values of the 240non-centrality parameter, so in that case the infinite sum is 241evaluated using the method of Benton and Krishnamoorthy, and 242the usual recurrence relations for successive terms. 243 244The quantile functions are computed by numeric inversion of the CDF. 245An improve starting guess is from 246Thomas Luu, 247[@http://discovery.ucl.ac.uk/1482128/, Fast and accurate parallel computation of quantile functions for random number generation, Doctoral Thesis, 2016]. 248 249There is no [@http://en.wikipedia.org/wiki/Closed_form closed form] 250for the mode of the noncentral chi-squared 251distribution: it is computed numerically by finding the maximum 252of the PDF. Likewise, the median is computed numerically via 253the quantile. 254 255The remaining non-member functions use the following formulas: 256 257[equation nc_chi_squ_ref7] 258 259Some analytic properties of noncentral distributions 260(particularly unimodality, and monotonicity of their modes) 261are surveyed and summarized by: 262 263Andrea van Aubel & Wolfgang Gawronski, Applied Mathematics and Computation, 141 (2003) 3-12. 264 265[endsect] [/section:nc_chi_squared_dist] 266 267[/ nc_chi_squared.qbk 268 Copyright 2008 John Maddock and Paul A. Bristow. 269 Distributed under the Boost Software License, Version 1.0. 270 (See accompanying file LICENSE_1_0.txt or copy at 271 http://www.boost.org/LICENSE_1_0.txt). 272] 273 274