Computation of the expected value of a function of a chi-distributed random variable

TL;DR

This paper introduces an efficient numerical method using a transformation and trapezoidal rule for accurately computing the expected value of functions of chi-distributed variables, relevant in statistical inference and confidence region analysis.

Contribution

It presents a novel application of Mori's transformation combined with exponentially convergent trapezoidal rule and an error bound for evaluating expectations involving chi-distributions.

Findings

Method achieves exponential convergence for suitable integrands.

Provides a new upper bound for the infinite sum approximation error.

Demonstrates suitability for computing coverage probabilities and volumes of confidence regions.

Abstract

We consider the problem of numerically evaluating the expected value of a smooth bounded function of a chi-distributed random variable, divided by the square root of the number of degrees of freedom. This problem arises in the contexts of simultaneous inference, the selection and ranking of populations and in the evaluation of multivariate t probabilities. It also arises in the assessment of the coverage probability and expected volume properties of the some non-standard confidence regions. We use a transformation put forward by Mori, followed by the application of the trapezoidal rule. This rule has the remarkable property that, for suitable integrands, it is exponentially convergent. We use it to create a nested sequence of quadrature rules, for the estimation of the approximation error, so that previous evaluations of the integrand are not wasted. The application of the trapezoidal rule requires the approximation of an infinite sum by a finite sum. We provide a new easily computed upper bound on the error of this approximation. Our overall conclusion is that this method is a very suitable candidate for the computation of the coverage and expected volume properties of non-standard confidence regions.

Tables

Table 1. Table 1: The approximation error for the simple procedure described in Section 3 , with ϵ = 10 − 17 italic-ϵ superscript 10 17 \epsilon=10^{-17} and stopped after the computation of the approximation for n = 65 𝑛 65 n=65 for ν = 1 𝜈 1 \nu=1 and n = 33 𝑛 33 n=33 for ν = 2 , 3 , 4 , 5 , 10 , 100 𝜈 2 3 4 5 10 100 \nu=2,3,4,5,10,100 and 1000. Here α = 0.10 , 0.05 𝛼 0.10 0.05 \alpha=0.10,0.05 and 0.02 and ν = 1 , 2 , 3 , 4 , 5 , 10 , 100 𝜈 1 2 3 4 5 10 100 \nu=1,2,3,4,5,10,100 and 1000. We interpret an entry 0 in this table as | approximation error | < 1.11 × 10 − 16 approximation error 1.11 superscript 10 16 |\text{approximation error}|<1.11\times 10^{-16} .

	$\nu =1$	$\nu =2$	$\nu =3$	$\nu =4$	$\nu =5$
$\alpha =0.10$	$-1.11\times {10}^{-16}$	$-2.22\times {10}^{-16}$	$0$	$-2.22\times {10}^{-16}$	$-1.11\times {10}^{-16}$
$\alpha =0.05$	$9.66\times {10}^{-15}$	$2.10\times {10}^{-14}$	$-1.11\times {10}^{-16}$	$-1.11\times {10}^{-16}$	$-1.11\times {10}^{-16}$
$\alpha =0.02$	$-1.23\times {10}^{-12}$	$5.82\times {10}^{-11}$	$9.99\times {10}^{-16}$	$-1.11\times {10}^{-16}$	$0$

Table 2. Table 2: The approximation error for Generalized Gauss Laguerre quadrature for α = 0.10 , 0.05 𝛼 0.10 0.05 \alpha=0.10,0.05 and 0.02 and ν = 1 , 2 , 3 , 4 , 5 , 6 , 10 , 100 𝜈 1 2 3 4 5 6 10 100 \nu=1,2,3,4,5,6,10,100 and 300. The number of nodes m 𝑚 m is 65 65 65 for ν = 1 𝜈 1 \nu=1 and 33 33 33 for ν = 2 , 3 , 4 , 5 , 6 , 10 , 100 𝜈 2 3 4 5 6 10 100 \nu=2,3,4,5,6,10,100 and 300.

	$\nu =1$	$\nu =2$	$\nu =3$	$\nu =4$	$\nu =5$
$\alpha =0.10$	$1.44\times {10}^{-2}$	$1.32\times {10}^{-3}$	$1.63\times {10}^{-4}$	$2.86\times {10}^{-5}$	$6.08\times {10}^{-6}$
$\alpha =0.05$	$3.25\times {10}^{-2}$	$2.04\times {10}^{-3}$	$2.26\times {10}^{-4}$	$3.77\times {10}^{-5}$	$7.84\times {10}^{-6}$
$\alpha =0.02$	$2.00\times {10}^{-2}$	$4.12\times {10}^{-3}$	$3.39\times {10}^{-4}$	$5.24\times {10}^{-5}$	$1.05\times {10}^{-5}$

Table 3. Table 3: The approximation error for the inverse cdf method, using Gauss Legendre quadrature, for α = 0.10 , 0.05 𝛼 0.10 0.05 \alpha=0.10,0.05 and 0.02 and ν = 1 , 2 , 3 , 4 , 5 , 6 , 10 , 100 𝜈 1 2 3 4 5 6 10 100 \nu=1,2,3,4,5,6,10,100 and 1000. The number of nodes m 𝑚 m is 65 65 65 for ν = 1 𝜈 1 \nu=1 and 33 33 33 for ν = 2 , 3 , 4 , 5 , 6 , 10 , 100 𝜈 2 3 4 5 6 10 100 \nu=2,3,4,5,6,10,100 and 1000.

	$\nu =1$	$\nu =2$	$\nu =3$	$\nu =4$	$\nu =5$
$\alpha =0.10$	$7.77\times {10}^{-16}$	$6.39\times {10}^{-6}$	$1.52\times {10}^{-5}$	$2.07\times {10}^{-5}$	$2.37\times {10}^{-5}$
$\alpha =0.05$	$8.88\times {10}^{-16}$	$9.43\times {10}^{-6}$	$2.06\times {10}^{-5}$	$2.70\times {10}^{-5}$	$2.96\times {10}^{-5}$
$\alpha =0.02$	$8.88\times {10}^{-16}$	$1.53\times {10}^{-5}$	$2.94\times {10}^{-5}$	$3.58\times {10}^{-5}$	$3.75\times {10}^{-5}$

