Fitting the Distribution of Linear Combinations of t− Variables with more than 2 Degrees of Freedom

Abstract

The linear combination of Student’s $t$ random variables (RVs) appears in many statistical applications. Unfortunately, the Student’s $t$ distribution is not closed under convolution, thus, deriving an exact and general distribution for the linear combination of $K$ Student’s $t$ RVs is infeasible, which motivates a fitting/approximation approach. Here, we focus on the scenario where the only constraint is that the number of degrees of freedom of each $t-$ RV is greater than two. Notice that since the odd moments/cumulants of the Student’s $t$ distribution are zero and the even moments/cumulants do not exist when their order is greater than the number of degrees of freedom, it becomes impossible to use conventional approaches based on moments/cumulants of order one or higher than two. To circumvent this issue, herein we propose fitting such a distribution to that of a scaled Student’s $t$ RV by exploiting the second moment together with either the first absolute moment or the characteristic function (CF). For the fitting based on the absolute moment, we depart from the case of the linear combination of $K=2$ Student’s $t$ RVs and then generalize to $K\ge 2$ through a simple iterative procedure. Meanwhile, the CF-based fitting is direct, but its accuracy (measured in terms of the Bhattacharyya distance metric) depends on the CF parameter configuration, for which we propose a simple but accurate approach. We numerically show that the CF-based fitting usually outperforms the absolute moment-based fitting and that both the scale and number of degrees of freedom of the fitting distribution increase almost linearly with $K$.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The Student’s $t$-distribution arises in numerous scenarios, e.g., when estimating the mean of a normally distributed population of unknown variance with relatively few samples and in Bayesian analysis of data from a normal family. Moreover, such a distribution plays a key role in many relevant statistical analyses, including Student’s $t$-test for assessing the statistical significance of the difference between two sample means, the construction of confidence intervals for the difference between two population means, and linear regression analysis [1–6].

One of the distinctive properties of the Student’s $t$ distribution is its heavy tail. This behavior is also seen in the famous family of stable distributions; however, the Student’s $t$ distribution is more analytically tractable, which allows, for example, to write down explicitly its likelihood function [1].

1.1. Main Statistics of the Student’s $t$ Distribution

The probability density function (PDF) and cumulative density function (CDF) of a Student’s $t$ random variable (RV) $T$ with $\nu \in {\mathbb{R}}^{+}$ degrees of freedom, i.e., $T\sim \mathcal{T}\nu$, is given by [6].$$\begin{array}{}\text{(1)}& f_{T}x={\alpha }_{\nu }{1+\frac{x^2}{\nu }}^{-\nu +1/2}\text{,}\text{(2)}& F_{T}x=1-\frac{1}{2}{\mathcal{I}}_{\nu /x^2+\nu }\frac{\nu }{2},\frac{1}{2}\text{,}\end{array}$$$\forall x\in \mathbb{R}$, where ${\mathcal{I}}_x\cdot ,\cdot$ is regularized incomplete beta function ([7], Equation (8.17.2)), and$$\begin{array}{}\text{(3)}& {\alpha }_{\nu }\triangleq \frac{1}{\sqrt{\nu }B\nu /2,1/2}=\frac{\mathrm{\Gamma }\nu +1/2}{\mathrm{\Gamma }\nu /2\sqrt{\nu \pi }}\text{.}\end{array}$$here, $B\cdot ,\cdot$ and $\mathrm{\Gamma }\cdot$ are the beta ([7], Equation (5.12.1)) and gamma ([7], Equation (5.2.1)) functions, respectively. Meanwhile, the integer moments are [6] as follows:$$\begin{array}{}\text{(4)}& \mathbb{E}{T}^m={\nu }^{m/2}{\displaystyle \prod }_{i=1}^{m/2}\frac{2i-1}{\nu -2i},m\text{even},m<\nu \text{,}\end{array}$$while $\mathbb{E}{T}^m=0$ for $m$ odd, and moments of order $\nu$ or higher do not exist. We observe that for the specific case of $m=2$ (second moment), (4) reduces to $\nu /\nu -2$.

Two statistics that play a key role in our proposed approach are the absolute moments and the characteristic function (CF). The latter is given by [6]$$\begin{array}{}\text{(5)}& {\text{CF}}_{T}r=\mathbb{E}{e}^{\dot{\iota }rT}=\frac{{\sqrt{\nu }r}^{\nu /2}{K}_{\nu /2}\sqrt{\nu }r}{2^{\nu /2-1}\mathrm{\Gamma }\nu /2},\nu >0\text{,}\end{array}$$where $K_{\nu }\cdot$ is the modified Bessel function of second kind and order $\nu$ ([7], Sec. 10.25) (Notice that ${\text{CF}}_{T}r\in \mathbb{R}$, while the moment generating function of $T$ does not exist). Meanwhile, the absolute moments can be obtained as follows:$$\begin{array}{c}\text{(6)}& \mathbb{E}{T}^m={\int }_{-\infty }^{\infty }{x}^m{f}_{T}x\mathrm{d}x\stackrel{a}{=}2{\alpha }_{\nu }{\int }_0^{\infty }{x}^m{1+\frac{x^2}{\nu }}^{-\nu +1/2}\mathrm{d}x\\ \stackrel{b}{=}\frac{{\alpha }_{\nu }{\nu }^{m+1/2}\mathrm{\Gamma }\nu -m/2\mathrm{\Gamma }m+1/2}{\mathrm{\Gamma }\nu +1/2}\text{for\hspace{0.17em}}m<\nu \\ \stackrel{c}{=}\frac{{\nu }^{m/2}\mathrm{\Gamma }\nu -m/2\mathrm{\Gamma }m+1/2}{\sqrt{\pi }\mathrm{\Gamma }\nu /2},m<\nu \text{,}\end{array}$$where (a) follows from leveraging the symmetry of $T$ around zero and from substituting (1), (b) exploits ([8], Equation (3.241.4)) to solve the definite integral, and (c) is attained after simple algebraic transformations after substituting (3).

1.2. On the Linear Combination of Student’s $t$ RVs

The linear combination of Student’s $t$ RVs, denoted as follows:$$\begin{array}{}\text{(7)}& Z\triangleq {\displaystyle \sum }_{i=1}^K{\sigma }_i{T}_i\text{,}\end{array}$$where $T_i\sim \mathcal{T}{\nu }_i$ and (without loss of generality) ${\sigma }_i>0$, appears in many statistical applications. For instance,

