206 lines
7.0 KiB
Plaintext
206 lines
7.0 KiB
Plaintext
[section:nc_beta_dist Noncentral Beta Distribution]
|
|
|
|
``#include <boost/math/distributions/non_central_beta.hpp>``
|
|
|
|
namespace boost{ namespace math{
|
|
|
|
template <class RealType = double,
|
|
class ``__Policy`` = ``__policy_class`` >
|
|
class non_central_beta_distribution;
|
|
|
|
typedef non_central_beta_distribution<> non_central_beta;
|
|
|
|
template <class RealType, class ``__Policy``>
|
|
class non_central_beta_distribution
|
|
{
|
|
public:
|
|
typedef RealType value_type;
|
|
typedef Policy policy_type;
|
|
|
|
// Constructor:
|
|
non_central_beta_distribution(RealType alpha, RealType beta, RealType lambda);
|
|
|
|
// Accessor to shape parameters:
|
|
RealType alpha()const;
|
|
RealType beta()const;
|
|
|
|
// Accessor to non-centrality parameter lambda:
|
|
RealType non_centrality()const;
|
|
};
|
|
|
|
}} // namespaces
|
|
|
|
The noncentral beta distribution is a generalization of the __beta_distrib.
|
|
|
|
It is defined as the ratio
|
|
[expression X = [chi][sub m][super 2]([lambda]) \/ ([chi][sub m][super 2]([lambda])
|
|
+ [chi][sub n][super 2])]
|
|
where [role serif_italic [chi][sub m][super 2]([lambda])]
|
|
is a noncentral [role serif_italic [chi][super 2]]
|
|
random variable with /m/ degrees of freedom, and [chi][sub n][super 2]
|
|
is a central [role serif_italic [chi][super 2] ] random variable with /n/ degrees of freedom.
|
|
|
|
This gives a PDF that can be expressed as a Poisson mixture
|
|
of beta distribution PDFs:
|
|
|
|
[equation nc_beta_ref1]
|
|
|
|
where P(i;[lambda]\/2) is the discrete Poisson probablity at /i/, with mean
|
|
[lambda]\/2, and I[sub x][super ']([alpha], [beta]) is the derivative of
|
|
the incomplete beta function. This leads to the usual form of the CDF
|
|
as:
|
|
|
|
[equation nc_beta_ref2]
|
|
|
|
The following graph illustrates how the distribution changes
|
|
for different values of [lambda]:
|
|
|
|
[graph nc_beta_pdf]
|
|
|
|
[h4 Member Functions]
|
|
|
|
non_central_beta_distribution(RealType a, RealType b, RealType lambda);
|
|
|
|
Constructs a noncentral beta distribution with shape parameters /a/ and /b/
|
|
and non-centrality parameter /lambda/.
|
|
|
|
Requires a > 0, b > 0 and lambda >= 0, otherwise calls __domain_error.
|
|
|
|
RealType alpha()const;
|
|
|
|
Returns the parameter /a/ from which this object was constructed.
|
|
|
|
RealType beta()const;
|
|
|
|
Returns the parameter /b/ from which this object was constructed.
|
|
|
|
RealType non_centrality()const;
|
|
|
|
Returns the parameter /lambda/ from which this object was constructed.
|
|
|
|
[h4 Non-member Accessors]
|
|
|
|
Most of the [link math_toolkit.dist_ref.nmp usual non-member accessor functions]
|
|
are supported: __cdf, __pdf, __quantile, __mean, __variance, __sd,
|
|
__median, __mode, __hazard, __chf, __range and __support.
|
|
|
|
Mean and variance are implemented using hypergeometric pfq functions and relations given in
|
|
[@http://reference.wolfram.com/mathematica/ref/NoncentralBetaDistribution.html Wolfram Noncentral Beta Distribution].
|
|
|
|
However, the following are not currently implemented:
|
|
__skewness, __kurtosis and __kurtosis_excess.
|
|
|
|
The domain of the random variable is \[0, 1\].
|
|
|
|
[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_beta_CDF]
|
|
|
|
[table_non_central_beta_CDF_complement]
|
|
|
|
Error rates for the PDF, the complement of the CDF and for the quantile
|
|
functions are broadly similar.
|
|
|
|
[h4 Tests]
|
|
|
|
There are two sets of test data used to verify this implementation:
|
|
firstly we can compare with a few sample values generated by the
|
|
[@http://www.r-project.org/ R library].
|
|
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, the crossover point
|
|
is taken to be the mean of the distribution: for this we use the
|
|
approximation due to: R. Chattamvelli and R. Shanmugam,
|
|
"Algorithm AS 310: Computing the Non-Central Beta Distribution Function",
|
|
Applied Statistics, Vol. 46, No. 1. (1997), pp. 146-156.
|
|
|
|
[equation nc_beta_ref3]
|
|
|
|
Then either the CDF or its complement is computed using the
|
|
relations:
|
|
|
|
[equation nc_beta_ref4]
|
|
|
|
The summation is performed by starting at i = [lambda]/2, and then recursing
|
|
in both directions, using the usual recurrence relations for the Poisson
|
|
PDF and incomplete beta functions. This is the "Method 2" described
|
|
by:
|
|
|
|
Denise Benton and K. Krishnamoorthy,
|
|
"Computing discrete mixtures of continuous
|
|
distributions: noncentral chisquare, noncentral t
|
|
and the distribution of the square of the sample
|
|
multiple correlation coefficient",
|
|
Computational Statistics & Data Analysis 43 (2003) 249-267.
|
|
|
|
Specific applications of the above formulae to the noncentral
|
|
beta distribution can be found in:
|
|
|
|
Russell V. Lenth,
|
|
"Algorithm AS 226: Computing Noncentral Beta Probabilities",
|
|
Applied Statistics, Vol. 36, No. 2. (1987), pp. 241-244.
|
|
|
|
H. Frick,
|
|
"Algorithm AS R84: A Remark on Algorithm AS 226: Computing Non-Central Beta
|
|
Probabilities", Applied Statistics, Vol. 39, No. 2. (1990), pp. 311-312.
|
|
|
|
Ming Long Lam,
|
|
"Remark AS R95: A Remark on Algorithm AS 226: Computing Non-Central Beta
|
|
Probabilities", Applied Statistics, Vol. 44, No. 4. (1995), pp. 551-552.
|
|
|
|
Harry O. Posten,
|
|
"An Effective Algorithm for the Noncentral Beta Distribution Function",
|
|
The American Statistician, Vol. 47, No. 2. (May, 1993), pp. 129-131.
|
|
|
|
R. Chattamvelli,
|
|
"A Note on the Noncentral Beta Distribution Function",
|
|
The American Statistician, Vol. 49, No. 2. (May, 1995), pp. 231-234.
|
|
|
|
Of these, the Posten reference provides the most complete overview,
|
|
and includes the modification starting iteration at [lambda]/2.
|
|
|
|
The main difference between this implementation and the above
|
|
references is the direct computation of the complement when most
|
|
efficient to do so, and the accumulation of the sum to -1 rather
|
|
than subtracting the result from 1 at the end: this can substantially
|
|
reduce the number of iterations required when the result is near 1.
|
|
|
|
The PDF is computed using the methodology of Benton and Krishnamoorthy
|
|
and the relation:
|
|
|
|
[equation nc_beta_ref1]
|
|
|
|
Quantiles are computed using a specially modified version of
|
|
__bracket_solve,
|
|
starting the search for the root at the mean of the distribution.
|
|
(A Cornish-Fisher type expansion was also tried, but while this gets
|
|
quite close to the root in many cases, when it is wrong it tends to
|
|
introduce quite pathological behaviour: more investigation in this
|
|
area is probably warranted).
|
|
|
|
[endsect] [/section:nc_beta_dist]
|
|
|
|
[/ nc_beta.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).
|
|
]
|
|
|