Table 4. Table 4: The approximation error for Gauss Legendre quadrature applied to the evaluation of ( 17 ), for α = 0.10 , 0.05 𝛼 0.10 0.05 \alpha=0.10,0.05 and 0.02 and ν = 1 , 2 , 3 , 4 , 5 , 10 , 100 𝜈 1 2 3 4 5 10 100 \nu=1,2,3,4,5,10,100 and 1000. The number of nodes m 𝑚 m is 65 65 65 for ν = 1 𝜈 1 \nu=1 and 33 33 33 for ν = 2 , 3 , 4 , 5 , 10 , 100 𝜈 2 3 4 5 10 100 \nu=2,3,4,5,10,100 and 1000.

	$\nu =1$	$\nu =2$	$\nu =3$	$\nu =4$	$\nu =5$
$\alpha =0.10$	$-2.84\times {10}^{-13}$	$-2.02\times {10}^{-10}$	$-1.26\times {10}^{-13}$	$-2.47\times {10}^{-13}$	$-6.75\times {10}^{-14}$
$\alpha =0.05$	$-1.67\times {10}^{-10}$	$-2.01\times {10}^{-8}$	$-6.44\times {10}^{-11}$	$8.88\times {10}^{-13}$	$1.47\times {10}^{-12}$
$\alpha =0.02$	$2.20\times {10}^{-7}$	$-5.48\times {10}^{-7}$	$-3.31\times {10}^{-9}$	$1.70\times {10}^{-10}$	$1.15\times {10}^{-12}$

Equations

$${\int }_0^{\infty }a(x)f_{\nu }(x)dx,$$

$${\int }_{-\infty }^{\infty }g(y)dy,$$

$$G(\omega )={\int }_{-\infty }^{\infty }g(y)\exp (-i\omega y)dy,$$

$$\left|h\sum _{j=-\infty }^{\infty }g(jh+\delta )-{\int }_{-\infty }^{\infty }g(y)dy\right|\le 2\sum _{j=1}^{\infty }\left|G\left(\frac{2\pi j}{h}\right)\right|,$$

$$h\sum _{j=-\infty }^{\infty }g(jh+\delta )$$

$$h\sum _{j=M}^{N}g(jh+\delta ),$$

$$f_{\nu}(x)=\begin{cases}\tau_{\nu}\,x^{\nu-1}\,\exp\big{(}-\nu\,x^{2}/2\big{)}&\text{for}\ \ x>0\\ 0&\text{otherwise,}\end{cases}$$

$${\tau }_{\nu }=\frac{{\nu }^{\nu /2}}{\Gamma (\nu /2)2^{(\nu /2)-1}}.$$

$$x(y)=\exp \left(\frac{1}{2}y-e^{-y}\right),$$

$$\frac{dx(y)}{dy}=\exp \left(\frac{1}{2}y-e^{-y}\right)\left(\frac{1}{2}+e^{-y}\right)$$

$$\int_{0}^{\infty}a(x)\,f_{\nu}(x)\,dx=\int_{-\infty}^{\infty}a\big{(}x(y)\big{)}\,\psi_{\nu}(y)\,dy,$$

$$\psi_{\nu}(y)=f_{\nu}\big{(}x(y)\big{)}\frac{dx(y)}{dy}.$$

$$|\psi_{\nu}(y)|\sim c_{1}\exp\big{(}-c_{2}\exp(c_{3}|y|)\big{)},\ y\rightarrow\pm\infty.$$

$$|G_{\nu}(\omega)|\leq c_{4}\exp\big{(}-c_{FT}|\omega|\big{)}$$

$$h\,\sum_{j=0}^{n-1}a\big{(}x(y_{\ell}+hj)\big{)}\,\psi_{\nu}(y_{\ell}+hj),$$

$$u_{\nu}(y,d)=Q_{\nu}\Big{(}\nu\,x^{2}(y)\Big{)}+1-Q_{\nu}\Big{(}\nu\,x^{2}\big{(}y+d\big{)}\Big{)}$$

$$\left|h\,\sum_{j=-\infty}^{-1}a\big{(}x(y_{\ell}+hj)\big{)}\,\psi_{\nu}(y_{\ell}+hj)\,+\,h\,\sum_{j=n}^{\infty}a\big{(}x(y_{\ell}+hj)\big{)}\,\psi_{\nu}(y_{\ell}+hj)\right|$$

$$h\sum _{j=-\infty }^{-1}{\psi }_{\nu }(y_{\ell }+hj)+h\sum _{j=n}^{\infty }{\psi }_{\nu }(y_{\ell }+hj),$$

$$h\sum _{j=n}^{\infty }{\psi }_{\nu }(y_{\ell }+hj)=h\sum _{j=1}^{\infty }{\psi }_{\nu }(y_u+hj),$$

$$h\sum _{j=1}^{\infty }{\psi }_{\nu }(y_u+hj)\le {\int }_{y_u}^{\infty }{\psi }_{\nu }(t)dt$$

$$\displaystyle=P\big{(}R>\nu^{1/2}x(y_{u})\big{)},$$

$$=1-Q_{\nu }(\nu x^2(y_u)).$$

$$h\,\sum_{j=-\infty}^{-1}\,\psi_{\nu}(y_{\ell}+hj)\leq\int_{-\infty}^{y_{\ell}}\psi_{\nu}(t)\,dt=Q_{\nu}\big{(}\nu\,x^{2}(y_{\ell})\big{)}.$$

$$\min_{y}u_{\nu}\big{(}y,d\big{)}=10^{-3}\,\epsilon.$$

