# PDF of $x_1^2+x_2^2-x_3^2-x_4^2$ with $x_i \sim N(\mu_i,1)$

What is the probability distribution function of the random variable $$X$$

$$X=x_1^2+x_2^2-x_3^2-x_4^2, \tag{1}$$

where $$x_i$$ are independent normally distributed variables $$\mathcal{N}(\mu_i,1)$$?

What I tried

$$X$$ can also be expressed as the difference of two noncentral $$\chi^2(k,\lambda)$$ distributed variables:

$$X=\chi^2(2,x_1^2+x_2^2)-\chi^2(2,x_3^2+x_4^2)=\chi^2(2,\lambda_1)-\chi^2(2,\lambda_2).\tag{2}$$

If we try to derive the pdf by the characteristic function of the noncentral $$\chi^2$$-distribution we encounter an integral of the form

$$f(z)=\int_{-\infty}^{\infty}\frac{{\rm exp}\left(i t z\right)}{3t^2+1}{\rm exp}\left(\frac{i t a+ t^2 b}{3t^2+1}\right){\rm d}t\\ {\rm with }\,\, z\in \mathbb{R},a\in \mathbb{R},b\in \mathbb{R}_{< 0}, i^2=-1,\tag{3}$$

that I could not solve.

In this forum I found only this related post however it is unclear if the answer given there is applicable here, or if the answer is right.

• I think because we have $y_i^2 \sim -x_i^2$ and these are independent then the distribution should be noncentral $\chi^2(4, \sum_i \mu_i^2)$, but i'm not 100% sure on that Commented Jan 24, 2022 at 21:34
• @bdeonovic That cannot possibly be correct, because $X$ has positive probability of being negative whereas that is never the case for any $\chi^2$ distribution.
– whuber
Commented Jan 24, 2022 at 21:40
• yep I just realized my own stupidity, need to get more coffee. I appreciate the bluntness of your correction though @whuber, in contrast to your typical gentle leading question :) Commented Jan 24, 2022 at 21:40
• @bdeonovic Sorry. I make those kinds of mistakes all the time--I blame haste ;-).
– whuber
Commented Jan 24, 2022 at 21:43
• You can numerically compute it / plot it in R using distr::plot(distr::Chisq(2, 2) - distr::Chisq(2, 2), to.draw.arg='d'). Commented Jan 28, 2022 at 8:01

As noted by the OP, this problem can be simplified to the difference in two independent noncentral $$\chi^2$$ random variables:

$$X=X_1^2+X_2^2-X_3^2-X_4^2=Z_1-Z_2$$

where $$Z_i\sim \chi^2(2,\lambda_i)$$, $$\lambda_1=\mu_1^2+\mu_2^2$$, $$\lambda_2=\mu_3^2+\mu_4^2$$ and $$Z_1$$ and $$Z_2$$ are independent.

Because your profile shows you have access to Mathematica, here is how one can find any desired moment of $$Z$$:

dist = TransformedDistribution[z1 - z2,
{z1 \[Distributed] NoncentralChiSquareDistribution[2, λ1],
z2 \[Distributed] NoncentralChiSquareDistribution[2, λ2]}];
Mean[dist]
(* λ1-λ2 *)
Variance[dist]
(* 4 (2+λ1+λ2) *)
Skewness[dist]
(* (3 (λ1-λ2))/(2+λ1+λ2)^(3/2) *)
Kurtosis[dist]
(* 3 (8+8 λ1+λ1^2+8 λ2+2 λ1 λ2+ *)
CentralMoment[dist, 5]
(* 960 (4 λ1+λ1^2-4 λ2-λ2^2) *)
CentralMoment[dist, 6]
(* 960 (48+72 λ1+24 λ1^2+λ1^3+72 λ2+24 λ1 λ2+3 λ1^2 λ2+24 λ2^2+3 λ1 λ2^2+λ2^3) *)
CentralMoment[dist, 7]
(* 40320 (24 λ1+12 λ1^2+λ1^3-24 λ2+λ1^2 λ2-12 λ2^2-λ1 λ2^2-λ2^3) *)


Other than maybe some special cases, it is unlikely that there's a nice closed-form for the density of $$Z$$. But to estimate the density numerically (and perform a check with a large random sample), consider the following for $$\lambda_1=\mu_1^2+\mu_2^2=1^2+2^2=5$$ and $$\lambda_2=\mu_3+\mu_4=1^2+\sqrt{2}^2=3$$:

pdf1 = PDF[NoncentralChiSquareDistribution[2, λ1], z1][[1, 1, 1]] // FunctionExpand


pdf2 = PDF[NoncentralChiSquareDistribution[2, λ2], z1 - z][[1, 1, 1]] // FunctionExpand


pdfz[z_?NumericQ, λ1_?NumericQ, λ2_?NumericQ] :=
NIntegrate[1/2 E^(1/2 (-z1 - λ1)) BesselI[0, Sqrt[z1] Sqrt[λ1]]*
1/2 E^(1/2 (z - z1 - λ2)) BesselI[0, Sqrt[-z + z1] Sqrt[λ2]],
{z1, Max[0, z], ∞}]

(* Check on a specific example *)
parms = {μ1 -> 1, μ2 -> 2, μ3 -> 1, μ4 -> Sqrt[2]};

(* Distribution of Z *)
dist = TransformedDistribution[z1 - z2,
{z1 \[Distributed] NoncentralChiSquareDistribution[2, μ1^2 + μ2^2],
z2 \[Distributed] NoncentralChiSquareDistribution[2, μ3^2 + μ4^2]}] /. parms;

(* Random sample from Z *)
zz = RandomVariate[dist, 1000000];

(* Plot density and density histogram *)
m = Mean[dist];
s = StandardDeviation[dist];
Show[Histogram[zz, "FreedmanDiaconis", "PDF", Frame -> True, FrameLabel -> {"Z", "Density"}]],
Plot[pdfz[z, μ1^2 + μ2^2 /. parms, μ3^2 + μ4^2 /. parms], {z, m - 4 s, m + 4 s}]]


