<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Noncentral Beta Distribution</title>
<link rel="stylesheet" href="../../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.74.0">
<link rel="home" href="../../../../index.html" title="Math Toolkit">
<link rel="up" href="../dists.html" title="Distributions">
<link rel="prev" href="negative_binomial_dist.html" title="Negative Binomial Distribution">
<link rel="next" href="nc_chi_squared_dist.html" title="Noncentral Chi-Squared Distribution">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="negative_binomial_dist.html"><img src="../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../dists.html"><img src="../../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="nc_chi_squared_dist.html"><img src="../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h5 class="title">
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist"></a><a class="link" href="nc_beta_dist.html" title="Noncentral Beta Distribution"> Noncentral
Beta Distribution</a>
</h5></div></div></div>
<p>
</p>
<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">non_central_beta</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span></pre>
<p>
</p>
<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">RealType</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">,</span>
<span class="keyword">class</span> <a class="link" href="../../../policy.html" title="Policies">Policy</a> <span class="special">=</span> <a class="link" href="../../../policy/pol_ref/pol_ref_ref.html" title="Policy Class Reference">policies::policy<></a> <span class="special">></span>
<span class="keyword">class</span> <span class="identifier">non_central_beta_distribution</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">non_central_beta_distribution</span><span class="special"><></span> <span class="identifier">non_central_beta</span><span class="special">;</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">RealType</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../../policy.html" title="Policies">Policy</a><span class="special">></span>
<span class="keyword">class</span> <span class="identifier">non_central_beta_distribution</span>
<span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
<span class="keyword">typedef</span> <span class="identifier">RealType</span> <span class="identifier">value_type</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">Policy</span> <span class="identifier">policy_type</span><span class="special">;</span>
<span class="comment">// Constructor:
</span> <span class="identifier">non_central_beta_distribution</span><span class="special">(</span><span class="identifier">RealType</span> <span class="identifier">alpha</span><span class="special">,</span> <span class="identifier">RealType</span> <span class="identifier">beta</span><span class="special">,</span> <span class="identifier">RealType</span> <span class="identifier">lambda</span><span class="special">);</span>
<span class="comment">// Accessor to shape parameters:
</span> <span class="identifier">RealType</span> <span class="identifier">alpha</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
<span class="identifier">RealType</span> <span class="identifier">beta</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
<span class="comment">// Accessor to non-centrality parameter lambda:
</span> <span class="identifier">RealType</span> <span class="identifier">non_centrality</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
<span class="special">};</span>
<span class="special">}}</span> <span class="comment">// namespaces
</span></pre>
<p>
The noncentral beta distribution is a generalization of the <a class="link" href="beta_dist.html" title="Beta Distribution">Beta
Distribution</a>.
</p>
<p>
It is defined as the ratio X = χ<sub>m</sub><sup>2</sup>(λ) / (χ<sub>m</sub><sup>2</sup>(λ) + χ<sub>n</sub><sup>2</sup>) where χ<sub>m</sub><sup>2</sup>(λ) is a noncentral
χ<sup>2</sup>
random variable with <span class="emphasis"><em>m</em></span> degrees of freedom, and
χ<sub>n</sub><sup>2</sup>
is a central χ<sup>2</sup> random variable with <span class="emphasis"><em>n</em></span> degrees of
freedom.
</p>
<p>
This gives a PDF that can be expressed as a Poisson mixture of beta distribution
PDFs:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../equations/nc_beta_ref1.png"></span>
</p>
<p>
where P(i;λ/2) is the discrete Poisson probablity at <span class="emphasis"><em>i</em></span>,
with mean λ/2, and I<sub>x</sub><sup>'</sup>(α, β) is the derivative of the incomplete beta function.
This leads to the usual form of the CDF as:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../equations/nc_beta_ref2.png"></span>
</p>
<p>
The following graph illustrates how the distribution changes for different
values of λ:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../graphs/nc_beta_pdf.png" align="middle"></span>
</p>
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist.member_functions"></a><h5>
<a name="id1047231"></a>
<a class="link" href="nc_beta_dist.html#math_toolkit.dist.dist_ref.dists.nc_beta_dist.member_functions">Member
Functions</a>
</h5>
<pre class="programlisting"><span class="identifier">non_central_beta_distribution</span><span class="special">(</span><span class="identifier">RealType</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">RealType</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">RealType</span> <span class="identifier">lambda</span><span class="special">);</span>
</pre>
<p>
Constructs a noncentral beta distribution with shape parameters <span class="emphasis"><em>a</em></span>
and <span class="emphasis"><em>b</em></span> and non-centrality parameter <span class="emphasis"><em>lambda</em></span>.
</p>
<p>
Requires a > 0, b > 0 and lambda >= 0, otherwise calls <a class="link" href="../../../main_overview/error_handling.html#domain_error">domain_error</a>.
</p>
<pre class="programlisting"><span class="identifier">RealType</span> <span class="identifier">alpha</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
</pre>
<p>
Returns the parameter <span class="emphasis"><em>a</em></span> from which this object was
constructed.
</p>
<pre class="programlisting"><span class="identifier">RealType</span> <span class="identifier">beta</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
</pre>
<p>
Returns the parameter <span class="emphasis"><em>b</em></span> from which this object was
constructed.
</p>
<pre class="programlisting"><span class="identifier">RealType</span> <span class="identifier">non_centrality</span><span class="special">()</span><span class="keyword">const</span><span class="special">;</span>
</pre>
<p>
Returns the parameter <span class="emphasis"><em>lambda</em></span> from which this object
was constructed.
</p>
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist.non_member_accessors"></a><h5>
<a name="id1047417"></a>
<a class="link" href="nc_beta_dist.html#math_toolkit.dist.dist_ref.dists.nc_beta_dist.non_member_accessors">Non-member
Accessors</a>
</h5>
<p>
Most of the <a class="link" href="../nmp.html" title="Non-Member Properties">usual non-member
accessor functions</a> are supported: <a class="link" href="../nmp.html#math.dist.cdf">Cumulative
Distribution Function</a>, <a class="link" href="../nmp.html#math.dist.pdf">Probability
Density Function</a>, <a class="link" href="../nmp.html#math.dist.quantile">Quantile</a>,
<a class="link" href="../nmp.html#math.dist.median">median</a>, <a class="link" href="../nmp.html#math.dist.mode">mode</a>,
<a class="link" href="../nmp.html#math.dist.hazard">Hazard Function</a>, <a class="link" href="../nmp.html#math.dist.chf">Cumulative
Hazard Function</a>, <a class="link" href="../nmp.html#math.dist.range">range</a>
and <a class="link" href="../nmp.html#math.dist.support">support</a>.
</p>
<p>
However, the following are not currently implemented: <a class="link" href="../nmp.html#math.dist.mean">mean</a>,
<a class="link" href="../nmp.html#math.dist.variance">variance</a>, <a class="link" href="../nmp.html#math.dist.sd">standard
deviation</a>, <a class="link" href="../nmp.html#math.dist.skewness">skewness</a>,
<a class="link" href="../nmp.html#math.dist.kurtosis">kurtosis</a> and <a class="link" href="../nmp.html#math.dist.kurtosis_excess">kurtosis_excess</a>.
</p>
<p>
The domain of the random variable is [0, 1].
</p>
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist.accuracy"></a><h5>
<a name="id1047518"></a>
<a class="link" href="nc_beta_dist.html#math_toolkit.dist.dist_ref.dists.nc_beta_dist.accuracy">Accuracy</a>
</h5>
<p>
The following table shows the peak errors (in units of <a href="http://en.wikipedia.org/wiki/Machine_epsilon" target="_top">epsilon</a>)
found on various platforms with various floating point types. No comparison
to the <a href="http://www.r-project.org/" target="_top">R-2.5.1 Math library</a>,
or to the FORTRAN implementations of AS226 or AS310 are given since these
appear to only guarantee absolute error: this would causes our test harness
to assign an <span class="emphasis"><em>"infinite"</em></span> error to these
libraries for some of our test values when measuring <span class="emphasis"><em>relative
error</em></span>. Unless otherwise specified any floating-point type
that is narrower than the one shown will have <a class="link" href="../../../backgrounders/relative_error.html#zero_error">effectively
zero error</a>.
</p>
<div class="table">
<a name="id1047553"></a><p class="title"><b>Table 13. Errors In CDF of the Noncentral Beta</b></p>
<div class="table-contents"><table class="table" summary="Errors In CDF of the Noncentral Beta">
<colgroup>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
α, β,λ < 200
</p>
</th>
<th>
<p>
α,β,λ > 200
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32, Visual C++ 8
</p>
</td>
<td>
<p>
Peak=620 Mean=22
</p>
</td>
<td>
<p>
Peak=8670 Mean=1040
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
RedHat Linux IA32, gcc-4.1.1
</p>
</td>
<td>
<p>
Peak=825 Mean=50
</p>
</td>
<td>
<p>
Peak=2.5x10<sup>4</sup> Mean=4000
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Redhat Linux IA64, gcc-3.4.4
</p>
</td>
<td>
<p>
Peak=825 Mean=30
</p>
</td>
<td>
<p>
Peak=1.7x10<sup>4</sup> Mean=2500
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HPUX IA64, aCC A.06.06
</p>
</td>
<td>
<p>
Peak=420 Mean=50
</p>
</td>
<td>
<p>
Peak=9200 Mean=1200
</p>
</td>
</tr>
</tbody>
</table></div>
</div>
<br class="table-break"><p>
Error rates for the PDF, the complement of the CDF and for the quantile
functions are broadly similar.
</p>
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist.tests"></a><h5>
<a name="id1047759"></a>
<a class="link" href="nc_beta_dist.html#math_toolkit.dist.dist_ref.dists.nc_beta_dist.tests">Tests</a>
</h5>
<p>
There are two sets of test data used to verify this implementation: firstly
we can compare with a few sample values generated by the <a href="http://www.r-project.org/" target="_top">R
library</a>. 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.
</p>
<a name="math_toolkit.dist.dist_ref.dists.nc_beta_dist.implementation"></a><h5>
<a name="id1047783"></a>
<a class="link" href="nc_beta_dist.html#math_toolkit.dist.dist_ref.dists.nc_beta_dist.implementation">Implementation</a>
</h5>
<p>
The CDF and its complement are evaluated as follows:
</p>
<p>
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.
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../equations/nc_beta_ref3.png"></span>
</p>
<p>
Then either the CDF or its complement is computed using the relations:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../equations/nc_beta_ref4.png"></span>
</p>
<p>
The summation is performed by starting at i = λ/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:
</p>
<p>
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.
</p>
<p>
Specific applications of the above formulae to the noncentral beta distribution
can be found in:
</p>
<p>
Russell V. Lenth, "Algorithm AS 226: Computing Noncentral Beta Probabilities",
Applied Statistics, Vol. 36, No. 2. (1987), pp. 241-244.
</p>
<p>
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.
</p>
<p>
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.
</p>
<p>
Harry O. Posten, "An Effective Algorithm for the Noncentral Beta
Distribution Function", The American Statistician, Vol. 47, No.
2. (May, 1993), pp. 129-131.
</p>
<p>
R. Chattamvelli, "A Note on the Noncentral Beta Distribution Function",
The American Statistician, Vol. 49, No. 2. (May, 1995), pp. 231-234.
</p>
<p>
Of these, the Posten reference provides the most complete overview, and
includes the modification starting iteration at λ/2.
</p>
<p>
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.
</p>
<p>
The PDF is computed using the methodology of Benton and Krishnamoorthy
and the relation:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../../equations/nc_beta_ref1.png"></span>
</p>
<p>
Quantiles are computed using a specially modified version of <a class="link" href="../../../toolkit/internals1/roots2.html" title="Root Finding Without Derivatives">bracket_and_solve_root</a>,
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).
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006 , 2007, 2008, 2009 John Maddock, Paul A. Bristow,
Hubert Holin, Xiaogang Zhang, Bruno Lalande, Johan Råde, Gautam Sewani
and Thijs van den Berg<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="negative_binomial_dist.html"><img src="../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../dists.html"><img src="../../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="nc_chi_squared_dist.html"><img src="../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