$$\underset{y}{\min }{u}_{\nu }\left(y,d\right)\text{is equal to some specified small positive number}.$$

$$\frac{10{\tau }_{\nu }}{9\nu }\left(\exp \left(-\frac{\nu }{2}\exp \left(\frac{9y_{u0}}{10}\right)2^{c_{T}k}\right)+\exp \left(-\nu \exp \left(-\frac{9y_{\ell 0}}{10}\right)2^{c_{T}k}\right)\right)$$

$$+2c_4\exp \left(-\left(\frac{2\pi c_{FT}}{h_0}\right)2^k\right),$$

$${\int }_{y_{u0}+k{h}_0}^{\infty }{\psi }_{\nu }(t)dt+{\int }_{-\infty }^{y_{\ell 0}-k{h}_0}{\psi }_{\nu }(t)dt.$$

$${\psi }_{\nu }(t)\le \frac{{\tau }_{\nu }}{2}\exp \left(-\frac{\nu }{2}\exp \left(\frac{9}{10}t\right)\right)\text{for all}t\ge t_1$$

$${\psi }_{\nu }(t)\le {\tau }_{\nu }\exp \left(-\nu \exp \left(-\frac{9}{10}t\right)\right)\text{for all}t\le t_2.$$

$${\int }_y^{\infty }{\psi }_{\nu }(t)dt\le \frac{10{\tau }_{\nu }}{9\nu }\exp \left(-\frac{\nu }{2}\exp \left(\frac{9}{10}y\right)\right)\text{for all}y\ge t_1$$

$${\int }_{-\infty }^y{\psi }_{\nu }(t)dt\le \frac{10{\tau }_{\nu }}{9\nu }\exp \left(-\nu \exp \left(-\frac{9}{10}y\right)\right)\text{for all}y\le t_2.$$

Full text

Keywords: Confidence interval; Coverage probability; Mixed rule transformation; Numerical integration; Trapezoidal rule

1 Introduction

Consider the problem of finding an accurate and efficient method of numerically computing an integral of the form

[TABLE]

where $a$ is a smooth bounded real-valued function, $\nu$ is a positive integer and $f_{\nu}$ is the probability density function (pdf) of a random variable with the same distribution as $R/\nu^{1/2}$, where $R$ has a $\chi_{\nu}$ distribution (i.e. $R^{2}$ has a $\chi_{\nu}^{2}$ distribution). Note that \eqref{FormofExpectedValue}=E\big{(}a(R/{\nu}^{1/2})\big{)}, which is the expected value of a smooth bounded function of $R/{\nu}^{1/2}$. We suppose that a computer program for the accurate and efficient evaluation of $a(x)$, for any given $x>0$, is either already available or can be easily written. In other words, our focus is solely on the numerical evaluation of the integral (1).

The evaluation of an integral of this form occurs in the context of simultaneous statistical inference and the selection and ranking of populations (Miller 1981; Hochberg and Tamhane 1987; Gupta and Panchapakesan 2002) and in the evaluation of central and non-central (Kshirsagar definition) multivariate t probabilities (Dunnett and Sobel 1955; Dunnett 1989; Genz and Bretz 2009), when the method of Miwa et al. (2003), briefly described in Mi et al. (2009), is used to compute $a(x)$.

The evaluation of an integral of the form $\eqref{FormofExpectedValue}$ also occurs in the computation of the coverage probabilities of post-model-selection confidence intervals, frequentist model averaged confidence intervals and other non-standard confidence regions (Farchione and Kabaila 2008; Kabaila and Giri, 2009a; Kabaila and Giri, 2009b; Kabaila and Farchione 2012; Kabaila et al. 2016; Kabaila et al. 2017; Abeysekera and Kabaila 2017; Kabaila and Giri, 2013; Kabaila and Tissera, 2014; Kabaila 2018). In all of these papers, this evaluation has previously been carried out by first truncating the integral (the truncation error is easily bounded) and then applying an adaptive numerical integration method.

Our search for a better method for the evaluation of an integral of the form $\eqref{FormofExpectedValue}$ has led us to seek out an appropriate transformation of the variable of integration, followed by the application trapezoidal rule over the real line. As noted by Schwartz (1969), “The real artistry of numerical integration lies in learning to make changes of the variable” appropriate for the problem at hand and that “this must be studied separately for every problem”. The literature on various initial changes of the variable of integration for the purpose of efficient numerical integration is very large, with early references including Davis and Rabinowitz (1984), Sag and Szekeres (1964) and Imhof (1963). Some simple illustrations of the power of appropriate changes variable of integration prior to numerical integration are provided by Avery and Soler (1988).

We use the transformation (2.6) of Mori (1988), followed by application of the trapezoidal rule. This transformation belongs to a family of transformations proposed and investigated by Takahasi and Mori (1973), Mori (1985) and others. The trapezoidal rule has the remarkable property that, for suitable integrands, it is exponentially convergent (Trefethen and Weideman, 2014). There are several well-known explanations for this remarkable property, including the Euler-Maclaurin summation formula and Fourier transform methods. A historical review of the these explanations is provided in Section 11 of Trefethen and Weideman (2014). The trapezoidal rule also has the great advantage that it can be used to create a nested sequence of quadrature rules, used for the estimation of the approximation error, so that previous evaluations of the function $a$ are not wasted.

For our purposes, the best description of the properties of the trapezoidal rule is found using the Fourier transform of the integrand and the Poisson summation formula. For the reader’s convenience, this well-known description is recounted in Section 2. The application of the transformation (2.6) of Mori (1988) to the integral $\eqref{FormofExpectedValue}$, followed by the application of the trapezoidal rule is described in Section 3. In this section, we describe a method of carrying out the required ‘trimming’ of the infinite sum approximation to the integral that leads to an easily-computed upper bound on the resulting trimming error. In subsection 3.1, we describe a simple and effective procedure, similar to that described by Mori (1988, pp.370–371), for evaluating the integral (1) that leads to a nested sequence of quadrature rules. In subsection 3.2, we describe an extension of this procedure that we prove to be exponentially convergent under the appropriate regularity condition.