• What is the advantage of this method compared to eq.3? Commented Dec 31, 2023 at 4:36
• (a) As there is likely no nice closed-form for (3) both approaches need to use a numerical approximation. (b) (3) requires dealing with $i$ and not all software packages can deal with that, (c) Even in Mathematica you'd need to use the Chop function to get rid of "little imaginary bits" left over from numerical imprecision, (d) $a$ and $b$ are not defined as of yet. (And, yes, (3) is a "compact closed-form" but not easily evaluated in this case. I wish it were otherwise.)
– JimB
Commented Dec 31, 2023 at 4:58
• Actually, the issue with $i$ in (3) can be removed and things become much more numerically stable. I'll add that as a separate answer in the next few days - unless someone else wants to do that. Still the density needs to be determined numerically.
– JimB
Commented Dec 31, 2023 at 6:38

Here is the inverse Fourier transform approach you've labeled as equation (3). (Again, I've used Mathematica as your profile indicates you have access to it and it does all of the algebraic manipulations for me.)

(* Characteristic functions of Z1 and -Z2 *)
cf1 = CharacteristicFunction[NoncentralChiSquareDistribution[2, λ1], t]
cf2 = CharacteristicFunction[NoncentralChiSquareDistribution[2, λ2], -t]
(* Characteristic function of Z = Z1 - Z2 *)
cf = cf1 cf2 // Simplify


$$\frac{e^{t \left(\frac{\text{\lambda 2}}{-2 t+i}-\frac{\text{\lambda 1}}{2 t+i}\right)}}{4 t^2+1}$$

(* Determine integrand for the inverse Fourier transform *)
(* https://stats.stackexchange.com/questions/104491/how-to-find-a-density-from-a-characteristic-function *)
(* Drop the imaginary term as it will integrate to zero *)
integrand = ComplexExpand[ReIm[Exp[-I z t] cf1 cf2]][[1]]/(2 π) // Together // FullSimplify


$$\frac{e^{-\frac{2 t^2 (\text{\lambda 1}+\text{\lambda 2})}{4 t^2+1}} \cos \left(\frac{t \left(-\text{\lambda 1}+\text{\lambda 2}+4 t^2 z+z\right)}{4 t^2+1}\right)}{8 \pi t^2+2 \pi }$$

(* Define pdf function using numerical integration as there is likely no closed-form solution *)
pdf[z_?NumericQ, λ1_?NumericQ, λ2_?NumericQ] :=
NIntegrate[(E^(-((2 t^2 (λ1 + λ2))/(1 + 4 t^2))) Cos[(t (z + 4 t^2 z - λ1 + λ2))/(1 + 4 t^2)])/
(2 π + 8 π t^2), {t, -∞, ∞}]

(* Set parameters *)
parms = {λ1 -> 5, λ2 -> 3};

(* Determine mean and standard deviation of Z for setting plotting limits *)
dist = TransformedDistribution[z1 - z2,
{z1 \[Distributed] NoncentralChiSquareDistribution[2, λ1],
z2 \[Distributed] NoncentralChiSquareDistribution[2, λ2]}] /. parms;
m = Mean[dist];
s = StandardDeviation[dist];

(* Generate a large random sample as a check *)
zz = RandomVariate[dist, 100000];

(* Plot density and histogram *)
Show[Histogram[zz, "FreedmanDiaconis", "PDF",
Frame -> True, FrameLabel -> {"Z", "Density"}],
Plot[Evaluate[pdf[z, λ1 /. parms, λ2 /. parms]], {z, m - 4 s, m + 4 s}]]