1
$\begingroup$

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)$.

$\endgroup$
8
  • $\begingroup$ 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.) $\endgroup$
    – whuber
    Mar 10, 2022 at 23:25
  • $\begingroup$ I edited the question. $\endgroup$ Mar 11, 2022 at 6:51
  • $\begingroup$ 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$. $\endgroup$
    – JimB
    Mar 11, 2022 at 7:03
  • $\begingroup$ 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)$ ? $\endgroup$
    – JimB
    Mar 11, 2022 at 19:40
  • $\begingroup$ @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)$. $\endgroup$ Mar 11, 2022 at 20:01

1 Answer 1

1
$\begingroup$

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:

Code to obtain joint distribution

The result is

joint distribution of X/S and Y/S

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]]

Again, a more readable version:

Code to find pdf of z

The results follow:

pdfs for various parameters

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]]

Results of simulations

There seems to be a match.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.