136 lines
4.8 KiB
Plaintext
136 lines
4.8 KiB
Plaintext
[section:owens_t Owen's T function]
|
|
|
|
[h4 Synopsis]
|
|
|
|
``
|
|
#include <boost/math/special_functions/owens_t.hpp>
|
|
``
|
|
|
|
namespace boost{ namespace math{
|
|
|
|
template <class T>
|
|
``__sf_result`` owens_t(T h, T a);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` owens_t(T h, T a, const ``__Policy``&);
|
|
|
|
}} // namespaces
|
|
|
|
[h4 Description]
|
|
|
|
Returns the
|
|
[@http://en.wikipedia.org/wiki/Owen%27s_T_function Owens_t function]
|
|
of ['h] and ['a].
|
|
|
|
[optional_policy]
|
|
|
|
[sixemspace][sixemspace][equation owens_t]
|
|
|
|
[$../graphs/plot_owens_t.png]
|
|
|
|
The function `owens_t(h, a)` gives the probability
|
|
of the event ['(X > h and 0 < Y < a * X)],
|
|
where ['X] and ['Y] are independent standard normal random variables.
|
|
|
|
For h and a > 0, T(h,a),
|
|
gives the volume of an uncorrelated bivariate normal distribution
|
|
with zero means and unit variances over the area between
|
|
['y = ax] and ['y = 0] and to the right of ['x = h].
|
|
|
|
That is the area shaded in the figure below (Owens 1956).
|
|
|
|
[graph owens_integration_area]
|
|
|
|
and is also illustrated by a 3D plot.
|
|
|
|
[$../graphs/plot_owens_3d_xyp.png]
|
|
|
|
This function is used in the computation of the __skew_normal_distrib.
|
|
It is also used in the computation of bivariate and
|
|
multivariate normal distribution probabilities.
|
|
The return type of this function is computed using the __arg_promotion_rules:
|
|
the result is of type `double` when T is an integer type, and type T otherwise.
|
|
|
|
Owen's original paper (page 1077) provides some additional corner cases.
|
|
|
|
[expression ['T(h, 0) = 0]]
|
|
|
|
[expression ['T(0, a) = [frac12][pi] arctan(a)]]
|
|
|
|
[expression ['T(h, 1) = [frac12] G(h) \[1 - G(h)\]]]
|
|
|
|
[expression ['T(h, [infin]) = G(|h|)]]
|
|
|
|
where G(h) is the univariate normal with zero mean and unit variance integral from -[infin] to h.
|
|
|
|
[h4 Accuracy]
|
|
|
|
Over the built-in types and range tested,
|
|
errors are less than 10 * std::numeric_limits<RealType>::epsilon().
|
|
|
|
[table_owens_t]
|
|
|
|
[h4 Testing]
|
|
|
|
Test data was generated by Patefield and Tandy algorithms T1 and T4,
|
|
and also the suggested reference routine T7.
|
|
|
|
* T1 was rejected if the result was too small compared to `atan(a)` (ie cancellation),
|
|
* T4 was rejected if there was no convergence,
|
|
* Both were rejected if they didn't agree.
|
|
|
|
Over the built-in types and range tested,
|
|
errors are less than 10 std::numeric_limits<RealType>::epsilon().
|
|
|
|
However, that there was a whole domain (large ['h], small ['a])
|
|
where it was not possible to generate any reliable test values
|
|
(all the methods got rejected for one reason or another).
|
|
|
|
There are also two sets of sanity tests: spot values are computed using __Mathematica and __R.
|
|
|
|
[h4 Implementation]
|
|
|
|
The function was proposed and evaluated by
|
|
[@http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.aoms/1177728074
|
|
Donald. B. Owen, Tables for computing bivariate normal probabilities,
|
|
Ann. Math. Statist., 27, 1075-1090 (1956)].
|
|
|
|
The algorithms of Patefield, M. and Tandy, D.
|
|
"Fast and accurate Calculation of Owen's T-Function", Journal of Statistical Software, 5 (5), 1 - 25 (2000)
|
|
are adapted for C++ with arbitrary RealType.
|
|
|
|
The Patefield-Tandy algorithm provides six methods of evalualution (T1 to T6);
|
|
the best method is selected according to the values of ['a] and ['h].
|
|
See the original paper and the source in
|
|
[@../../../../boost/math/special_functions/owens_t.hpp owens_t.hpp] for details.
|
|
|
|
The Patefield-Tandy algorithm is accurate to approximately 20 decimal places, so for
|
|
types with greater precision we use:
|
|
|
|
* A modified version of T1 which folds the calculation of ['atan(h)] into the T1 series
|
|
(to avoid subtracting two values similar in magnitude), and then accelerates the
|
|
resulting alternating series using method 1 from H. Cohen, F. Rodriguez Villegas, D. Zagier,
|
|
"Convergence acceleration of alternating series", Bonn, (1991). The result is valid everywhere,
|
|
but doesn't always converge, or may become too divergent in the first few terms to sum accurately.
|
|
This is used for ['ah < 1].
|
|
* A modified version of T2 which is accelerated in the same manner as T1. This is used for ['h > 1].
|
|
* A version of T4 only when both T1 and T2 have failed to produce an accurate answer.
|
|
* Fallback to the Patefiled Tandy algorithm when all the above methods fail: this happens not at all
|
|
for our test data at 100 decimal digits precision. However, there is a difficult area when
|
|
['a] is very close to 1 and the precision increases which may cause this to happen in very exceptional
|
|
circumstances.
|
|
|
|
Using the above algorithm and a 100-decimal digit type, results accurate to 80 decimal places were obtained
|
|
in the difficult area where ['a] is close to 1, and greater than 95 decimal places elsewhere.
|
|
|
|
[endsect] [/section:owens_t The owens_t Function]
|
|
|
|
[/
|
|
Copyright 2012 Bejamin Sobotta, 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).
|
|
]
|
|
|
|
|