While researching a topic in probability, I encountered a special mathematical function that I had not previously known about. The function is called Owen's T function (Owen, 1956, Ann. Math. Stat.). It is useful for computing multivariate probabilities, for defining the skew-normal distribution, and for computing probabilities associated with certain hypothesis tests. This article shows how to compute and visualize Owen's T function in SAS. I will show applications of Owen's T function in a subsequent article.

Owen's T function is defined as a one-dimensional integral, so this article uses the QUAD subroutine in SAS IML to evaluate it. However, it is worth mentioning that Patefield, M. and Tandy, D. (2000, JSS) provide a faster and more accurate algorithm. Their algorithm is more complex. It evaluates the function by selecting from among six different numerical methods, depending on the specified values of (x, α).

The definition of Owen's T function

The definition of Owen's T function, T(x, α), uses α as the upper limit of an integral and X as a parameter to the integrand:
$T(x, \alpha) = \frac{1}{2\pi} \int_0^\alpha \frac{e^{-\frac{x^2}{2}(1+t^2)}}{1+t^2} dt$

From this definition, it is clear that T is an even function in x and an odd function in α:

• T is symmetric in x: T(-x, α) = T(x, α)
• T is anti-symmetric in α: T(x, -α) = -T(x, α)

It is less obvious that Owen's T function simplifies if x=0 or α=1:

• T(0,α) = 1/(2 π) tan^-1(α)
• T(x,1) = 1/2 * Φ(x) (1 - Φ(x)), where Φ(x) is the standard univariate normal CDF.

Define Owen's T function in SAS

The QUAD subroutine in SAS IML is a general-purpose integration routine. The domain of integration in the formula is [0, α]. However, α can be negative, so it might be necessary to reverse the limits of integration. To handle that possibility, I specify the domain of integration as [0, |α|] and negate the answer if α < 0. Here is one way to define Owen's T function:

proc iml;
/* Define the Integrand for Owen's T function:
   The QUAD subroutine passes a scalar value, t, during the integration process.
   Use the GLOBAL statement to pass the parameter value, x (via 'g_x').
*/
start OwenT_Integrand(t) global(g_x);
   t2 = 1 + t##2;
   return ( exp(-0.5 * g_x##2 * t2) / t2 );
finish;
 
/* Define Owen's T Function by integrating OwenT_Integrand on [0, alpha].
   Loop over the elements of x to allow vector inputs. alpha is a scalar.
*/
start OwenT(x, alpha) global(g_x);
   pi = constant("pi");
   n = nrow(x);
   y = j(nrow(x), 1, .);
   /* make sure limits of integration are increasing. If alpha < 0, reverse sign of answer */
   limits = 0 || abs(alpha);
   do i = 1 to n;
      g_x = x[i];   /* set parameter, 'x' */
      call quad(val, "OwenT_Integrand", limits);
      y[i] = val;
   end;
   return( choose(alpha>=0, y/(2*pi), -y/(2*pi)) ); /* negate answer if alpha < 0 */
finish;
 
/* test evaluating the function for x in [-3,3] */
x = T(do(-3, 3, 6/50));
alpha = 1;
T = OwenT(x, alpha);
title "Graph of Owen's T Function: alpha=1";
call series(x, T) grid={x y};

The plot shows the graph of Owen's T function when the parameter α is fixed at α=1 and the parameter x is in the interval [-3, 3]. You can see that the graph looks somewhat bell-shaped.

The following panel of graphs shows the T(x, α) function for various values of α. You can see that the shape of the function is similar for all α values, but the amplitude of the function depends on α. You can show analytically that the magnitude of the T function is greatest at x=0. The height of the function at x=0 is 1/(2 π) tan^-1(α).

If you like 3-D graphs, the Boost C++ library contains a graph of T(x, α) over a rectangular domain.

Summary

Owen's T function is defined as a one-dimensional integral with two parameters. This article uses the QUAD subroutine in SAS IML to evaluate the integral for any parameter values. A subsequent article will discuss applications of Owen's T function in probability and statistics.