123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273 |
- [section:nc_chi_squared_dist Noncentral Chi-Squared Distribution]
- ``#include <boost/math/distributions/non_central_chi_squared.hpp>``
- namespace boost{ namespace math{
- template <class RealType = double,
- class ``__Policy`` = ``__policy_class`` >
- class non_central_chi_squared_distribution;
- typedef non_central_chi_squared_distribution<> non_central_chi_squared;
- template <class RealType, class ``__Policy``>
- class non_central_chi_squared_distribution
- {
- public:
- typedef RealType value_type;
- typedef Policy policy_type;
- // Constructor:
- non_central_chi_squared_distribution(RealType v, RealType lambda);
- // Accessor to degrees of freedom parameter v:
- RealType degrees_of_freedom()const;
- // Accessor to non centrality parameter lambda:
- RealType non_centrality()const;
- // Parameter finders:
- static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p);
- template <class A, class B, class C>
- static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c);
- static RealType find_non_centrality(RealType v, RealType x, RealType p);
- template <class A, class B, class C>
- static RealType find_non_centrality(const complemented3_type<A,B,C>& c);
- };
- }} // namespaces
- The noncentral chi-squared distribution is a generalization of the
- __chi_squared_distrib. If ['X[sub i]] are /[nu]/ independent, normally
- distributed random variables with means /[mu][sub i]/ and variances
- ['[sigma][sub i][super 2]], then the random variable
- [equation nc_chi_squ_ref1]
- is distributed according to the noncentral chi-squared distribution.
- The noncentral chi-squared distribution has two parameters:
- /[nu]/ which specifies the number of degrees of freedom
- (i.e. the number of ['X[sub i])], and [lambda] which is related to the
- mean of the random variables ['X[sub i]] by:
- [equation nc_chi_squ_ref2]
- (Note that some references define [lambda] as one half of the above sum).
- This leads to a PDF of:
- [equation nc_chi_squ_ref3]
- where ['f(x;k)] is the central chi-squared distribution PDF, and
- ['I[sub v](x)] is a modified Bessel function of the first kind.
- The following graph illustrates how the distribution changes
- for different values of [lambda]:
- [graph nccs_pdf]
- [h4 Member Functions]
- non_central_chi_squared_distribution(RealType v, RealType lambda);
- Constructs a Chi-Squared distribution with /v/ degrees of freedom
- and non-centrality parameter /lambda/.
- Requires v > 0 and lambda >= 0, otherwise calls __domain_error.
- RealType degrees_of_freedom()const;
- Returns the parameter /v/ from which this object was constructed.
- RealType non_centrality()const;
- Returns the parameter /lambda/ from which this object was constructed.
- static RealType find_degrees_of_freedom(RealType lambda, RealType x, RealType p);
- This function returns the number of degrees of freedom /v/ such that:
- `cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p`
- template <class A, class B, class C>
- static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c);
- When called with argument `boost::math::complement(lambda, x, q)`
- this function returns the number of degrees of freedom /v/ such that:
- `cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`.
- static RealType find_non_centrality(RealType v, RealType x, RealType p);
- This function returns the non centrality parameter /lambda/ such that:
- `cdf(non_central_chi_squared<RealType, Policy>(v, lambda), x) == p`
- template <class A, class B, class C>
- static RealType find_non_centrality(const complemented3_type<A,B,C>& c);
- When called with argument `boost::math::complement(v, x, q)`
- this function returns the non centrality parameter /lambda/ such that:
- `cdf(complement(non_central_chi_squared<RealType, Policy>(v, lambda), x)) == q`.
- [h4 Non-member Accessors]
- All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions]
- that are generic to all distributions are supported: __usual_accessors.
- The domain of the random variable is \[0, +[infin]\].
- [h4 Examples]
- There is a
- [link math_toolkit.stat_tut.weg.nccs_eg worked example]
- for the noncentral chi-squared distribution.
- [h4 Accuracy]
- The following table shows the peak errors
- (in units of [@http://en.wikipedia.org/wiki/Machine_epsilon epsilon])
- found on various platforms with various floating point types.
- The failures in the comparison to the [@http://www.r-project.org/ R Math library],
- seem to be mostly in the corner cases when the probablity would be very small.
- Unless otherwise specified any floating-point type that is narrower
- than the one shown will have __zero_error.
- [table_non_central_chi_squared_CDF]
- [table_non_central_chi_squared_CDF_complement]
- Error rates for the quantile
- functions are broadly similar. Special mention should go to
- the `mode` function: there is no closed form for this function,
- so it is evaluated numerically by finding the maxima of the PDF:
- in principal this can not produce an accuracy greater than the
- square root of the machine epsilon.
- [h4 Tests]
- There are two sets of test data used to verify this implementation:
- firstly we can compare with published data, for example with
- Table 6 of "Self-Validating Computations of Probabilities for
- Selected Central and Noncentral Univariate Probability Functions",
- Morgan C. Wang and William J. Kennedy,
- Journal of the American Statistical Association,
- Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
- Secondly, we have tables of test data, computed with this
- implementation and using interval arithmetic - this data should
- be accurate to at least 50 decimal digits - and is the used for
- our accuracy tests.
- [h4 Implementation]
- The CDF and its complement are evaluated as follows:
- First we determine which of the two values (the CDF or its
- complement) is likely to be the smaller: for this we can use the
- relation due to Temme (see "Asymptotic and Numerical Aspects of the
- Noncentral Chi-Square Distribution", N. M. Temme, Computers Math. Applic.
- Vol 25, No. 5, 55-63, 1993) that:
- F([nu],[lambda];[nu]+[lambda]) [asymp] 0.5
- and so compute the CDF when the random variable is less than
- [nu]+[lambda], and its complement when the random variable is
- greater than [nu]+[lambda]. If necessary the computed result
- is then subtracted from 1 to give the desired result (the CDF or its
- complement).
- For small values of the non centrality parameter, the CDF is computed
- using the method of Ding (see "Algorithm AS 275: Computing the Non-Central
- #2 Distribution Function", Cherng G. Ding, Applied Statistics, Vol. 41,
- No. 2. (1992), pp. 478-482). This uses the following series representation:
- [equation nc_chi_squ_ref4]
- which requires just one call to __gamma_p_derivative with the subsequent
- terms being computed by recursion as shown above.
- For larger values of the non-centrality parameter, Ding's method can take
- an unreasonable number of terms before convergence is achieved. Furthermore,
- the largest term is not the first term, so in extreme cases the first term may
- be zero, leading to a zero result, even though the true value may be non-zero.
- Therefore, when the non-centrality parameter is greater than 200, the method due
- to Krishnamoorthy (see "Computing discrete mixtures of continuous distributions:
- noncentral chisquare, noncentral t and the distribution of the
- square of the sample multiple correlation coefficient",
- Denise Benton and K. Krishnamoorthy, Computational Statistics &
- Data Analysis, 43, (2003), 249-267) is used.
- This method uses the well known sum:
- [equation nc_chi_squ_ref5]
- Where ['P[sub a](x)] is the incomplete gamma function.
- The method starts at the [lambda]th term, which is where the Poisson weighting
- function achieves its maximum value, although this is not necessarily
- the largest overall term. Subsequent terms are calculated via the normal
- recurrence relations for the incomplete gamma function, and iteration proceeds
- both forwards and backwards until sufficient precision has been achieved. It
- should be noted that recurrence in the forwards direction of P[sub a](x) is
- numerically unstable. However, since we always start /after/ the largest
- term in the series, numeric instability is introduced more slowly than the
- series converges.
- Computation of the complement of the CDF uses an extension of Krishnamoorthy's
- method, given that:
- [equation nc_chi_squ_ref6]
- we can again start at the [lambda]'th term and proceed in both directions from
- there until the required precision is achieved. This time it is backwards
- recursion on the incomplete gamma function Q[sub a](x) which is unstable.
- However, as long as we start well /before/ the largest term, this is not an
- issue in practice.
- The PDF is computed directly using the relation:
- [equation nc_chi_squ_ref3]
- Where ['f(x; v)] is the PDF of the central __chi_squared_distrib and
- ['I[sub v](x)] is a modified Bessel function, see __cyl_bessel_i.
- For small values of the
- non-centrality parameter the relation in terms of __cyl_bessel_i
- is used. However, this method fails for large values of the
- non-centrality parameter, so in that case the infinite sum is
- evaluated using the method of Benton and Krishnamoorthy, and
- the usual recurrence relations for successive terms.
- The quantile functions are computed by numeric inversion of the CDF.
- An improve starting quess is from
- Thomas Luu,
- [@http://discovery.ucl.ac.uk/1482128/, Fast and accurate parallel computation of quantile functions for random number generation, Doctorial Thesis, 2016].
- There is no [@http://en.wikipedia.org/wiki/Closed_form closed form]
- for the mode of the noncentral chi-squared
- distribution: it is computed numerically by finding the maximum
- of the PDF. Likewise, the median is computed numerically via
- the quantile.
- The remaining non-member functions use the following formulas:
- [equation nc_chi_squ_ref7]
- Some analytic properties of noncentral distributions
- (particularly unimodality, and monotonicity of their modes)
- are surveyed and summarized by:
- Andrea van Aubel & Wolfgang Gawronski, Applied Mathematics and Computation, 141 (2003) 3-12.
- [endsect] [/section:nc_chi_squared_dist]
- [/ nc_chi_squared.qbk
- Copyright 2008 John Maddock and Paul A. Bristow.
- Distributed under the Boost Software License, Version 1.0.
- (See accompanying file LICENSE_1_0.txt or copy at
- http://www.boost.org/LICENSE_1_0.txt).
- ]
|