(i) Fairweather [9] proposed a method based on the pivotal quantity $Z$ to obtain an accurate confidence interval for the common mean of several normal populations. Notice that the problem of characterizing the distribution of independent samples that are collected from different normal populations with a common mean but possibly with different variances appears in many practical applications, e.g., when different instruments/methods/laboratories are used to measure substances or products to assess their average quality [10];

(ii) The Behrens–Fisher distribution of the test statistic for testing the equality of the means of two normal populations with unknown variances is that of a linear combination of two independent Student’s $t$ RVs. The problem appears in many traditional statistical problems, e.g., check [3, 11–13];

(iii) The distribution of RVs $T_i$ can approximate other heavy-tailed symmetric distributions, e.g., $X/Y|Y\ge y_0$, where $X$ and $Y$ are respectively Gaussian and Rayleigh-distributed. In such scenarios, the distribution of their linear combination $Z$ may be extremely valuable. Interestingly, the sum of random variables in the form $X/Y|Y\ge y_0$ appears in the scenario proposed in [14], where the goal is to determine the number of active devices in a machine-type wireless communication network by relying on coordinated pilot transmissions without much signaling overhead, which facilitates the posterior data decoding procedures.

Unfortunately, the Student’s $t$ distribution is not closed under convolution [15], thus, deriving the exact distribution of $Z$ has been shown to be a cumbersome task, especially for an arbitrary number of degrees of freedom ${\nu }_i$ and number of addends $K$. For instance, the methods proposed in [2–5, 15] are restricted to the case of all $T_i$s having an odd number of degrees of freedom. Meanwhile, the PDF of $Z$ is given in [16] as an infinite series, but only for the specific case of $K=2$.

In general, approximation methods are often more tractable and appealing, which motivates our work in this paper. Specifically, we aim to accurately approximate the distribution of $Z$ in closed-form given that $T_i$ are independently distributed with ${\nu }_i>2$, $\forall i$, and no other constraints (the assumption that RVs $T_i$ are independent is common in the literature, e.g., [2, 4–6, 16]). To the best of our knowledge, this is the first work to (satisfactorily) address this.

1.3. Our Approach

For the approximation, we resort to a Student’s $t$ distribution fitting. Specifically, we aim to accurately fit $Z\sim {\sigma }_z\mathcal{T}{\nu }_z$ with ${\nu }_z>2$, which should hold, at least intuitively, as both distributions share the same symmetric and bell-shaped form. However, what might appear to be a simple and straightforward approach is not when considering that the Student’s t distributions have more than two degrees of freedom and no other constraints. We elaborate on this as follows:

The distribution fitting approaches commonly rely on moments (including $L-$ moments [17]) or cumulants matching. Specifically, at least two moments and/or cumulants of $Z$ are needed to match those of a scaled Student’s t distribution since such a distribution is characterized only by the scale ${\sigma }_z$, and the number of degrees of freedom ${\nu }_z$. However, the challenge lies in that ${\nu }_z$ (and each ${\nu }_i$) must be greater than the moment/cumulant order, while (i) the odd moments/cumulants cannot be used since they are zero, and (ii) the negative moments do not converge since $f_{Z}0>0$. This implies that (i) we cannot fit moments/cumulants of order higher than 2 in order to allow ${\nu }_z\in 2,\infty$, and (ii) we cannot rely on the first moment/cumulant. Meanwhile, fractional moments could be used, but they are complex and difficult to compute in general.

In this work, we resort to a fitting based on the second moment matching together with absolute moment matching, or CF matching, to circumvent the above issues. Note that using second moment matching is a natural choice given its simplicity, while exploiting the absolute moment also seems appealing. However, although absolute moments of any order, $\mathbb{E}{\cdot }^m$, with $m\in \mathbb{R},m<\nu$, could be used, they are cumbersome to derive due to the limited separability of the absolute value of a sum, thus, we focus on the simplest $m=1$ case. Finally, a CF matching performance is intriguing as it is not a commonly adopted approach in the literature for distribution fitting problems, especially because moments or other simpler statistics are often available, so we adopt it here given the special characteristics/challenges of the considered problem.

2. Computation of the Relevant Statistics of $Z$

2.1. Second Moment

Therefore, the second moment of $Z$ is given by$$\begin{array}{}\text{(8)}& \mathbb{E}{Z}^2={\displaystyle \sum }_{i=1}^K\mathbb{E}{T}_i^2={\displaystyle \sum }_{i=1}^K{\sigma }_i^2\frac{{\nu }_i}{{\nu }_i-2}\text{,}\end{array}$$which comes from leveraging the independence and zero-mean features of $T_i$ and from using (4) with $m=2$.

2.2. Characteristic Function

The CF of the sum of independent RVs matches the product of their independent CFs, thus$$\begin{array}{c}\text{(9)}& {\text{CF}}_{Z}r={\displaystyle \prod }_{i=1}^K{\text{CF}}_{{\sigma }_i{T}_i}r,={\displaystyle \prod }_{i=1}^K{\text{CF}}_{T_i}{\sigma }_{i}r,\\ ={\displaystyle \prod }_{i=1}^K\frac{{\sqrt{{\nu }_i}{\sigma }_{i}r}^{{\nu }_i/2}{K}_{{\nu }_i/2}\sqrt{{\nu }_i}{\sigma }_{i}r}{2^{{\nu }_i/2-1}\mathrm{\Gamma }{\nu }_i/2},\forall {\nu }_i>0,\\ =2^K{r/2}^{{\sum }_{i=1}^K{\nu }_i/2}{\displaystyle \prod }_{i=1}^K\frac{{\sqrt{{\nu }_i}{\sigma }_i}^{{\nu }_i/2}{K}_{{\nu }_i/2}\sqrt{{\nu }_i}{\sigma }_{i}r}{\mathrm{\Gamma }{\nu }_i/2},\forall {\nu }_i>0\text{.}\end{array}$$

2.3. Absolute Moment

The absolute moment of $Z$ obeys$$\begin{array}{}\text{(10)}& \mathbb{E}Z=\mathbb{E}{\displaystyle \sum }_{i=1}^K{\sigma }_i{T}_i\text{.}\end{array}$$

Remarkably, further simplifying (10) is not a trivial task. Furthermore, its computational complexity scales with $K$. Therefore, we focus on the case $K=2$, but leverage the corresponding results for the distribution fitting of the linear combination of any $K\ge 2$ RVs in Section 4.