In Section 4 we use the simple test scenario that consists of evaluating a known univariate t probability (i.e the value of $\eqref{FormofExpectedValue}$ is known). We compare the performance of the method described in subection 3.1 with the following two methods:

Generalized Gauss Laguerre quadrature

Change the variable of integration in (1) to $y=\nu\,x^{2}/2$. In effect, we express the expectation of interest, E\big{(}a(R/{\nu}^{1/2})\big{)}, as E\big{(}a(2^{1/2}\,V^{1/2}/\nu^{1/2})\big{)}, where $V=R^{2}/2$ has a gamma$(\nu/2,1)$ distribution. We then apply Generalized Gauss Laguerre quadrature. This method has been widely applied in the literature on simultaneous inference and the selection and ranking of populations. 2. 2.

Inverse cdf method

Change the variable of integration in (1) to $y=F_{\nu}(x)$, where $F_{\nu}$ denotes the cumulative distribution function (cdf) that corresponds to the pdf $f_{\nu}$. This transforms the integral $\eqref{FormofExpectedValue}$ into an integral over the interval $[0,1]$. In effect, we express the expectation of interest, E\big{(}a(R/{\nu}^{1/2})\big{)}, as E\big{(}a(F_{\nu}^{-1}(U))\big{)}, where $U=F_{\nu}(R/\nu^{1/2})$ is uniformly distributed on $(0,1)$. We then apply Gauss Legendre quadrature. This method has been applied to the evaluation of central and non-central (Kshirsagar definition) multivariate t probabilities.

The purpose of this comparison is to illustrate the factors that may lead to a relatively poor performance of these two methods. The computations for this paper were carried out using the R computer language.

Finally, in Section 5 we discuss the application of the procedures described in Section 3 to the computation of the coverage probability and scaled expected length of of post-model-selection and frequentist model averaged confidence intervals. We also consider the application of these procedures to the computation of the coverage probability and scaled expected volume of other non-standard confidence regions.

2 Properties of the trapezoidal rule found using the Fourier transform of the integrand

Suppose that we wish to evaluate

[TABLE]

where $g$ is a real-valued absolutely integrable function. Let $G$ denote that Fourier transform of $g$. This transform is defined by

[TABLE]

where $i=\sqrt{-1}$ and the angular frequency $\omega\in\mathbb{R}$. Since $g$ is real-valued, $G(\omega)$ is an even function of $\omega$ (see e.g. p.11 of Papoulis 1962). It follows from the Poisson summation formula (see e.g. p.47 of Papoulis 1962) that

[TABLE]

for all $\delta\in[0,h)$. The left-hand side is the discretization error. This error is small when $|G(\omega)|$ decays rapidly as $\omega\rightarrow\infty$ and $h$ is sufficiently small.

We approximate the infinite sum

[TABLE]

by the finite sum

[TABLE]

for appropriately chosen integers $M$ and $N$ ($M<N$). The “trapezoidal rule” approximation to (2) is (5). The absolute value of the difference (5) $-$ (4) is called the trimming error. For (5) to be a good approximation to (2), we require that both the discretization error and the trimming error are small.

3 Application of the the transformation (2.6) of Mori (1988),

followed by the application of the trapezoidal rule

The pdf $f_{\nu}$ is given by

[TABLE]

where

[TABLE]

Throughout this section we suppose that $\nu$ is given. To evaluate (1), we first apply the transformation (2.6) of Mori (1988), namely

[TABLE]

so that

[TABLE]

and

[TABLE]

where

[TABLE]

As noted by Mori (1985) this transformation leads to $\psi_{\nu}(y)$ having double exponential decay as $y\rightarrow\pm\infty$, i.e. there exist positive numbers $c_{1}$, $c_{2}$ and $c_{3}$ such that

[TABLE]

This implies that g_{\nu}(y)=a\big{(}x(y)\big{)}\,\psi_{\nu}(y) also has double exponential decay as $y\rightarrow\pm\infty$. Computational results show that the function $\psi_{\nu}$ is unimodal for all positive integers $\nu$. Let $y_{\nu}^{*}$ denote the value of $y$ at which $\psi_{\nu}(y)$ is maximized. The value of $y_{\nu}^{*}$ is roughly 0.85 for all positive integers $\nu$. We suppose, without loss of generality, that $|a(x)|\leq 1$ for all $x\in\mathbb{R}$.

Let $G_{\nu}$ denote the Fourier transform of $g_{\nu}(y)$. We now introduce the following assumption:

Assumption FT:  There exist positive numbers $c_{4}$ and $c_{FT}$ such that

[TABLE]

for all $\omega\in\mathbb{R}$. In other words, $G_{\nu}(\omega)$ has single exponential decay as $\omega\rightarrow\pm\infty$.

Theorem 5.1 of Trefethen and Weideman (2014) provides conditions on the function $g_{\nu}(y)$ that imply that this assumption holds.

We will approximate (7) by

[TABLE]

where $n$ denotes the number of evaluations of the integrand a\big{(}x(y)\big{)}\psi_{\nu}(y), $h$ denotes the step length and the first evaluation of this integrand is at $y_{\ell}$. Let $d=nh$. Of course, our aim is to choose \big{(}n,h,y_{\ell}\big{)} such that (9) provides an efficient and accurate approximation.

We will use the following result which provides an easily computed upper bound on the trimming error.

Lemma 1. * Suppose that $y_{\ell}<y_{\nu}^{*}$ and that $y_{\ell}+d>y_{\nu}^{*}$. Then, when we approximate (7) by (9), the trimming error is bounded above by u_{\nu}\big{(}y_{\ell},d\big{)}, where*

[TABLE]

*and $Q_{\nu}$ denotes the $\chi_{\nu}^{2}$ cdf. *

Proof. Suppose that $y_{\ell}<y_{\nu}^{*}$ and that $y_{\ell}+d>y_{\nu}^{*}$. The trimming error is

