mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-10-26 02:20:20 -04:00
274 lines
10 KiB
Plaintext
274 lines
10 KiB
Plaintext
[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).
|
|
]
|
|
|