The absolute moment for the case of $K=2$ can be computed as follows:$$\begin{array}{c}\text{(11)}& \mathbb{E}Z=\mathbb{E}{\sigma }_1{T}_1+{\sigma }_2{T}_2,\stackrel{a}{=}2\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\underset{-{\sigma }_2{x}_2/{\sigma }_1}{\overset{\infty }{{\displaystyle \int }}}{\sigma }_1{x}_1+{\sigma }_2{x}_2{f}_{T_1}{x}_1{f}_{T_2}{x}_2\mathrm{d}{x}_1\mathrm{d}{x}_2,\\ \stackrel{b}{=}2{\sigma }_1\underset{I_1}{\underbrace{\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\underset{{\sigma }_2{x}_2/{\sigma }_1}{\overset{\infty }{{\displaystyle \int }}}{x}_1{f}_{T_1}{x}_1{f}_{T_2}{x}_2\mathrm{d}{x}_1\mathrm{d}{x}_2}}+2{\sigma }_2\underset{I_2}{\underbrace{\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\underset{-{\sigma }_2{x}_2/{\sigma }_1}{\overset{\infty }{{\displaystyle \int }}}{x}_2{f}_{T_1}{x}_1{f}_{T_2}{x}_2\mathrm{d}{x}_1\mathrm{d}{x}_2}}\text{,}\end{array}$$where in (a), we exploit the fact that ${\sigma }_1{T}_1+{\sigma }_2{T}_2$ is symmetric around 0, thus, it can adopt positive and negative values with probability 0.5. Then, (b) comes after applying the integral operator to each of the integrand’s addends and leveraging the symmetry of ${\int }_{-a}^a{x}_1{f}_{T_1}{x}_1\mathrm{d}{x}_1=0$. Now, observe that $\mathbb{E}Z=2{\sigma }_1{I}_1+{\sigma }_2{I}_2$, and we are concerned with computing $I_1$ and $I_2$.

In the case of $I_1$, we have that$$\begin{array}{c}\text{(12)}& I_1={\int }_{-\infty }^{\infty }{\int }_{{\sigma }_2{x}_2/{\sigma }_1}^{\infty }{x}_1{f}_{T_1}{x}_1{f}_{T_2}{x}_2\mathrm{d}{x}_1\mathrm{d}{x}_2\stackrel{a}{=}\frac{2{\alpha }_{{\nu }_1}{\alpha }_{{\nu }_2}{\nu }_1}{{\nu }_1-1}{\int }_0^{\infty }{1+\frac{x_2^2{\sigma }_2^2}{{\nu }_1{\sigma }_1^2}}^{1-{\nu }_1/2}{1+\frac{x_2^2}{{\nu }_2}}^{-{\nu }_2+1/2}\mathrm{d}{x}_2\\ \stackrel{b}{=}\frac{{\alpha }_{{\nu }_1}{\alpha }_{{\nu }_2}B1/2,{\nu }_1+{\nu }_2-1/2}{{\nu }_1-1{\nu }_1^{-3/2}{\sigma }_2/{\sigma }_1}{}_2{}^{}{F}_1\frac{{\nu }_2+1}{2},\frac{1}{2},\frac{{\nu }_1+{\nu }_2}{2},1-\frac{{\nu }_1{\sigma }_1^2}{{\nu }_2{\sigma }_2^2}\text{,}\end{array}$$where (a) comes from solving the inner integral via [8], Equation (2.27.7)], while we leverage ([8], Equation (3.259.3)) to solve the remaining integral in (b).

In the case of $I_2$, we have that$$\begin{array}{c}\text{(13)}& I_2\stackrel{a}{=}{\int }_{-\infty }^{\infty }{x}_{2}1-F_{T_1}-\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2\stackrel{b}{=}-{\int }_{-\infty }^0{x}_2{F}_{T_1}-\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2-{\int }_0^{\infty }{x}_2{F}_{T_1}-\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2\\ \stackrel{c}{=}{\int }_0^{\infty }{x}_2{F}_{T_1}\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2-{\int }_0^{\infty }{x}_{2}1-F_{T_1}\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2\\ \stackrel{d}{=}2\underset{I_{\mathrm{2,2}}}{\underbrace{{\int }_0^{\infty }{x}_2{F}_{T_1}{\sigma }_2{x}_2/{\sigma }_1{f}_{T_2}{x}_2\mathrm{d}{x}_2}}-\underset{I_{\mathrm{2,1}}}{\underbrace{{\int }_0^{\infty }{x}_2{f}_{T_2}{x}_2\mathrm{d}{x}_2}}\text{,}\end{array}$$where (a) comes from using the CDF definition, and (b) from leveraging ${\int }_{-\infty }^{\infty }{x}_2{f}_{T_2}{x}_2\mathrm{d}{x}_2=0$, followed by splitting the integration region such that the sign of $x_2$ can be fixed accordingly. The latter, together with the symmetry of $T_1$, is exploited to attain (c), while (d) is immediately obtained after simple algebraic transformations. Note that $I_2=2I_{\mathrm{2,2}}-I_{\mathrm{2,1}}$, where $I_1$ and $I_2$ require integral computations as shown in (13). Fortunately, their calculation can be further simplified as described next.

In the case of $I_{\mathrm{2,1}}$, we have that$$\begin{array}{}\text{(14)}& I_{\mathrm{2,1}}=\frac{\sqrt{{\nu }_2}\mathrm{\Gamma }{\nu }_2-1/2}{2\mathrm{\Gamma }{\nu }_2/2\sqrt{\pi }}\text{,}\end{array}$$which comes from using ([8], Equation (3.241.4)). Meanwhile$$\begin{array}{c}\text{(15)}& I_{\mathrm{2,2}}={\int }_0^{\infty }{x}_2{F}_{T_1}\frac{{\sigma }_2{x}_2}{{\sigma }_1}{f}_{T_2}{x}_2\mathrm{d}{x}_2,\stackrel{a}{=}{\int }_0^{\infty }{x}_{2}1-\frac{1}{2}{\mathcal{I}}_{\frac{{\nu }_1}{{\sigma }_2^2{x}_2^2/{\sigma }_1^2+{\nu }_1}}\frac{{\nu }_1}{2},\frac{1}{2}{f}_{T_2}{x}_2\mathrm{d}{x}_2,\\ \stackrel{b}{=}{\int }_0^{\infty }{x}_2{f}_{T_2}{x}_2\mathrm{d}{x}_2-\frac{{\alpha }_{{\nu }_2}}{2}\underset{I_{\mathrm{2,2}}^{\prime }}{\underbrace{{\int }_0^{\infty }{x}_2{\mathcal{I}}_{{\nu }_1/{\sigma }_2^2{x}_2^2/{\sigma }_1^2+{\nu }_1}\frac{{\nu }_1}{2},\frac{1}{2}{1+\frac{x_2^2}{{\nu }_2}}^{-{\nu }_2+1/2}\mathrm{d}{x}_2}},\\ \stackrel{c}{=}{I}_{\mathrm{2,1}}-\frac{{\alpha }_{{\nu }_2}}{2}{I}_{\mathrm{2,2}}^{\prime }\text{,}\end{array}$$where (a) comes from substituting (2), (b) follows from rearranging terms, while (c) is obtained by leveraging (14).

