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