1. The owen library for Haskell

2017-05-09
Source

The 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.469446951953614e-18
> owenT 1 0.5     - 0.04306469112078536563
6.938893903907228e-18
> owenT 0.5 0.1   - 0.01399302023628015424
1.734723475976807e-18
> owenT 0.5 0.9   - 0.10007270175061385845
0.0
> owenT 5 0.5     - 0.00000014192549621069
-2.6903883987164715e-18
> owenT 2 0.99999 - 0.01111626714677311317
-4.683753385137379e-17

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.932762696572854e-2
> (1 - normcdf 1)/2
7.932762696572854e-2

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.

Non-central Student distribution

The pStudent function of the owen library evaluates the cumulative distribution function of the non-central 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.259229246019515e-14
> pStudent 0.5 3 2.5 - 0.022741355255468675
1.8677073776451891e-13

Owen Q-function

The first Owen $Q$-function is defined by $Q_1(\nu, t, \delta, R) = \frac{1}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac12(\nu-2)}} \int_0^R \Phi\left(\frac{tx}{\sqrt{\nu}}-\delta\right) x^{\nu-1} 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 name owenQ1. 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.220446049250313e-16
> owenQ1 2 3 2 5  - 0.62938193306526904118
-1.1102230246251565e-16
> owenQ1 9 3 2 5  - 0.77030437685878389173
-1.1102230246251565e-16
> owenQ1 10 3 2 5 - 0.77383547873740988537
1.1102230246251565e-16

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