Now, to simplify $I_{\mathrm{2,2}}^{\prime }$, we introduce the following variable transformation$$\begin{array}{}\text{(16)}& z^2\triangleq \frac{{\nu }_1}{{\sigma }_2^2{x}_2^2/{\sigma }_1^2+{\nu }_1}⟶x_2^2=\frac{{\nu }_1{\sigma }_1^2}{{\sigma }_2^2}\frac{1}{z^2}-1\text{,}\end{array}$$thus, $x_2\mathrm{d}{x}_2=-{\nu }_1{\sigma }_1^2/{\sigma }_2^2{z}^3\mathrm{d}z$. Then$$\begin{array}{}\text{(17)}& I_{\mathrm{2,2}}^{\prime }={\omega }_1{\int }_0^1{z}^{-3}{\mathcal{I}}_{z^2}\frac{{\nu }_1}{2},\frac{1}{2}{{\omega }_2+\frac{1}{z^2}}^{-{\nu }_2+1/2}\mathrm{d}z\text{,}\end{array}$$where ${\omega }_1\triangleq {{\nu }_1{\sigma }_1^2/{\sigma }_2^2}^{-{\nu }_2-1/2}{\nu }_2^{{\nu }_2+1/2}$, ${\omega }_2\triangleq {\nu }_2{\sigma }_2^2/{\nu }_1{\sigma }_1^2-1$. Now, we leverage the Taylor series expansion of an incomplete regularized beta function [18] to write$$\begin{array}{}\text{(18)}& I_{z^2}\frac{{\nu }_1}{2},\frac{1}{2}=\frac{2}{\sqrt{\pi }B{\nu }_1/2,1/2}{\displaystyle \sum _{i=0}^{\infty }\frac{z^{{\nu }_1+2i}\Gamma 1/2+i}{{\nu }_1+2ii!}}.\end{array}$$

Then, by substituting (18) into (17), one gets$$\begin{array}{c}\text{(19)}& I_{\mathrm{2,2}}^{\prime }=\frac{2{\omega }_1/\sqrt{\pi }}{B\frac{{\nu }_1}{2},\frac{1}{2}}{\displaystyle \sum }_{i=0}^{\infty }\frac{\mathrm{\Gamma }\frac{1}{2}+i}{{\nu }_1+2ii!}{\int }_0^1\frac{{{\omega }_2+\frac{1}{z^2}}^{-\frac{{\nu }_2+1}{2}}}{z^{3-{\nu }_1-2i}}\mathrm{d}z\stackrel{a}{=}\frac{{\omega }_1/\sqrt{\pi }}{B\frac{{\nu }_1}{2},\frac{1}{2}}{\displaystyle \sum }_{i=0}^{\infty }\frac{\mathrm{\Gamma }\frac{1}{2}+i}{{\nu }_1+2ii!}{\int }_0^1\frac{{y{\omega }_2+1}^{-\frac{{\nu }_2+1}{2}}}{y^{\frac{3-{\nu }_1-{\nu }_2-2i}{2}}}\mathrm{d}y\\ \stackrel{b}{=}\frac{2{\omega }_1}{\sqrt{\pi }B\frac{{\nu }_1}{2},\frac{1}{2}}{\displaystyle \sum }_{i=0}^{\infty }\frac{\mathrm{\Gamma }i+\frac{1}{2}{}_2{F}_1\frac{{\nu }_2+1}{2},\frac{{\nu }_1+{\nu }_2+2i-1}{2},\frac{{\nu }_1+{\nu }_2+2i+1}{2},-{\omega }_2}{{\nu }_1+2i{\nu }_1+{\nu }_2+2i-1i!}\text{,}\end{array}$$where (a) comes from exploiting the variable transformation $y=z^2$, while the integral is solved in (b) by leveraging ([8], Equation (3.197.8)). Then $I_{\mathrm{2,2}}^{\prime }$ can be easily estimated by truncating the infinite sum in (19). As illustrated in Figure 1, the relative approximation error decreases following a power law decay, thus, a relatively small number of addends is needed, especially for small ${\nu }_2$, ${\nu }_1$, and large ${\omega }_2$.

[figure(s) omitted; refer to PDF]

By combining (12), (13), (14), (15), and (19), and substituting them into (11), we obtain$$\begin{array}{c}\text{(20)}& \mathbb{E}Z=\frac{B{1/2,{\nu }_1+{\nu }_2-1/2}_2{F}_1{\nu }_2+1/2,1/2,{\nu }_1+{\nu }_2/2,1-{\nu }_1{\sigma }_1^2/{\nu }_2{\sigma }_2^2}{{\nu }_1-1{\nu }_1^{-3/2}{\sigma }_2/2{\alpha }_{{\nu }_1}{\alpha }_{{\nu }_2}{\sigma }_1^2}+\frac{{\sigma }_2\mathrm{\Gamma }{\nu }_2-1/2}{\mathrm{\Gamma }{\nu }_2/2\sqrt{\pi /{\nu }_2}}-\frac{4{\alpha }_{{\nu }_2}{\sigma }_2{\omega }_1}{\sqrt{\pi }B{\nu }_1/2,1/2}\times {\displaystyle \sum }_{i=0}^{\infty }\frac{\mathrm{\Gamma }{i+1/2}_2{F}_1{\nu }_2+1/2,{\nu }_1+{\nu }_2+2i-1/2,{\nu }_1+{\nu }_2+2i+1/2,-{\omega }_2}{{\nu }_1+2i{\nu }_1+{\nu }_2+2i-1i!}\text{.}\end{array}$$

2.4. Special Case: Sum of i.i.d. Student’s $t$ RVs

Herein, we consider the special case of the sum of i.i.d. Student’s $t$ RVs, i.e., ${\sigma }_i=\sigma$, ${\nu }_i=\nu$, and $\forall i=1,\dots ,K$. With this in hand, the computation of $\mathbb{E}{Z}^2$, ${\text{CF}}_{Z}r$, and $\mathbb{E}Z$ (for $K=2$) can be more easily obtained as follows:

