
The `owen` library for Haskell
20170509
SourceThe
owen
library is available on Github. The functions it provides are described and illustrated in this post.Owen \(T\)function
The Owen \(T\)function is defined by \[ T(h,a) = \frac{1}{2\pi}\int_0^a \frac{e^{\frac12 h^2(1+x^2)}}{1+x^2}\mathrm{d}x \] for \(a, h \in \mathbb{R}\). It corresponds to a certain probability and then its value always lies between \(0\) and \(1\).
Below we compare some obtained values of \(T(h, a)\) to the ones given by Wolfram up to \(20\) digits:
> import OwenT > owenT 0.1 0.1  0.01578338051718359918 3.469446951953614e18 > owenT 1 0.5  0.04306469112078536563 6.938893903907228e18 > owenT 0.5 0.1  0.01399302023628015424 1.734723475976807e18 > owenT 0.5 0.9  0.10007270175061385845 0.0 > owenT 5 0.5  0.00000014192549621069 2.6903883987164715e18 > owenT 2 0.99999  0.01111626714677311317 4.683753385137379e17
The
owenT
function allows infinite values of \(h\) and \(a\):> import OwenT > import Data.Number.Erf >  zero: > owenT (1/0) 3 0.0 > owenT (1/0) 3 0.0 >  equality: > owenT 1 (1/0) 7.932762696572854e2 > (1  normcdf 1)/2 7.932762696572854e2
The algorithm runs as follows. If \(0 \leq a \leq 1\), and \(0 \leq h \leq 8\), a series expansion is used, the one given by Owen (1956). The series is truncated after the \(50\)th term. If \(0\leq a \leq 1\), and \(h \geq 8\), an asymptotic approximation is used. Otherwise, the properties of the Owen \(T\)function are exploited to come down to the case \(0 \leq a \leq 1\). The main property involves the Gaussian cumulative function.
As shown by Owen (1956), the Owen \(T\)function can be used to evaluate the cumulative distribution function of the bivariate Gaussian distribution.
Noncentral Student distribution
The
pStudent
function of theowen
library evaluates the cumulative distribution function of the noncentral Student distribution with an integer number of degrees of freedom. For odd values, the algorithm resorts to the Owen \(T\)function. This algorithm is the one given by Owen (1965).Below we compare some values to the ones given by R:
> import Student > pStudent 0.5 2 2.5  0.022741814853305026 6.259229246019515e14 > pStudent 0.5 3 2.5  0.022741355255468675 1.8677073776451891e13
Owen Qfunction
The first Owen \(Q\)function is defined by \[ Q_1(\nu, t, \delta, R) = \frac{1}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac12(\nu2)}} \int_0^R \Phi\left(\frac{tx}{\sqrt{\nu}}\delta\right) x^{\nu1} e^{\frac{x^2}{2}} \mathrm{d}x \] for \(\nu >0\), \(t, \delta \in \mathbb{R}\) and \(R>0\).
Its value always lies between \(0\) and \(1\).
It is implemented in the
owen
library for integer values of \(\nu\), under the nameowenQ1
. The algorithm used is the one given by Owen (1965).The results are theoretically exact for even values of \(\nu\), up to the fact that the function resorts to the Gaussian cumulative function. For odd values of \(\nu\), it also resorts to the Owen \(T\)function.
Below we compare some obtained values of \(Q_1(\nu, t, \delta, R)\) to the ones given by Wolfram up to \(20\) digits:
> import OwenQ > owenQ1 1 3 2 5  0.52485658843054409291 2.220446049250313e16 > owenQ1 2 3 2 5  0.62938193306526904118 1.1102230246251565e16 > owenQ1 9 3 2 5  0.77030437685878389173 1.1102230246251565e16 > owenQ1 10 3 2 5  0.77383547873740988537 1.1102230246251565e16
The
owenQ1
function also returns a value for \(t,\delta=\pm \infty\):> import OwenQ > import Math.Gamma as G >  zero: > owenQ1 3 (1/0) 2 2 0.0 > owenQ1 3 2 (1/0) 2 0.0 >  equalities: > owenQ1 5 (1/0) 9 3 0.8909358420502276 > owenQ1 5 (1/0) 99 3 0.8909358420502276 > owenQ1 5 99 (1/0) 3 0.8909358420502276 > owenQ1 5 9 (1/0) 3 0.8909358420502276 > G.p (5/2) (3^2/2) 0.8909358420502276
References
Owen, D. B. (1956). Tables for computing bivariate normal probabilities. Ann. Math. Statist. 27, 107590.
Owen, D. B. (1965). A Special Case of a Bivariate NonCentral \(t\)Distribution. Biometrika 52, 437446.