[TABLE]

The trimming error is bounded above by

[TABLE]

since, for all positive integers $\nu$, $\psi_{\nu}(y)>0$ for all $y\in\mathbb{R}$. Observe that

[TABLE]

where $y_{u}=y_{\ell}+d$. We now use the same reasoning as for the integral test for series convergence. Since $\psi_{\nu}(y_{u}+t)$ is a decreasing function of $t\geq y_{\nu}^{*}$,

[TABLE]

Similarly, since $\psi_{\nu}(y_{\ell}+t)$ is an increasing function of $t\in(-\infty,y_{\nu}^{*}]$,

[TABLE]

Therefore (10) is bounded above by u_{\nu}\big{(}y_{\ell},d\big{)}. $\square$

3.1 A simple and effective procedure for evaluating the integral (7)

Suppose that we are given the value $\epsilon>0$ of a desired upper bound on the absolute value of the approximation error that we will develop. We now describe a simple and effective procedure for evaluating the integral (7), to roughly this accuracy, that leads to a nested sequence of quadrature rules. This procedure, which is similar to that described by Mori (1988, pp.370–371), consists of the following steps:

Step 1: Choose $y_{\ell}$ and $d$ and an initial value of $n$

The upper bound (3) on the discretization error suggests that, for a given value of the upper bound on the trimming error, as given in Lemma 1, it makes sense to minimize $h$. This provides the motivation for the following choice of $d$. Choose $d$ such that

[TABLE]

Choose $y_{\ell}$ to be the value of $y$ minimizing u_{\nu}\big{(}y,d\big{)}. This will ensure that the magnitude of the approximation error will be dominated by the discretization error. This is not as wasteful of evaluations of the integrand $g_{\nu}(y)$ as it might seem at first since $g_{\nu}(y)$ has double exponential decay as $y\rightarrow\pm\infty$. We have chosen the initial value of $n$ to be 5. Proceed to the next step.

Step 2: For given $(n,h,y_{\ell})$, evaluate the approximation (9)

Evaluate the approximation (9) and store the result. Using the stored values of the approximations decide whether or not to stop the procedure. Because the magnitude of the approximation error is dominated by the discretization error, this stopping rule can depend simply on estimating the discretization errors, as in the procedure described by Mori (1988, pp.370–371). Proceed to the next step.

Step 3: Halve $h$ and go back to the previous step

3.2 An exponentially convergent procedure for evaluating the integral (7)

While the procedure described in the previous subsection is simple to program and effective (as evidenced by the numerical results presented in Section 4), it does not lead to exponential convergence. We now describe a procedure that results in a nested sequence of quadrature rules that, under Assumption FT, is exponentially convergent. The fact that $g_{\nu}(y)$ has double exponential decay as $y\rightarrow\pm\infty$, whereas its Fourier transform $G_{\nu}(\omega)$ has only single exponential decay as $\omega\rightarrow\pm\infty$, implies that, at each iteration, $d$ should be increased at a slower rate than $1/h$. By adding a given positive number $2b$ to $d$ and halving $h$ at each iteration, we obtain exponential convergence. For simplicity of exposition, we have not included a stopping rule in the description of this procedure.

Step 1: An initial choice of a reasonable value of $(y_{\ell},n,d)$

Choose an initial value of $n$, which we denote by $n_{0}$. The initial value of $h$, denoted by $h_{0}$, is the initial value of $d$ (to be specified shortly) divided by $n_{0}$. We choose $b$ to be some small positive integer multiple of $h_{0}$. For the sake of concretness, we have chosen $b=h_{0}$. The initial value of $d$ is such that

[TABLE]

The initial value of $y_{\ell}$, denoted by $y_{\ell 0}$, is the value of $y$ minimizing $u_{\nu}(y,d)$ for the chosen initial value of $d$, denoted by $d_{0}$. Let $y_{u0}=y_{\ell 0}+d_{0}$. Proceed to the next step.

Step 2: For given $(y_{\ell},n,d)$, evaluate the approximation (9)

Evaluate the approximation (9) and store the result. Proceed to the next step.

Step 3: Add $2b$ to $d$, halve $h$ and choose the new value of $y_{\ell}$

Add $2b$ to $d$ and halve $h$. Choose the new value of $y_{\ell}$ to be $y_{\ell}-b$. It will be convenient for the proof of exponential convergence to define the iteration number $k$ by $h=h_{0}/2^{k}$. Go back to the previous step.

The following theorem states that under Assumption FT this procedure is exponentially convergent. The type of convergence described in this theorem is consistent with that other double exponential types of quadrature formulas (Mori and Sugihara, 2001).

Theorem 1. * Suppose that Assumption FT holds. Then the magnitude of the approximation error is, for all sufficiently large iteration numbers $k$, bounded above by*

[TABLE]

*where c_{T}=9\,h_{0}\big{/}\big{(}10\,\log_{e}(2)\big{)}. Since, at iteration number $k$, $n=(n_{0}+2\,k)\,2^{k}$, the magnitude of the approximation error converges exponentially to 0 as $n\rightarrow\infty$. *

Proof. Suppose that Assumption FT holds. By the proof of Lemma 1, the trimming error for iteration number $k$, is bounded above by

[TABLE]

It may be shown that there exist $t_{1}<\infty$ and $t_{2}>-\infty$ such that

[TABLE]

and

[TABLE]

It follows from this that

[TABLE]

and

[TABLE]

Therefore, for all sufficiently large iteration numbers $k$, (11) is bounded above by

[TABLE]

where c_{T}=9h_{0}\big{/}\big{(}10\,\log_{e}(2)\big{)}.

It follows from the upper bound (3) on the discretization error and Assumption FT that, for all sufficiently large iteration numbers $k$, the discretization error is bounded above by

[TABLE]

$\square$

4 Comparison with two other methods of numerical integration