Under such conditions, (8) and (9) directly simplify to$$\begin{array}{c}\text{(21)}& \mathbb{E}{Z}^2=\frac{\sigma \nu K}{\nu -2}\text{,}{\text{CF}}_{Z}r={\text{CF}}_T^K\sigma r={\frac{{\sqrt{\nu }\sigma r}^{\nu /2}{K}_{\nu /2}\sqrt{\nu }\sigma r}{2^{\nu /2-1}\mathrm{\Gamma }\nu /2}}^K\text{.}\end{array}$$

As for $\mathbb{E}Z$ for $K=2$, for which $\mathbb{E}Z=2\sigma I_1+I_{\mathrm{2,1}}-{\alpha }_{\nu }{I}_{\mathrm{2,2}}^{\prime }$, one may depart from (a) in (12) to write$$\begin{array}{c}\text{(22)}& I_1=\frac{2{\alpha }_{\nu }^2\nu }{\nu -1}{\int }_0^{\infty }{1+x_2^2/\nu }^{-\nu }\mathrm{d}{x}_2=\frac{\mathrm{\Gamma }\nu -1/2\mathrm{\Gamma }\nu -1/2}{2^{\nu }\mathrm{\Gamma }{\nu /2}^3{\nu }^{-1/2}}\text{,}\end{array}$$which leverages ([8], Equation (3.251.2)) and (3) followed by some simple algebraic simplifications. Meanwhile, (14) can be directly used with $\nu$ instead of ${\nu }_2$. Now, (19) can be computed by departing from (17) as$$\begin{array}{c}\text{(23)}& I_{\mathrm{2,2}}^{\prime }=\nu {\int }_0^1{z}^{\nu -2}{\mathcal{I}}_{z^2}\frac{\nu }{2},\frac{1}{2}\mathrm{d}z\stackrel{a}{=}\frac{\mathrm{\Gamma }\nu +1/2}{\mathrm{\Gamma }\nu /2\sqrt{\pi }}{\int }_0^1{x}^{\nu -3/2}{}_2{F}_1{\nu }_2,\frac{1}{2},1+\frac{\nu }{2},x\mathrm{d}x\\ \stackrel{b}{=}\frac{\nu }{\nu -1}1-\frac{\sqrt{\pi }\mathrm{\Gamma }\nu -1/2}{2^{\nu -1}\mathrm{\Gamma }{\nu /2}^2}\text{,}\end{array}$$where (a) comes from stating the regularized incomplete beta function in terms of a hypergeometric function according to ([7], Equation (8.17.7)) and using $x\triangleq z^2⟶\mathrm{d}x=2z\mathrm{d}z$, while (b) follows from leveraging [19] together with iterative integration by parts. Finally, by combining the above results, we obtain$$\begin{array}{}\text{(24)}& \mathbb{E}Z=\frac{\sigma \sqrt{\nu }\mathrm{\Gamma }\nu -1/2\mathrm{\Gamma }\nu -1/2}{2^{\nu -2}\mathrm{\Gamma }{\nu /2}^3}\text{.}\end{array}$$

3. Distribution Fitting

Next, we follow two different distribution-fitting approaches. The first approach is based on matching the second and absolute moments, while the second approach relies on matching the second moment and the CF for a certain $r$. We also illustrate their accuracy.

3.1. Fitting Based on Second and Absolute Moments

According to (4) with $m=2$ and (6) with $m=1$, the set of equations to solve is as follows:$$\begin{array}{}\text{(25)}& \frac{{\sigma }_z^2{\nu }_z}{{\nu }_z-2}=\mathbb{E}{Z}^2,\frac{{\sigma }_z\sqrt{{\nu }_z}\mathrm{\Gamma }{\nu }_z-1/2}{\mathrm{\Gamma }{\nu }_z/2\sqrt{\pi }}=\mathbb{E}Z\text{,}\end{array}$$with variables ${\sigma }_z,{\nu }_z$. Recall that $\mathbb{E}{Z}^2$ is given by (8), while $\mathbb{E}Z$ is given in (20) for the case of $K=2$. By isolating ${\sigma }_z$ in the first equation, i.e., ${\sigma }_z=\sqrt{{\nu }_z-2\mathbb{E}{Z}^2/{\nu }_z}$, and substituting it into the second equation, the system of (25) transforms to$$\begin{array}{}\text{(26)}& h{\nu }_z\triangleq \frac{\mathrm{\Gamma }{\nu }_z-1/2\sqrt{{\nu }_z-2}}{\mathrm{\Gamma }{\nu }_z/2}=\frac{\sqrt{\pi }\mathbb{E}Z}{\sqrt{\mathbb{E}{Z}^2}}\text{.}\end{array}$$

We observe that attaining an exact closed-form solution for ${\nu }_z$ in (26) is not viable, thus, we resort to a low-complex approximation of $h{\nu }_z$. For this, we plot $h{\nu }_z$ vs. ${\nu }_z$ in Figure 2, and realize that $h{\nu }_z$ has approximately the form of a quotient of two linear functions. Hence, we state$$\begin{array}{}\text{(27)}& h{\nu }_z\approx \frac{p_1{\nu }_z+p_2}{{\nu }_z+p_3}\text{.}\end{array}$$

[figure(s) omitted; refer to PDF]

Here, we know that $h2=0$, and$$\begin{array}{c}\text{(28)}& \underset{{\nu }_z⟶\infty }{\lim }h{\nu }_z=\underset{{\nu }_z⟶\infty }{\lim }\frac{\mathrm{\Gamma }{\nu }_z-1/2\sqrt{{\nu }_z-2}}{\mathrm{\Gamma }{\nu }_z/2},\stackrel{a}{=}\underset{{\nu }_z⟶\infty }{\lim }\sqrt{\frac{2{\nu }_z-2}{{\nu }_z}},\\ =\sqrt{2}\text{,}\end{array}$$where (a) comes from using ([7], Equation 5.11.12). Then, in order to satisfy such conditions, we can directly set $p_1$ and $p_2$ as follows:$$\begin{array}{}\text{(29)}& \underset{{\nu }_z⟶\infty }{\lim }h{\nu }_z=p_1=\sqrt{2}\text{,}\text{(30)}& h2=\frac{2p_1+p_2}{2+p_3}=0⟶2p_1+p_2=0⟶p_2=-2\sqrt{2}\text{,}\end{array}$$where (30) uses the result in (29) in the last step. Finally, $p_3=-\sqrt{3}$ can be obtained easily by standard curve fitting. The accuracy of (27) is also depicted in Figure 2.

