[ceph.git] / ceph / src / boost / libs / math / doc / sf / owens_t.qbk

[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.

[: ['T(h, 0) = 0]]

[:['T(0, a) = [frac12][pi] arctan(a)]]

[:['T(h, 1) = [frac12] G(h) \[1 - G(h)\]]]

[:['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).
]
Commit	Line	Data
7c673cae FG	1	[section:owens_t Owen's T function]
	2
	3	[h4 Synopsis]
	4
	5	``
	6	#include <boost/math/special_functions/owens_t.hpp>
	7	``
	8
	9	namespace boost{ namespace math{
	10
	11	template <class T>
	12	``__sf_result`` owens_t(T h, T a);
	13
	14	template <class T, class ``__Policy``>
	15	``__sf_result`` owens_t(T h, T a, const ``__Policy``&);
	16
	17	}} // namespaces
	18
	19	[h4 Description]
	20
	21	Returns the
	22	[@http://en.wikipedia.org/wiki/Owen%27s_T_function Owens_t function]
	23	of ['h] and ['a].
	24
	25	[optional_policy]
	26
	27	[sixemspace][sixemspace][equation owens_t]
	28
	29	[$../graphs/plot_owens_t.png]
	30
	31	The function `owens_t(h, a)` gives the probability
	32	of the event ['(X > h and 0 < Y < a * X)],
	33	where ['X] and ['Y] are independent standard normal random variables.
	34
	35	For h and a > 0, T(h,a),
	36	gives the volume of an uncorrelated bivariate normal distribution
	37	with zero means and unit variances over the area between
	38	['y = ax] and ['y = 0] and to the right of ['x = h].
	39
	40	That is the area shaded in the figure below (Owens 1956).
	41
	42	[graph owens_integration_area]
	43
	44	and is also illustrated by a 3D plot.
	45
	46	[$../graphs/plot_owens_3d_xyp.png]
	47
	48	This function is used in the computation of the __skew_normal_distrib.
	49	It is also used in the computation of bivariate and
	50	multivariate normal distribution probabilities.
	51	The return type of this function is computed using the __arg_promotion_rules:
	52	the result is of type `double` when T is an integer type, and type T otherwise.
	53
	54	Owen's original paper (page 1077) provides some additional corner cases.
	55
	56	[: ['T(h, 0) = 0]]
	57
	58	[:['T(0, a) = [frac12][pi] arctan(a)]]
	59
	60	[:['T(h, 1) = [frac12] G(h) \[1 - G(h)\]]]
	61
	62	[:['T(h, [infin]) = G(\|h\|)]]
	63
	64	where G(h) is the univariate normal with zero mean and unit variance integral from -[infin] to h.
65
66	[h4 Accuracy]
67
68	Over the built-in types and range tested,
69	errors are less than 10 * std::numeric_limits<RealType>::epsilon().
70
71	[table_owens_t]
72
73	[h4 Testing]
74
75	Test data was generated by Patefield and Tandy algorithms T1 and T4,
76	and also the suggested reference routine T7.
77
78	* T1 was rejected if the result was too small compared to `atan(a)` (ie cancellation),
79	* T4 was rejected if there was no convergence,
80	* Both were rejected if they didn't agree.
81
82	Over the built-in types and range tested,
83	errors are less than 10 std::numeric_limits<RealType>::epsilon().
84
85	However, that there was a whole domain (large ['h], small ['a])
86	where it was not possible to generate any reliable test values
87	(all the methods got rejected for one reason or another).
88
89	There are also two sets of sanity tests: spot values are computed using __Mathematica and __R.
90
91
92	[h4 Implementation]
93
94	The function was proposed and evaluated by
95	[@http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.aoms/1177728074
96	Donald. B. Owen, Tables for computing bivariate normal probabilities,
97	Ann. Math. Statist., 27, 1075-1090 (1956)].
98
99	The algorithms of Patefield, M. and Tandy, D.
100	"Fast and accurate Calculation of Owen's T-Function", Journal of Statistical Software, 5 (5), 1 - 25 (2000)
101	are adapted for C++ with arbitrary RealType.
102
103	The Patefield-Tandy algorithm provides six methods of evalualution (T1 to T6);
104	the best method is selected according to the values of ['a] and ['h].
105	See the original paper and the source in
106	[@../../../../boost/math/special_functions/owens_t.hpp owens_t.hpp] for details.
107
108	The Patefield-Tandy algorithm is accurate to approximately 20 decimal places, so for
109	types with greater precision we use:
110
111	* A modified version of T1 which folds the calculation of ['atan(h)] into the T1 series
112	(to avoid subtracting two values similar in magnitude), and then accelerates the
113	resulting alternating series using method 1 from H. Cohen, F. Rodriguez Villegas, D. Zagier,
114	"Convergence acceleration of alternating series", Bonn, (1991). The result is valid everywhere,
115	but doesn't always converge, or may become too divergent in the first few terms to sum accurately.
116	This is used for ['ah < 1].
117	* A modified version of T2 which is accelerated in the same manner as T1. This is used for ['h > 1].
118	* A version of T4 only when both T1 and T2 have failed to produce an accurate answer.
119	* Fallback to the Patefiled Tandy algorithm when all the above methods fail: this happens not at all
120	for our test data at 100 decimal digits precision. However, there is a difficult area when
121	['a] is very close to 1 and the precision increases which may cause this to happen in very exceptional
122	circumstances.
123
124	Using the above algorithm and a 100-decimal digit type, results accurate to 80 decimal places were obtained
125	in the difficult area where ['a] is close to 1, and greater than 95 decimal places elsewhere.
126
127	[endsect] [/section:owens_t The owens_t Function]
128
129	[/
130	Copyright 2012 Bejamin Sobotta, John Maddock and Paul A. Bristow.
131	Distributed under the Boost Software License, Version 1.0.
132	(See accompanying file LICENSE_1_0.txt or copy at
133	http://www.boost.org/LICENSE_1_0.txt).
134	]
135
136