In this section we use the simple test scenario that consists of evaluating a known univariate t probability. We compare the performance of the method described in the previous section with the two other methods described in the introduction.

Through the consideration of the coverage probability of a $1-\alpha$ t-interval, it may be shown that

[TABLE]

where

[TABLE]

with $\Phi$ the $N(0,1)$ cdf and the quantile $t_{\nu,a}$ defined by P\big{(}T\leq t_{\nu,a}\big{)}=a for $T\sim t_{\nu}$. Figure 1 provides an illustration of the fact that the $a_{\nu,\alpha}(x)$’s are smooth bounded functions of $x$ for the values of $\alpha$ considered and all positive integers $\nu$. This figure presents graphs of $a_{\nu,\alpha}(x)$ as a function of $x$ for $\alpha=0.05$ and $\nu=1,2$ and $\infty$. The graph labeled $\nu=\infty$ refers to the case that $t_{\nu,1-\alpha/2}$ is replaced by its limit, as $\nu\rightarrow\infty$.

4.1 The transformation

(2.6) of Mori (1988), followed by the application of the trapezoidal rule

Apply the transformation (2.6) of Mori (1988), so that

[TABLE]

where g_{\nu,\alpha}(y)=a_{\nu,\alpha}\big{(}x(y)\big{)}\,\psi_{\nu}(y).

We apply the simple procedure described in Section 3, with $\epsilon=10^{-17}$ and stopped after the computation of the approximation for $n=65$ for $\nu=1$ and $n=33$ for $\nu=2,3,4,5,10,100$ and 1000. The approximation error is defined to be this approximation minus $1-\alpha$. Table 1 presents the approximation error for $\alpha=0.10,0.05$ and 0.02 and $\nu=1,2,3,4,5,10,100$ and 1000. Due to the finite precision of our computations in R, we interpret an entry 0 in this table as $|\text{approximation error}|<1.11\times 10^{-16}$.

4.2 Generalized Gauss Laguerre quadrature

To apply Generalized Gauss Laguerre quadrature to the evaluation of (1), change the variable of integration to $y=\nu\,x^{2}/2$, so that

[TABLE]

where $c(y)=y^{(\nu/2)-1}\,\exp(-y)$ and d_{\nu}(y)=a\big{(}(2y/\nu)^{1/2}\big{)}. We then apply Generalized Gauss Laguerre quadrature, with $m$ nodes (samples), to approximate

[TABLE]

by

[TABLE]

for the appropriately chosen $w_{j}$’s (which are all positive) and $y_{j}$’s ($0<y_{1}<\dots<y_{m}<\infty$). We define the approximation error to be (14) minus $1-\alpha$.

Graphs of $d_{\nu}(y)$ as a function of $y$ are shown in Figure 2 for $\nu=1,2,3$ and 10 and $\alpha=0.10,0.05$ and 0.02. It should be noted that the horizontal scales in each of the four panels of this figure are very different. It is known that Generalized Gauss Laguerre quadrature with $m$ nodes will lead to the exact result if $d_{\nu}(y)$ is a polynomial in $y\in[0,\infty)$ of degree $2m-1$ (Chandrasekhar, 1960, p.65).

We assess how well $d_{\nu}(y)$ can be approximated by a polynomial, over the finite interval of values of $y$ such that $c(y)/\Gamma(\nu/2)$ is substantially greater than 0, as follows. Any polynomial $p$ of degree $u$ can be written as

[TABLE]

Set $a_{0}=1$ and require that $\sum_{j=1}^{u}a_{j}=1$, so that the functions $p$ and $d_{\nu}$ take the same values at both $y=0$ and $y=1$. A first approximation to $d_{\nu}(y)$ by $p(y)$ over the interval $y\in[0,1]$ is obtained by minimizing a measure of distance between $d_{\nu}(y)$ and $1-(1-y)^{j}$, over $j\in\{1,\dots,u\}$. A better approximation is obtained by minimizing a measure of distance between $d_{\nu}(y)$ and $1-\sum_{j=1}^{u}a_{j}(1-y)^{j}$, over $a_{1},\dots,a_{u}$, subject to $\sum_{j=1}^{u}a_{j}=1$. It follows from the shapes of the graphs in Figure 3 that to approximate $d_{\nu}(y)$ well by a polynomial, over the finite interval of values of $y$ such that $c(y)/\Gamma(\nu/2)$ is substantially greater than 0, we would require this polynomial to be of very high degree, particularly for small $\nu$. This suggests that Generalized Gauss Laguerre quadrature, with a given number of nodes $m$, will be most inaccurate for $\nu=1$ and will have increasing accuracy as $\nu$ increases.

This suggested result is borne out by Table 2 , which lists the approximation error for Generalized Gauss Laguerre quadrature for $\alpha=0.10,0.05$ and 0.02 and $\nu=1,2,3,4,5,6,10,100$ and 300. We have chosen the number of nodes $m$ to be the same as the number of integrand evaluations in Table 1. In other words, the number of nodes $m$ is $65$ for $\nu=1$ and $33$ for $\nu=2,3,4,5,6,10,100$ and 300.

Further confirmation of the unsuitability of Generalized Gauss Laguerre quadrature, in the scenario under consideration, for $\nu=1$ and $\nu=2$ is provided by Figure 3. The top and bottom panels of this figure are scatterplots of the $(y_{j},w_{j})$’s for $(\nu,m)=(1,65)$ and $(\nu,m)=(2,33)$, respectively ($y_{j}\leq 50$). For $(\nu,m)=(1,65)$ and $(\nu,m)=(2,33)$ there are 30 values of $y_{j}>50$ and 9 values of $y_{j}>50$, respectively. When we compare the top panel of Figure 3 with the top left panel (the case $\nu=1$) of Figure 2, we observe the following. Generalized Gauss Laguerre quadrature uses very few samples for the values of $y$ where the function $d_{\nu}(y)$ is changing rapidly with increasing $y$, while using a large number of samples for values of $y$ at which this function hardly changes with increasing $y$. Indeed, for $(\nu,m)=(1,65)$ there are only 2 nodes in the interval $[0,0.1]$. A similar conclusion results from comparing the bottom panel of Figure 3 with the top right panel (the case $\nu=2$) of Figure 2. For $(\nu,m)=(2,33)$ there are only 3 nodes in the interval $[0,1]$.