With (27) in place, one can estimate ${\nu }_z$ and ${\sigma }_z$ as follows:$$\begin{array}{}\text{(31)}& {\nu }_z^{\star }\approx \frac{\pi \mathbb{E}Z-2\sqrt{2\mathbb{E}{Z}^2}}{\sqrt{\pi }\mathbb{E}Z-\sqrt{2\mathbb{E}{Z}^2}}\text{,}\text{(32)}& {\sigma }_z^{\star }\approx \sqrt{\frac{\pi -2\sqrt{\pi }\mathbb{E}Z\mathbb{E}{Z}^2}{\pi \mathbb{E}Z-2\sqrt{2\mathbb{E}{Z}^2}}}\text{,}\end{array}$$and use the approximation $Z\sim {\sigma }_z^{\star }\mathcal{T}{\nu }_z^{\star }$. Such a distribution fit is illustrated in Figure 3 and evinces the appropriateness of our approach. All in all, the fitting accuracy only seems to be affected by the distribution tails, which is expected considering that only two (low-order) features of the child distributions are used for the fitting. Nevertheless, even in the tails region, the accuracy is surprisingly good, being only critically affected when the child distributions have significantly diverging degrees of freedom and scaling factors in opposite directions, e.g., small ${\nu }_1,{\sigma }_2$ and large ${\nu }_2,{\sigma }_1$.

[figure(s) omitted; refer to PDF]

3.1.1. Linear Combination of $K>2$ Student’s $t$ RVs

Notice that if a scaled Student’s $t$ distribution fits accurately the distribution of the linear combination of $K=2$ Student’s $t$ RVs, then such a Student’s $t$ fitting approach applies for the linear combination of any $K\ge 2$ Student’s $t$ RVs. This can be easily shown by induction as follows:

According to the previous subsection’s results, we can state that$$\begin{array}{}\text{(33)}& {\sigma }_1{T}_1+{\sigma }_2{T}_2\sim {\sigma }_{z_2}\mathcal{T}{\nu }_{z_2}\end{array}$$holds approximately. Here, ${\sigma }_{z_n}=g_1{\sigma }_1,\dots ,{\sigma }_n,{\nu }_1,\dots ,{\nu }_n$ and ${\nu }_{z_n}=g_2{\sigma }_1,\dots ,{\sigma }_n,{\nu }_1,\dots ,{\nu }_n$, where $g_1\cdot$ and $g_2\cdot$ are transformation functions (given by (31) and (32) in the particular case of $n=2$). Then, assume that ${\sum }_{i=1}^K{\sigma }_i{T}_i\sim {\sigma }_{z_K}\mathcal{T}{\nu }_{z_K}$, and observe that$$\begin{array}{c}\text{(34)}& {\displaystyle \sum }_{i=1}^{K+1}{\sigma }_i{T}_i={\sigma }_{K+1}\underset{\mathcal{T}{\nu }_{K+1}}{\underbrace{T_{K+1}}}+\underset{{\sigma }_{z_K}\mathcal{T}{\nu }_{z_K}}{\underbrace{{\displaystyle \sum }_{i=1}^K{\sigma }_i{T}_i}}\sim \underset{{\sigma }_{z_{K+1}}}{\underbrace{g_1{\sigma }_{z_K},{\sigma }_{K+1},{\nu }_{z_K},{\nu }_{K+1}}}\mathcal{T}\underset{{\nu }_{z_{K+1}}}{\underbrace{g_1{\sigma }_{z_K},{\sigma }_{K+1},{\nu }_{z_K},{\nu }_{K+1}}}\text{.}\end{array}$$

Thus, proving the hypothesis.

Wrapping up, the linear combination of $K\ge 2$ Student’s $t$ RVs is approximately distributed as ${\sigma }_{z_K}\mathcal{T}{\nu }_{z_K}$, where ${\sigma }_{z_K}$, ${\nu }_{z_K}$ can be iteratively obtained as follows:$$\begin{array}{c}\text{(35)}& {\sigma }_{z_{i+1}}=g_1{\sigma }_{z_i},{\sigma }_{i+1},{\nu }_{z_i},{\nu }_{i+1}\text{,}{\nu }_{z_{i+1}}=g_2{\sigma }_{z_i},{\sigma }_{i+1},{\nu }_{z_i},{\nu }_{i+1}\text{,}\end{array}$$$\forall i\ge 1$, where ${\sigma }_{z_1}={\sigma }_1$, ${\nu }_{z_1}={\nu }_1$, and $g_1{\sigma }_1,{\sigma }_2,{\nu }_1,{\nu }_2$, $g_2{\sigma }_1,{\sigma }_2,{\nu }_1,{\nu }_2$ are respectively given by (31) and (32), which can be computed by leveraging (8) and (20). The accuracy of such a procedure has been corroborated by several simulation campaigns, and it is illustrated here in Figure 4 for an example set of distribution parameters.

[figure(s) omitted; refer to PDF]

3.2. Fitting Based on Second Moment and Characteristic Function

According to (4) with $m=2$ and (5), the set of equations to solve is as follows:$$\begin{array}{}\text{(36)}& \frac{{\sigma }_z^2{\nu }_z}{{\nu }_z-2}=\mathbb{E}{Z}^2,\frac{{\sqrt{{\nu }_z}{\sigma }_{z}r}^{{\nu }_z/2}{K}_{{\nu }_z/2}\sqrt{{\nu }_z}{\sigma }_{z}r}{2^{{\nu }_z/2-1}\mathrm{\Gamma }{\nu }_z/2}={\text{CF}}_{Z}r,\end{array}$$with variables ${\sigma }_z,{\nu }_z$. Recall that $\mathbb{E}{Z}^2$ is given by (8), while ${\text{CF}}_{Z}r$ is given in (9). As in Section 3.1, we first isolate ${\sigma }_z$ in the first equation, i.e., ${\sigma }_z=\sqrt{{\nu }_z-2\mathbb{E}{Z}^2/{\nu }_z}$, and then substitute it into the second equation. By doing this, the system of (36) transforms to$$\begin{array}{}\text{(37)}& g{\nu }_z\triangleq \frac{{r\sqrt{{\nu }_z-2\mathbb{E}{Z}^2}}^{{\nu }_z/2}{K}_{{\nu }_z/2}r\sqrt{{\nu }_z-2\mathbb{E}{Z}^2}}{2^{{\nu }_z/2-1}\mathrm{\Gamma }{\nu }_z/2}={\text{CF}}_{Z}r\text{.}\end{array}$$

