Owen cumulative functions

The Owen distribution

Let $Z \sim {\cal N}(0,1)$ and X ∼ χν2 be two independent random variables. For real numbers δ1 and δ2, define the two random variables $$ T_1 = \frac{Z+\delta_1}{\sqrt{X/\nu}} \quad \text{and} \quad\; T_2 = \frac{Z+\delta_2}{\sqrt{X/\nu}}. $$

Both T1 and T2 follow a non-central Student distribution. The number of degrees of freedom is ν for each of them, and their respective non-centrality parameters are δ1 and δ2 respectively.

Owen (1965) studied the distribution of the pair (T1, T2).

The four Owen cumulative functions are $$ \begin{align} O_1(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \leq t_1, T_2 \leq t_2), \\ O_2(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \leq t_1, T_2 \geq t_2), \\ O_3(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \geq t_1, T_2 \geq t_2), \\ O_4(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \geq t_1, T_2 \leq t_2). \end{align} $$

Owen provided an efficient way to evaluate these functions when ν is an integer number. Owen’s algorithms are implemented in the OwenQ package.

For δ1 > δ2, these four functions are implemented in the OwenQ package under the respective names powen1, powen2, powen3 and powen4. For general values of δ1 and δ2, they are implemented under the respective names psbt1, psbt2, psbt3 and psbt4.

Non-central Student distribution

Owen (1965) also provided an algorithm to evaluate the cumulative distribution function of a univariate non-central Student distribution with an integer number of degrees of freedom. This evaluation is performed by the function ptOwen of the OwenQ package.

ptOwen(q=1, nu=3, delta=2)
## [1] 0.1573494
pt(q=1, df=3, ncp=2)
## [1] 0.1573494

It is known that the pt function is not reliable when the non-centrality parameter ncp is large. Below we compare the values given by ptOwen and pt to the value given by Wolfram|Alpha (returned by the command N[CDF[NoncentralStudentTDistribution[4,70],80]]):

p1 <- pt(q=80, df=4, ncp=70)
p2 <- ptOwen(q=80, nu=4, delta=70)
wolfram <- 0.54742763380700947685
p1 - wolfram
## [1] 0.02268711
p2 - wolfram
## [1] 3.330669e-16

Limitations

When q=delta, the value of ptOwen(q, nu, delta) should go to 0.5 as nu increases to infinity. The examples below show the failure of this expectation when nu is too large.

ptOwen(q=50, nu=3500, delta=50)
## [1] 0.4986697
ptOwen(q=50, nu=3600, delta=50)
## [1] 0.4989289
ptOwen(q=50, nu=3650, delta=50)
## [1] 0.4522949
ptOwen(q=50, nu=3660, delta=50)
## [1] 0.3349762
ptOwen(q=50, nu=3670, delta=50)
## [1] 0.7607795
ptOwen(q=50, nu=3680, delta=50)
## [1] 0

Since all Owen’s algorithms are somehow similar to the algorithm evaluating ptOwen, we can suspect that the other ones also suffer from certain limitations. In the other vignette, we investigate the reliability and the limitations of OwenQ.

In order to do some comparisons, and thanks to the BH package, we have ported the boost implementation of the cumulative Student distribution to OwenQ. It is evaluated by the internal function pt_boost. We concluded that this function is highly reliable. In particular it does not suffer from the failures of ptOwen we have just shown:

OwenQ:::pt_boost(q=50, nu=3500, delta=50)
## [1] 0.4986697
OwenQ:::pt_boost(q=50, nu=3600, delta=50)
## [1] 0.4987041
OwenQ:::pt_boost(q=50, nu=3650, delta=50)
## [1] 0.4987206
OwenQ:::pt_boost(q=50, nu=3660, delta=50)
## [1] 0.4987238
OwenQ:::pt_boost(q=50, nu=3670, delta=50)
## [1] 0.4987271
OwenQ:::pt_boost(q=50, nu=3680, delta=50)
## [1] 0.4987303

Application to equivalence testing

The Owen distribution intervenes in the calculation of the power of equivalence tests.

One sample

Assume a statistical model given by a sample $y_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2)$ for i = 1, …, n. We want to demonstrate that, up to a given confidence level, the mean μ belongs to a certain interval [Δ1, Δ2]. In other words, we are interested in the alternative hypothesis H1: {Δ1 ≤ μ ≤ Δ2}.

Consider the usual 100(1 − 2α)%-confidence interval about μ: $$ \left[\bar y - t^\ast_{n-1}(\alpha)\frac{\hat\sigma}{\sqrt{n}}, \, \bar y + t^\ast_{n-1}(\alpha)\frac{\hat\sigma}{\sqrt{n}} \right]. $$

The H1 hypothesis is accepted at level α when this interval falls into the interval [Δ1, Δ2].

This can be written as follows: $$ T_1 := \frac{\bar y - \Delta_1}{\hat\sigma/\sqrt{n}} \geq t^\ast_{n-1}(\alpha) \quad \text{and} \quad T_2 := \frac{\bar y - \Delta_2}{\hat\sigma/\sqrt{n}} \leq - t^\ast_{n-1}(\alpha). $$

Observe that $$ T_1 = \frac{z - \delta_1}{\dfrac{\sqrt{n-1}\hat\sigma/\sigma}{\sqrt{n-1}}} $$ where $$ z = \frac{\sqrt{n}}{\sigma}(\bar y - \mu) \sim {\cal N}(0,1) \quad \text{and} \quad \delta_1 = \frac{\mu - \Delta_1}{\sigma/\sqrt{n}}. $$