Of course, one could greatly increase the number of nodes $m$ and then approximate (14) by $\sum_{j=1}^{q}w_{j}\,d_{\nu}(y_{j})$, where $q$ is much less than $m$. This is unsatisfactory for the following two reasons. Firstly, the raison d’etre of Generalized Gauss Laguerre quadrature is that with $m$ nodes it leads to the exact result for polynomials of degree $2m-1$. This fundamental property is lost when this approximation is carried out. Secondly, this is a rather ad hoc way of forcing more samples of the function $d_{\nu}(y)$ into the quadrature formula for the values of $y$ for which this function changes rapidly with increasing $y$.

4.3 Inverse cdf method, using Gauss Legendre quadrature

Change the variable of integration to $y=F_{\nu}(x)$, where $F_{\nu}$ denotes the cdf corresponding to the pdf $f_{\nu}$, so that

[TABLE]

A similar transformation is used, for example, by Genz and Bretz (2009, p.32). If desired, we can compute $F_{\nu}^{-1}(y)$ using either F_{\nu}^{-1}(y)=\big{(}Q_{\nu}^{-1}(y)\big{/}\nu\big{)}^{1/2} or F_{\nu}^{-1}(y)=F_{R}^{-1}(y)\big{/}\nu^{1/2}, where $F_{R}$ denotes the $\chi_{\nu}$ cdf of $R$. We then change the variable of integration to $z=2y-1$ to obtain

[TABLE]

where b_{\nu}(z)=a\left(F_{\nu}^{-1}((z+1)/2)\right)\big{/}2. We then approximate the right-hand side, using Gauss Legendre quadrature with $m$ nodes, by

[TABLE]

for the appropriately chosen $\widetilde{w}_{j}$’s (which are all positive) and $z_{j}$’s ($-1<z_{1}<\dots<z_{m}<1$). We define the approximation error to be (15) minus $1-\alpha$.

Graphs of $b_{\nu}(z)$ as a function of $z$ are shown in Figure 4 for for $\nu=1,3,10$ and 100 and $\alpha=0.10,0.05$ and 0.02. It should be noted that the horizontal scale for the $\nu=1$ panel is different from the horizontal scale for the $\nu=3$, $\nu=10$ and $\nu=100$ panels (which are the same). It is known that Gauss Legendre quadrature with $m$ nodes will lead to the exact result if $b_{\nu}(z)$ is a polynomial in $z\in[-1,1]$ of degree $2m-1$. When interpreting Figure 4, it is important to remember that $b_{\nu}(-1)=0$ and that $b_{\nu}(z)$ is an increasing continuous function of $z\in[-1,1]$. It is evident, then, from this figure that $b_{\nu}(z)$ increases very rapidly as $z$ increases from zero for $\nu=10$ and $\nu=100$.

It follows from Figure 4 and the same kinds of considerations as in subsection 4.2 that the degree of the polynomial in $z$ needed to approximate $b_{\nu}(z)$ well in the interval $z\in[-1,1]$ increases with increasing $\nu$. This suggests that the inverse cdf method, using Gauss Legendre quadrature with a given number of nodes $m$, will be most accurate for $\nu=1$ and will have decreasing accuracy as $\nu$ increases.

This suggested result is borne out by the first 7 columns (the columns labelled $\nu=1$ to $\nu=10$) of Table 4, which lists the approximation error for Gauss Legendre quadrature for $\alpha=0.10,0.05$ and 0.02 and $\nu=1,2,3,4,5,6,10,100$ and 1000. We have chosen the number of nodes $m$ to be the same as the number of integrand evaluations in Table 1. In other words, the number of nodes $m$ is $65$ for $\nu=1$ and $33$ for $\nu=2,3,4,5,6,10,100$ and 1000.

5 Application to the computation of the coverage probability and scaled expected volume of non-standard confidence regions

To assess the coverage probability and expected volume properties of the non-standard confidence regions considered in the references co-authored with Kabaila, one needs to evaluate, for given $\nu$, integrals of the form (1) for hundreds, or thousands or even tens of thousands of different functions $a$. Each of these functions is smooth and bounded and the evaluation of $a(y)$ for any given $y$ is computationally expensive. In this case, the following “set-up costs” are negligible:

For the simple procedure described in Section 3 (the transformation (2.6) of Mori (1988), followed by application of the trapezoidal rule), the “set-up cost” is computing $y_{\ell}$ and $d$. 2. 2.

For Generalized Gauss Laguerre quadrature, the “set-up cost” is computing the weights, $w_{j}$’s, and nodes, $y_{j}$’s, for this quadrature, followed by the computation of the $(2y_{j}/\nu)^{1/2}$’s 3. 3.

For the Inverse cdf method, the “set-up cost” consists of computing the weights, $w_{j}$’s, and nodes, $z_{j}$’s, for Gauss Legendre quadrature, followed by the computation of the $F_{\nu}^{-1}((z_{j}+1)/2)$’s.

In other words, the number of evaluations of the function $a$ provides a reasonable guide to the computational effort for each of these methods.

We now consider in detail the evaluations of integrals of the form (1) in the references co-authored with Kabaila. Kabaila and Giri (2009a), Kabaila and Giri (2009b), Kabaila and Tissera (2014) and Abeysekera and Kabaila (2017) need to evaluate integrals of the form

[TABLE]

where $\xi$ and $\kappa$ are a positive integers and $\lambda:[0,\infty)\rightarrow\mathbb{R}$ is a smooth bounded function. This integral can be converted into the form (1) by changing the variable of integration to $y=c(\kappa,\xi)\,x$, where c(\kappa,\xi)=\big{(}\kappa/(\kappa+\xi)\big{)}^{1/2}, so that