It can be shown that $g{\nu }_z$ is a decreasing function and that ${\lim }_{{\nu }_z⟶2^{+}}g{\nu }_z=1$ and ${\lim }_{{\nu }_z⟶\infty }g{\nu }_z=lr\sqrt{\mathbb{E}{Z}^2}\ge 0$ for some decreasing function $l\cdot$. This is illustrated in Figure 5 and implies that there is at most one real solution ${\nu }_z^{\star }$ for (37), which can be easily found via the bisection method. After this, one sets ${\sigma }_z^{\star }=\sqrt{{\nu }_z^{\star }-2\mathbb{E}{Z}^2/{\nu }_z^{\star }}$ and uses the approximation $Z\sim {\sigma }_z^{\star }\mathcal{T}{\nu }_z^{\star }$.

[figure(s) omitted; refer to PDF]

The main challenge with this approach is the proper setting of $r$. On one hand, (37) may not have solution for certain values of $r$. On the other hand, different feasible values of $r$ may lead to significantly different fitting accuracy figures. We investigate these issues in Section 4.3. For now, let us conveniently set$$\begin{array}{}\text{(38)}& r=\mathbb{E}{Z^2}^{-1/2},\end{array}$$so that $g{\nu }_z$ simplifies to$$\begin{array}{}\text{(39)}& g{\nu }_z\triangleq \frac{{{\nu }_z-2}^{{\nu }_z/4}{K}_{{\nu }_z/2}\sqrt{{\nu }_z-2}}{2^{{\nu }_z/2-1}\mathrm{\Gamma }{\nu }_z/2}\text{.}\end{array}$$

Observe that attaining an exact closed-form solution for ${\nu }_z$ in $g{\nu }_z={\text{CF}}_Z\mathbb{E}{Z^2}^{-1/2}$ using (39) is still not viable. Fortunately, $g{\nu }_z$ can be accurately approximated by a very tractable fractional function of the form$$\begin{array}{}\text{(40)}& g{\nu }_z\approx \frac{p_1{\nu }_z+p_2}{{\nu }_z+p_3}\text{.}\end{array}$$

Meanwhile, we know that ${\lim }_{{\nu }_z⟶2^{+}}g{\nu }_z=1$, and ${\lim }_{{\nu }_z⟶\infty }g{\nu }_z\approx 0.6070$ as shown in Figure 5. Then, in order to satisfy such conditions, we can directly set $p_1$ and $p_2$ as follows:$$\begin{array}{}\text{(41)}& \underset{{\nu }_z⟶\infty }{\lim }g{\nu }_z=p_1=0.607\text{,}\text{(42)}& g2\approx \frac{2p_1+p_2}{2+p_3}=1⟶2p_1+p_2=2+p_3⟶p_3=2p_1+p_2-2=p_2-0.786\text{,}\end{array}$$where (42) leverages the result in (41) in the last step. Finally, $p_2=-0.7606$ can be obtained easily by standard curve fitting, and then one can set $p_3=-0.7606-0.786=-1.5466$. The accuracy of (40) is depicted in Figure 6.

[figure(s) omitted; refer to PDF]

With (40) in place, one can estimate ${\nu }_z$ and ${\sigma }_z$ as$$\begin{array}{c}\text{(43)}& {\nu }_z^{\star }\approx \frac{0.7606-1.5466\times {\text{CF}}_Z\mathbb{E}{Z^2}^{-1/2}}{0.607-{\text{CF}}_Z\mathbb{E}{Z^2}^{-1/2}}\text{,}{\sigma }_z^{\star }\approx \sqrt{\frac{{\text{CF}}_Z\mathbb{E}{Z^2}^{-1/2}-1\mathbb{E}{Z}^2}{1.6775-3.4111\times {\text{CF}}_Z\mathbb{E}{Z^2}^{-1/2}}}\text{,}\end{array}$$and use the approximation $Z\sim {\sigma }_z^{\star }\mathcal{T}{\nu }_z^{\star }$. Such a distribution fitting is illustrated in Figures 7 and 8, and evinces the appropriateness of our approach. Finally, notice that the results here agree also with our previous observations around Figure 3.

[figure(s) omitted; refer to PDF]

4. Fitting Accuracy Analysis

In this section, we assess the accuracy of the fitting methods discussed in Section 3 by adopting the Bhattacharyya distance metric [20]. This metric measures the similarity of two probability distributions, which in this case are the true distribution of $Z$, $f_{Z}z$, and the approximate scaled Student’s $t$ distribution $f_{\widehat{Z}}z$ based on one of the proposed fitting approaches. Notice that since the true/exact distribution of $Z$, $f_{Z}z$, is unknown and difficult to compute for a general $K$, we leverage a Monte Carlo approach to estimate the Bhattacharyya distance. Specifically, such a metric is computed as follows:$$\begin{array}{}\text{(44)}& d_{B}Z,\widehat{Z}=-\ln {\displaystyle \sum }_{n=1}^N\sqrt{f_Z{z}_n{f}_{\widehat{Z}}{z}_n}\text{,}\end{array}$$where $z_n$ is the $n-$ th sample taken from $\sigma {\sum }_{i=1}^K{T}_i$, and $f_Z{z}_n$ is estimated using its histogram. As $N⟶\infty$, (44) approaches the exact Bhattacharyya distance between the continuous probability distributions of $Z$ and $\widehat{Z}$. Finally, we focus on the special case of the sum of i.i.d. RVs for simplicity, adopt $N={10}^6$, and set $\sigma =1$ without loss of generality.

4.1. Absolute Moment vs. CF-Based Fitting

Figure 9 shows the fitting accuracy of the two methods proposed in this work together with that of a benchmark approach based on the second and fourth moments matching. As commented earlier, such an approach can only be used when $\nu \ge 4$, which was one of the key motivations for our work. From the figure, we can tell that (i) the proposed fitting approaches are more accurate than the benchmark based on second and fourth moments matching, which is also restricted to the cases where $\nu \ge 4$; (ii) the fitting based on the CF matching is generally more accurate than the one based on absolute moment matching, although both approaches tend to converge as $\nu$ increases; and (iii) for relatively small $\nu$, the proposed fitting approaches are more accurate for a smaller $K$, while this behavior might be reverted as $\nu$ increases.

[figure(s) omitted; refer to PDF]

4.2. Scaling Laws of the Fitting Parameters

Herein, we leverage the i.i.d. Assumption to illustrate in Figure 10 how the parameters of the fitting distribution scale with $K$. For this, we adopt only the CF-based fitting, as our previous results indicate that it is the most accurate. We observe that both parameters of the fitting distribution, the scale ${\sigma }_z$ and the number of degrees of freedom ${\nu }_z$, increase with $K$. Specifically, ${\sigma }_z,{\nu }_z\sim {\gamma }_1{K}^{{\gamma }_2}+{\gamma }_3$ with ${\gamma }_1,{\gamma }_2>0$ and ${\gamma }_2<1$ (concave increase). Moreover, ${\sigma }_z$ and ${\nu }_z$ are respectively increasing and decreasing functions of $\nu$.

