# The distribution of the difference between two correlated non-central t distribution

Suppose a binormal population $$\{X, Y\}$$ with means $$\mathbf{\mu} = \{\mu_1,\mu_2\} \ne \{0,0\}$$ and covariance $$\Sigma= \sigma^2\begin{bmatrix}1 & \rho\\ \rho &1 \end{bmatrix}$$. Let $$S^2$$ be an estimator of $$\sigma^2$$ with $$f$$ degrees of freedom.

It is known that the variates $$\{X/S, Y/S\}$$ follows a bivariate non-central $$t$$ distribution (e.g., Kshirsagar, 1961).

Walgren (1980) derived the distribution of the product $$X/S \times Y/S$$. Is there a derivation for their difference, $$X/S - Y/S$$?

Edit: In the present context, I use the pooled standard deviation, that is the mean of the separate standard deviations,

$$S_X = \sum_{i=1}^n (X_i-\bar{X})^2 / (n-1)$$ $$S_Y = \sum_{i=1}^n (Y_i-\bar{Y})^2 / (n-1)$$

with $$S = \sqrt{ (S_X^2 + S_Y^2)/2 }$$. This estimate of $$\sigma$$ is independent of both $$X$$ and $$Y$$, and as shown here is a chi-square distribution with degrees of freedom $$2(n-1)/(1+\rho^2)$$.

• Please state which estimator of $\sigma^2$ you have in mind. (At least three would be natural in this setting and there are many others one could use as well.)
– whuber
Mar 10, 2022 at 23:25
• I edited the question. Mar 11, 2022 at 6:51
• I don't have a "derivation" but I'm getting $\frac{\Gamma \left(\frac{v+1}{2}\right) \left(2 (1-\rho ) \sigma ^2 v\right)^{v/2} \left(2 (1-\rho ) \sigma ^2 v+(-\mu_1 +\mu_2+z)^2\right)^{-\frac{1}{2} (v+1)}}{\sqrt{\pi } \Gamma \left(\frac{v}{2}\right)}$ for the pdf using Mathematica with $Z=X/S-Y/S$.
– JimB
Mar 11, 2022 at 7:03
• Do you mean that $(n-1)S^2/\sigma^2$ has a chi-square distribution with degrees of freedom $2(n-1)/(1+\rho^2)$ ?
– JimB
Mar 11, 2022 at 19:40
• @JimB The result from Ben in the link provided indicates that $2(n-1)/(1+\rho^2) S^2 / \sigma^2 \sim (1+\rho^2)/(2(n-1)) \chi^2_v$ in which $v = 2(n-1)/(1+\rho^2)$. Mar 11, 2022 at 20:01

Here is an approach that shows how to obtain a numerical approximation to the probability density function of $$Z=X/S-Y/S$$. (I haven't been successful in finding an analytic solution.)

Using the joint density of $$X/S$$ and $$Y/S$$ found in Kshirsagar 1961 (as given in the question):

r = {{1, \[Rho]}, {\[Rho], 1}}; (* Correlation matrix *)
\[Mu] = {\[Mu]x, \[Mu]y};
t = {x, y};
f = 2 (n - 1)/(1 + \[Rho]^2); (* Degrees of freedom for estimate of S^2 *)
jointPDF = (Exp[-\[Mu] . Inverse[r] . \[Mu]/(2 \[Sigma]^2)]/(\[Pi] f Sqrt[ Det[r]] Gamma[f/2]))*
Sum[(2^(\[Alpha]/2) (t . Inverse[r] . \[Mu])^\[Alpha]  Gamma[(f + 2 + \[Alpha])/2])/
(\[Sigma]^\[Alpha] f^(\[Alpha]/2) \[Alpha]!
(1 + t . Inverse[r] . t/f)^((f + 2 + \[Alpha])/2)), {\[Alpha], 0, \[Infinity]}];
jointPDF = FullSimplify[jointPDF, Assumptions -> {-1 < \[Rho] < 1, \[Sigma] > 0, n > 1,
n \[Element] Integers, \[Mu]x \[Element] Reals, \[Mu]y \[Element]
Reals, x \[Element] Reals, y \[Element] Reals}]


A more readable version of the code is below:

The result is

The pdf of the difference $$Z = X/S - Y/S$$ can be found numerically by replacing $$y$$ with $$x-z$$ and then integrating over $$x$$:

(* Numerical estimate of pdf of X/S - Y/S for a few values of n \
(sample size for estimating \[Sigma])*)
pdfz100 =
Table[{z,
NIntegrate[
jointPDF /. {y -> x - z,
n -> 100, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y ->
3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
pdfz4 = Table[{z,
NIntegrate[
jointPDF /. {y -> x - z,
n -> 4, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y ->
3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
pdfz2 = Table[{z,
NIntegrate[
jointPDF /. {y -> x - z,
n -> 2, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y ->
3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
ListPlot[{pdfz100, pdfz4, pdfz2}, Joined -> True, ImageSize -> Large,
PlotLegends -> {"n = 100", "n = 4", "n = 2"},
PlotLabel ->
Style["\[Sigma] = 2, \[Rho] = 1/2, \!$$\*SubscriptBox[\(\[Mu]$$, \
$$x$$]\) = 1, \!$$\*SubscriptBox[\(\[Mu]$$, $$y$$]\) = 3", Bold, 18]]


The results follow:

As a check one should perform some simulations.

parms = {\[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y -> 3, n -> 2};
nsim = 100000; (* Number of simulations *)
(* Data for x and y *)
data = RandomVariate[
BinormalDistribution[{\[Mu]x, \[Mu]y}, {\[Sigma], \[Sigma]}, \[Rho]] /. parms, nsim];
(* Data to for estimating S *)
xy = RandomVariate[
BinormalDistribution[{mu1, mu2}, {\[Sigma], \[Sigma]}, \[Rho]] /.
parms, {nsim, n /. parms}];
s = Sqrt[(Variance[#[[All, 1]]]/2 + Variance[#[[All, 2]]]/2) & /@ xy];
(* Z = X/S - Y/S *)
zz = data[[All, 1]]/s - data[[All, 2]]/s;

(* Numerically estimate the pdf of z *)
pdfz = Table[{z, NIntegrate[jointPDF /. y -> x - z /. parms,
{x, -\[Infinity], \[Infinity]}]},
{z, Quantile[zz, 0.005], Quantile[zz, 0.995], (Quantile[zz, 0.995] - Quantile[zz, 0.005])/200}];

(* Plot the results *)
Show[Histogram[zz, "FreedmanDiaconis", "PDF"],
ListPlot[pdfz, Joined -> True, PlotRange -> All]]


There seems to be a match.