[TABLE]

where $\nu=\kappa+\xi$ and a(y)=\lambda\big{(}y\big{/}c(\kappa,\xi)\big{)} is a smooth bounded function of $y\geq 0$.

An important measure of the performance of a confidence interval is its coverage probability function. The assessment of the coverage probability functions of (a) the post-model-selection confidence intervals considered by Kabaila and Giri (2009b) and Kabaila and Farchione (2012) and (b) the frequentist model averaged confidence intervals considered by Kabaila et al. (2016) is carried out by plotting the graphs of these functions. This requires the evaluation, for some given $\nu$, of an expression of the form (1) for, say, 200 different functions $a$. Each of these functions is smooth and bounded and the evaluation of $a(y)$ for any given $y$ is computationally epxensive.

Abeysekera and Kabaila (2017), Kabaila and Giri (2009a), Kabaila and Giri (2013) and Kabaila and Tissera (2014) construct non-standard confidence regions with guaranteed coverage using the following computations. They numerically optimize a criterion related to the expected volume of a parametric family of non-standard confidence regions, subject to a coverage probabililty constraint. The coverage probabililty, for a particular true parameter value, of a member of this family is given by an expression of the form (1), for some given $\nu$. The function $a$ is smooth and bounded and the evaluation of $a(y)$ for any given $y$ is computationally expensive. As this numerical constrained optimization proceeds, the evaluation of an expression of the form (1) needs to be carried out for thousands or even tens of thousands of different functions $a$, for the same given $\nu$.

6 Remarks

For the computation of $y_{\ell}$ and $d$ in Step 1 of our procedure, we have evaluated $Q_{\nu}$ using the R function pchisq. This evaluation of $Q_{\nu}$ is carried out using well-established methods for the evaluation of the incomplete gamma integral. These methods include the series expansion described by Shea (1988), as well as a continued fraction expansion due to Gauss, which greatly simplifies for $\nu$ an even positive integer. As already noted, for the types of problems considered in the references co-authored with Kabaila, the “set-up cost” of Step 1 is negligible. However, if one really needed to reduce the computation time for Step 1 then one could do so by replacing the exact evaluation of the tail probabilities of the $\chi_{\nu}^{2}$ distribution by upper bounds (such as Chernoff bounds) on these probabilities and by simplifying the minimization and root finding steps needed to evaluate $y_{\ell}$ and $d$.

A reviewer has suggested, following Step 1 of our procedure, the application of Gauss-Legendre quadrature (instead of the application of the trapezoidal rule) to the evaluation of the truncated integral

[TABLE]

This application is made in the usual way by first carrying out a straight line transformation of the interval $[y_{\ell},y_{u}]$ to $[-1,1]$. The resulting approximation errors are shown in Table 4. The magnitudes of these approximation errors are all larger than the magnitudes of the corresponding approximation errors for the trapezoidal rule, with the same number of integrand evaluations, reported in Table 1. In other words, for the same number of evaluations of the integrand, the trapezoidal rule outperforms Gauss-Legendre quadrature applied to the evaluation of (17), in terms of magnitude of approximation error. This result may be explained by the fact that the Gauss-Legendre quadrature nodes, which lie in the interval $[-1,1]$, cluster near the values $-1$ and 1, where the transformed integrand takes values very close to zero. This clustering also leads to Gauss-Legendre quadrature nodes not far from 0, where the transformed integrand differs most from zero, being more widely spaced than for the trapezoidal rule, with the same number of evaluations of the integrand.

An additional advantage of the trapezoidal rule is that, unlike Gauss-Legendre quadrature, it leads to a nested sequence of quadrature rules that can be used for the estimation of the approximation error.

7 Discussion

In Section 4, for both Generalized Gauss Laguerre quadrature and the Inverse cdf method, we present graphs whose features accurately predict their performance in terms of accuracy for a given number of evaluations of the function $a$. As noted in Section 5, the number of evaluations of the function $a$ is a reasonable measure of computational effort when the “set-up costs” are negligible, as in the situations considered in the references co-authored with Kabaila.

Our findings for the test scenario considered in Section 4 are as follows. The Generalized Gauss Laguerre quadrature method performs worst for $\nu=1$, and has performance that improves with increasing $\nu$. It has the worst performance of the three methods for $\nu\in\{1,2,3,4\}$. The inverse cdf method, using Gauss Legendre quadrature, has the best performance of the three methods $\nu=1$ and $\alpha\in\{0.05,0.02\}$, and has performance that decreases as $\nu$ increases through the values $2,3,4,5,6$ and 10. The method described in Section 3 (application of the the transformation (2.6) of Mori, 1988) has the best performance for $\nu=1$ and $\alpha=0.1$, $\nu\in\{1,2,3,4,5,10\}$, has very close to the best performance for $\nu=100$ and has the best performance for $\nu=1000$. For many of the situations considered in the references co-authored with Kabaila, the smallest possible value of $\nu$, in the evaluation of integrals of the form (1), is 2.

The procedures described in Section 3 use a nested sequence of quadrature rules, for the estimation of the approximation error, so that previous evaluations of the integrand are not wasted. This nested sequence can be implemented in a very simple computer program. This is an important advantage of this method over the other two methods.

Taken together, the results presented in this paper show that the simple procedure described in subsection 3.1 is a very suitable candidate for the computation of the coverage and expected volume properties of non-standard confidence regions considered in the references co-authored with Kabaila. The work presented in the paper is motivated by a need to compute these properties. However, it is clear that the application of transformations, such as those put forward by Schwartz (1969), Takahasi and Mori (1973) and Mori (1988), followed by the application of the trapezoidal rule will be useful in computing the expected values of functions of continuous random variables for a wide range of probability distributions of these random variables.

Acknowledgement

This work was supported by an Australian Government Research Training Program Scholarship.