[figure(s) omitted; refer to PDF]

4.3. Performance Impact of the CF Parameter Setting

As commented earlier in Section 3.2, different choices of $r$ in (36) (or directly in (37)) may lead to significantly different fitting accuracy figures. Here, we investigate these issues by illustrating the estimated Bhattacharyya distance as a function of $r$ in Figure 11. Observe that for relatively large values of $\nu$, the accuracy depends little on the specific value of $r$. This situation changes drastically as $\nu$ decreases approaches two, under which proper setting $r$ becomes more critical. Indeed, there is an accuracy-optimum value of $r$, especially noticeable when $\nu$ is small. Notice that although the optimum $r$ may be cumbersome to determine beforehand in practice. Nevertheless, and as a rule of dumb, relatively small values of $r$ are usually preferred, and our proposal in (38) appears to be a valid (and simple) configuration approach.

[figure(s) omitted; refer to PDF]

5. Conclusion

In this work, we proposed two fitting approaches for the distribution of linear combinations of the Student’s $t$ RVs with more than two degrees of freedom. They leverage the second moment together with either the first absolute moment or the CF to fit the distribution to that of a scaled Student’s $t$ RV. For the former, we first analytically obtained the absolute moment of a linear combination of $K=2$ Student’s $t$ RVs and then generalized to $K\ge 2$ through a simple iterative procedure, while the fitting is direct for the latter but its accuracy depends on the CF parameter. Notably, we proposed a simple CF parameter configuration and showed that it can lead to high fitting accuracy. We resorted to Monte Carlo simulations and adopted the Bhattacharyya distance metric for numerically quantifying the fitting accuracy. We showed that the CF-based fitting can usually outperform the absolute moment -based fitting, although the accuracy provided by both approaches converges when the $t-$ RVs have a sufficiently large number of degrees of freedom. Interestingly, both proposed approaches outperform benchmark fitting based on second and fourth moments matching, which is only applicable when all $t-$ RVs have at least four degrees of freedom. Interestingly, both the scale and number of degrees of freedom of the fitting distribution were shown to increase almost linearly with $K$.

Finally, notice that our work opens the path to even more general fitting approaches. Indeed, by discarding the second moment and leveraging only the absolute moments and CFs, one may be able to accurately characterize the distribution of the linear combinations of $t-$ RVs with more than one degree of freedom. Moreover, by solely relying on CFs matching with different parameters $r$, one may completely remove the constraint on the number of degrees of freedom and arbitrarily fit the distribution of any linear combination of $t-$ RVs to that of a scaled Student’s t RV. However, further studies are required on how to achieve this in an optimal and/or simple manner.

Acknowledgments

This work was supported by the Academy of Finland (6G Flagship Program under Grant 346208) and the Finnish Foundation for Technology Promotion.

References

[1] V. E. Bening, V. Y. Korolev, "On an application of the Student distribution in the theory of probability and mathematical statistics," Theory of Probability and Its Applications, vol. 49 no. 3, pp. 377-391, 2005.

[2] G. A. Walker, J. G. Saw, "The distribution of linear combinations of t-variables," Journal of the American Statistical Association, vol. 73 no. 364, pp. 876-878, DOI: 10.1080/01621459.1978.10480117, 1978.

[3] V. Witkovský, "Exact distribution of positive linear combinations of inverted chi-square random variables with odd degrees of freedom," Statistics & Probability Letters, vol. 56 no. 1, pp. 45-50, DOI: 10.1016/s0167-7152(01)00165-1, 2002.

[4] C. Berg, C. Vignat, "Linearization coefficients of bessel polynomials and properties of student t-distributions," Constructive Approximation, vol. 27 no. 1, pp. 15-32, DOI: 10.1007/s00365-006-0643-6, 2008.

[5] S. Nadarajah, D. Dey, "Convolutions of the t distribution," Computers & Mathematics with Applications, vol. 49 no. 5-6, pp. 715-721, 2005.

[6] M. Ahsanullah, B. G. Kibria, M. Shakil, Normal and Student’s T Distributions and Their Applications, 2014.

[7] F. W. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, NIST Handbook of Mathematical Functions Hardback and CD-ROM, 2010.

[8] I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series, and Products, 2014.

[9] W. R. Fairweather, "A method of obtaining an exact confidence interval for the common mean of several normal populations," Applied Statistics, vol. 21 no. 3, pp. 229-233, DOI: 10.2307/2346273, 1972.

[10] K. Krishnamoorthy, Y. Lu, "Inferences on the common mean of several normal populations based on the generalized variable method," Biometrics, vol. 59 no. 2, pp. 237-247, DOI: 10.1111/1541-0420.00030, 2003.

[11] S.-H. Kim, A. S. Cohen, "On the behrens-fisher problem: a review," Journal of Educational and Behavioral Statistics, vol. 23 no. 4, pp. 356-377, DOI: 10.3102/10769986023004356, 1998.

[12] S. Weerahandi, Exact Statistical Methods for Data Analysis, 2003.

[13] M. Hollander, D. A. Wolfe, E. Chicken, Nonparametric Statistical Methods, 2013.

[14] O. L. López, G. Brante, R. D. Souza, M. Juntti, M. Latva-aho, "Coordinated pilot transmissions for detecting the signal sparsity level in a massive iot network under rayleigh fading," 2022. https://arxiv.org/abs/2205.00406

[15] N. C. Petroni, "Mixtures in nonstable lévy processes," Journal of Physics A: Mathematical and Theoretical, vol. 40 no. 10, pp. 2227-2250, DOI: 10.1088/1751-8113/40/10/001, 2007.

[16] C. Berg, C. Vignat, "On the density of the sum of two independent Student t-random vectors," Statistics & Probability Letters, vol. 80 no. 13-14, pp. 1043-1055, DOI: 10.1016/j.spl.2010.02.019, 2010.

[17] T. C. Headrick, "A characterization of power method transformations through L-moments," Journal of Probability and Statistics, vol. 2011, pp. 2011-2022, DOI: 10.1155/2011/497463, 2011.

[18] "Regularized incomplete beta function: series representations," 1998. https://functions.wolfram.com/GammaBetaErf/BetaRegularized/06/01/03/01/01/

[19] "Gauss hypergeometric function 2F1: Integration," 1998. https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/21/ShowAll.html

[20] T. Kailath, "The divergence and bhattacharyya distance measures in signal selection," IEEE Transactions on Communications, vol. 15 no. 1, pp. 52-60, DOI: 10.1109/tcom.1967.1089532, 1967.