By reasoning in the same way for T2, we find that the pair (T1, T2) follows the Owen distribution with degrees of freedom ν = n − 1, and non-centrality parameters δ1 given above and $\delta_2 = \tfrac{\mu - \Delta_2}{\sigma/\sqrt{n}}$.

Therefore the power of the test - i.e. the probability to accept H1 - is given by O4(n − 1, tn − 1*(α), −tn − 1*(α), δ1, δ2), and then can be evaluated thanks to powen4.

Inconclusive equivalence test

The result of the equivalence test is said to be inconclusive when only one of the bounds of the confidence interval falls into the interval [Δ1, Δ2].

The probability to get an inconclusive result can be obtained with the OwenQ package. We show and check that below with the help of simulations.

Delta1 <- -2; Delta2 <- 2 
mu <- 1; sigma <- 6; n <- 30L
alpha <- 0.05
nsims <- 1e6L
equivalence <- inconclusive <- numeric(nsims)
for (i in 1L:nsims) {
  y <- rnorm(n, mu, sigma)
  CI <- t.test(x = y, conf.level = 1-2*alpha)$conf.int
  equivalence[i] <- (CI[1] > Delta1) && (CI[2] < Delta2)
  inconclusive[i] <- ((CI[1] < Delta1) && (CI[2] > Delta1)) ||
    ((CI[1] < Delta2) && (CI[2] > Delta2))
}
dof <- n-1
q <- qt(1-alpha, dof)
se <- sqrt(1/n)*sigma
delta1 <- (mu-Delta1)/se; delta2 <- (mu-Delta2)/se
# probability to get equivalence
mean(equivalence)
## [1] 0.093195
powen4(dof, q, -q, delta1, delta2)
## [1] 0.09300963
# probability to get inconclusive
mean(inconclusive)
## [1] 0.901164
ptOwen(q, dof, delta2) - ptOwen(-q, dof, delta1) - powen4(dof, q, -q, delta1, delta2)
## [1] 0.90139

Parallel design

Now consider two independent samples $x_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2)$ for i = 1, …, n1. and $y_i \sim_{\text{iid}} {\cal N}(\nu, \sigma^2)$ for i = 1, …, n2 and the alternative hypothesis H1: {|μ − ν| ≤ Δ}. The power of this test is evaluated by the function powerTOST below. The parameter delta0 corresponds to the difference μ − ν.

powerTOST <- function(alpha, delta0, Delta, sigma, n1, n2) {
  se <- sqrt(1/n1 + 1/n2) * sigma
  delta1 <- (delta0 + Delta) / se
  delta2 <- (delta0 - Delta) / se
  dof <- n1 + n2 - 2
  q <- qt(1 - alpha, dof)
  powen4(dof, q, -q, delta1, delta2)
}

The Owen T-function

The OwenQ package also provides an implementation of the Owen T-function, under the name OwenT. This is a port of the function owens_t of boost, the peer-reviewed collection of C++ libraries.

The Owen Q-functions

The Owen cumulative functions are based on the two Owen Q-functions $$ 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 $$ and $$ Q_2(\nu, t, \delta, R) = \frac{1}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac12(\nu-2)}} \int_R^\infty \Phi\left(\frac{tx}{\sqrt{\nu}}-\delta\right) x^{\nu-1} e^{-\frac{x^2}{2}} \mathrm{d}x. $$

They are implemented in the OwenQ package under the respective names OwenQ1 and OwenQ2 (for integer values of ν).

Application: Equal-tailed tolerance interval

Consider the statistical model given by a sample $y_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2)$ for i = 1, …, n.

An equal-tailed (p, α)-tolerance interval is an interval containing at least 100p% of the “center data” with 100(1 − α)% confidence.

The natural choice for such an interval is, as shown by Owen (1965), [ȳ − keσ̂, ȳ + keσ̂] where ke satisfies $$ O_2(n-1, k_e\sqrt{n}, -k_e\sqrt{n}, \delta, -\delta) = 1-\alpha $$ where $\delta = \sqrt{n}z_{(1+p)/2}$.

Therefore, the tolerance factor ke can be determined with the help of the powen2 function of the OwenQ package. But it is more efficient to use the function spowen2 for this purpose; spowen2(nu, t, delta) has the same value as powen2(nu, t, -t, delta, -delta) but it is evaluated more efficiently.

p <- 0.9; alpha <- 0.05
n <- 100
delta <- sqrt(n)*qnorm((1+p)/2)
uniroot(function(ke) spowen2(n-1, ke*sqrt(n), delta) - (1-alpha), 
        lower=qnorm(1-alpha), upper=5, extendInt = "upX", tol=1e-9)$root
## [1] 1.981513

The ke factors are tabulated in Krishnamoorthy & Mathew’s book.

Internal functions: numerical integration

The OwenQ package provides three internal functions, iOwenQ1, iOwenQ2 and ipowen4, which respectively perform the evaluation of Q1, Q2 and O4 by numerical integration using the RcppNumerical package. They can be called with a positive non-integer value of ν. The evaluation of O4 with ipowen4 is correct only for t1 > t2 and δ1 > δ2.

References

  • Owen, D. B. (1965). A special case of a bivariate noncentral t-distribution. Biometrika 52, 437-446.

  • Krishnamoorthy & Mathew (2009). Statistical Tolerance Regions. Wiley.