Given are 4 independent normal variables $x_i \sim \mathcal{N}(\mu_i,\sigma_i)_{i=1,\ldots,4}$, that are used to define the random variable $$Z\sim (x_1 - x_2 + x_3 + x_4)^2 - (-x_1 + x_2 + x_3 + x_4)^2.\tag{1}$$

Each squared term in eq.(1) is a linear combination of independent normal variables, whence each of those squares has a noncentral $\chi^2(k,M)$-distribution with $k=1$ degree of freedom and noncentrality parameter $M$. The two linear combinations are also orthogonal, $\begin{pmatrix} 1 & -1 & 1&1 \end{pmatrix} \cdot \begin{pmatrix} -1 & 1 & 1&1 \end{pmatrix}^T =0$ , whence uncorrelated; and because they are uncorrelated joint random variables, they are independent (see here).

Thus the squared linear combinations in eq.(1) can be written as a difference of $\chi^2$-variables $$Z\sim s \left( \chi^2\left(1,M_1\right)- \chi^2\left(1,M_2\right)\right)\tag{2}$$

with noncentrality parameters $M_1=\dfrac{1}{s}\left(\mu_1-\mu_2+\mu_3+\mu_4\right)^2, M_2=\dfrac{1}{s}\left(-\mu_1+\mu_2+\mu_3+\mu_4\right)^2$ and total variance $s=\sum_{i=1}^4\sigma^2_i$.

However this approach only works, if all $\sigma_i$ are identical. How to express eq.(1) by $\chi^2$-variables for different $\sigma_i$?

Graphical examples for PDF

All $\sigma_i$ are identical.

$\sigma=(5,5,5,5), \mu=(2,-11,-1,4)$, black: eq(1), red: eq(2) (both plots almost overlap)

All $\sigma_i$ are different.

$\sigma=(5,4,3,2), \mu=(2,-11,-1,4)$, black: eq(1), red: eq(2)

  • $\begingroup$ Short answer: $Z$ can be expressed as a quadratic form in four iid standard Normal variables. Diagonalize it, reducing the question to a weighted sum of squares of non-central chi-squared variates: you're done. $\endgroup$
    – whuber
    Commented Dec 26, 2023 at 15:18
  • $\begingroup$ Even shorter. Let $Z_1=4(X_1-X_2)$ and $Z_2=X_3+X_4$. $Z_1$ and $Z_2$ are independent normals with respective means $4(\mu_1-\mu_2)$ and $\mu_3+\mu_4$ and respective variances $16(\sigma_1^2+\sigma_2^2)$ and $\sigma_3^2+\sigma_4^2$. Now $Z=(Z_1/4+Z_2)^2-(-Z_1+Z_2)^2=Z_1 Z_2$. Now you're just dealing with the product of two independent normals. $\endgroup$
    – JimB
    Commented Dec 26, 2023 at 18:46
  • $\begingroup$ @JimB: Right -- but it's "even shorter" only for this special case. Generally, a quadratic form will be expressible as a linear combination of noncentral chi-squared variates, but not as a product of Normal variates. $\endgroup$
    – whuber
    Commented Dec 26, 2023 at 20:37
  • $\begingroup$ @whuber I will amend my answer with that in mind. $\endgroup$
    – JimB
    Commented Dec 26, 2023 at 20:46

3 Answers 3


This problem is initially written as the difference in the squares of two correlated normals.


where all of the $X_i$ are independent random variables with potentially different means and variances: $X_i \sim N(\mu_i,\sigma_i^2)$.

The request is to re-write $Z$ in terms of the difference in two independent noncentral chisquare random variables times some constant:


where $W_i\sim \chi^2(1,m_i)$. To do so we do something similar to the idea presented in @whuber 's comment: The characteristic functions are found for both characterizations of $Z$ and then we equate all of the coefficients for the characteristic function argument $t$.

To have a less cluttered expression for


we first re-write $Z$ as follows:

$$Z=(Z_{12}+Z_{34})^2 - (-Z_{12}+Z_{34})^2$$

where $Z_{12}=X_1-X_2$ and $Z_{34}=X_3+X_4$. $Z_{12}$ and $Z_{34}$ are independent normals with respective means and variances $\mu_{12}=\mu_1-\mu_2$, $\mu_{34}=\mu_3+\mu_4$, $\sigma_{12}^2=\sigma_1^2+\sigma_2^2$, and $\sigma_{34}^2=\sigma_3^2+\sigma_4^2$.

Written this way the characteristic function of $Z$ is found with

dist = TransformedDistribution[(z[12] + z[34])^2 - (-z[12] + z[34])^2,
    {z[12] \[Distributed] NormalDistribution[μ[12], σ[12]],
     z[34] \[Distributed] NormalDistribution[μ[34], σ[34]]}];
cf = CharacteristicFunction[dist, t]

$$\frac{e^ {-\frac{4 t \left(-i \mu _{12} \mu _{34}+2 \mu _{34}^2 \sigma _{12}^2 t+2 \mu _{12}^2 \sigma _{34}^2 t\right)}{16 \sigma _{12}^2 \sigma _{34}^2 t^2+1}}}{\sqrt{16 \sigma _{12}^2 \sigma _{34}^2 t^2+1}}=\frac{e^{\frac{-t^2 \left(8 \mu _{34}^2 \sigma _{12}^2+8 \mu _{12}^2 \sigma _{34}^2\right)+i 4 \mu _{12} \mu _{34} t}{16 \sigma _{12}^2 \sigma _{34}^2 t^2+1}}}{\sqrt{16 \sigma _{12}^2 \sigma _{34}^2 t^2+1}}$$

The characteristic function for the alternate formulation is found with the following:

dist2 = TransformedDistribution[s (w[1] - w[2]),
   {w[1] \[Distributed] NoncentralChiSquareDistribution[1, m[1]],
    w[2] \[Distributed] NoncentralChiSquareDistribution[1, m[2]]}];
cf2 = CharacteristicFunction[dist2, t];

$$\frac{e^{-s t \left(\frac{m_1}{2 s t+i}+\frac{m_2}{2 s t-i}\right)}}{\sqrt{4 s^2 t^2+1}}=\frac{e^{\frac{-2 \left(m_1+m_2\right) s^2 t^2+i \left(m_1-m_2\right) s t}{4 s^2 t^2+1}}}{\sqrt{4 s^2 t^2+1}}$$

We can find the matching terms of $t$ and $t^2$ in the exponents to determine $m_1$, $m_2$, and $s$. This means we need to simultaneously solve

$$16 \sigma^2_{12} \sigma^2_{34}=4s^2$$ $$8 \mu^2_{34}\sigma^2_{12}+8\mu^2_{12}\sigma^2_{34}=2s^2(m_1+m_2)$$ $$4\mu_{12}\mu_{34}=s(m_1-m_2)$$

sol = Solve[{16 σ[12]^2 σ[34]^2 == 4 s^2, 
  8 μ[34]^2 σ[12]^2 + 8 μ[12]^2 σ[34]^2 == 2 s^2 (m[1] + m[2]), 
  4 μ[12] μ[34] == s (m[1] - m[2])}, {s, m[1], m[2]}]

This results in two solutions:

$$s= \pm 2 \sigma _{12} \sigma _{34}$$ $$m_1= \frac{\left(\mu _{34} \sigma _{12}\pm\mu _{12} \sigma _{34}\right){}^2}{2 \sigma _{12}^2 \sigma _{34}^2}$$ $$m_2= \frac{\left(\mu _{34} \sigma _{12}\mp\mu _{12} \sigma _{34}\right){}^2}{2 \sigma _{12}^2 \sigma _{34}^2}$$

Finally, we can plug in the original parameters:

parms = {μ[12] -> μ[1] - μ[2], μ[34] -> μ[3] + μ[4], σ[12] -> Sqrt[σ[1]^2 + σ[2]^2],  σ[34] -> Sqrt[σ[3]^2 + σ[4]^2]}
sol = sol //. parms

$$s= \pm 2 \sqrt{\sigma _1^2+\sigma _2^2} \sqrt{\sigma _3^2+\sigma _4^2}$$ $$m_1= \frac{\left(\left(\mu _2-\mu _1\right) \sqrt{\sigma _3^2+\sigma _4^2}\pm\left(\mu _3+\mu _4\right) \sqrt{\sigma _1^2+\sigma _2^2}\right){}^2}{2 \left(\sigma _1^2+\sigma _2^2\right) \left(\sigma _3^2+\sigma _4^2\right)}$$ $$m_2= \frac{\left(\left(\mu _1-\mu _2\right) \sqrt{\sigma _3^2+\sigma _4^2}\mp \left(\mu _3+\mu _4\right) \sqrt{\sigma _1^2+\sigma _2^2}\right){}^2}{2 \left(\sigma _1^2+\sigma _2^2\right) \left(\sigma _3^2+\sigma _4^2\right)}$$

As a check consider one of the given examples.

parms = {μ[1] -> 2, μ[2] -> -11, μ[3] -> -1, μ[4] -> 4, σ[1] -> 5,
  σ[2] -> 4, σ[3] -> 3, σ[4] -> 2};
parms1 = {μ[12] -> μ[1] - μ[2], μ[34] -> μ[3] + μ[4], σ[12] -> Sqrt[σ[1]^2 + σ[2]^2],
  σ[34] -> Sqrt[σ[3]^2 + σ[4]^2]};
z1 = RandomVariate[dist /. parms1 /. parms, 1000000];

dist2 = TransformedDistribution[
  s (w[1] - w[2]) /. sol[[2]] /. parms1 /. parms, 
  {w[1] \[Distributed] NoncentralChiSquareDistribution[1, m[1]], 
   w[2] \[Distributed] NoncentralChiSquareDistribution[1, m[2]]} /. sol[[2]] /. parms1 /. parms];
z2 = RandomVariate[dist2 /. parms, 1000000];

Histogram[z1, "FreedmanDiaconis", "PDF"]

Histogram of samples from original distribution

Histogram[z2, "FreedmanDiaconis", "PDF"]

Histogram of transformed distribution

We can also check on the corresponding sample means and variances:

{Mean[#], Variance[#]} & /@ {z1, z2}
(* {{156.019, 49649.8}, {156.27, 49587.}} *)

The pdf can also be easily found numerically using the characteristic function for specified values of the parameters.

Answer to a comment below:

It was asked to find the distribution and characteristic function of a more complicated random variable:


Rather than directly using that definition, it appears that Mathematica needs a little help.

First determine the multivariate normal distribution of the 4 differences in normal random variables.

dist = TransformedDistribution[{z[5] - z[3], z[2] - z[4], z[1] - z[3], z[6] - z[4]},
  {z[1] \[Distributed] NormalDistribution[μ[1], σ[1]], 
   z[2] \[Distributed] NormalDistribution[μ[2], σ[2]], 
   z[3] \[Distributed] NormalDistribution[μ[3], σ[3]], 
   z[4] \[Distributed] NormalDistribution[μ[4], σ[4]], 
   z[5] \[Distributed] NormalDistribution[μ[5], σ[5]], 
   z[6] \[Distributed] NormalDistribution[μ[6], σ[6]]}]

(*  MultinormalDistribution[{-μ[3] + μ[5], μ[2] - μ[4], μ[1] - μ[3], -μ[4] + μ[6]}, 
  {{σ[3]^2 + σ[5]^2, 0,               σ[3]^2,           0}, 
   {0,               σ[2]^2 + σ[4]^2, 0,                σ[4]^2}, 
   {σ[3]^2,          0,               σ[1]^2 + σ[3]^2,  0}, 
   {0,               σ[4]^2,          0,                σ[4]^2 + σ[6]^2}}]*)

Now set up rules such that the parameters in the above distribution each have a single symbol. This seems to be necessary for Mathematica to complete in a reasonable amount of time.

parms1 = {-μ[3] + μ[5] -> μ1, μ[2] - μ[4] -> μ2, μ[1] - μ[3] -> μ3, -μ[4] + μ[6] -> μ4,
  σ[3]^2 + σ[5]^2 -> v1, σ[2]^2 + σ[4]^2 -> v2, σ[1]^2 + σ[3]^2 -> v3, 
  σ[4]^2 + σ[6]^2 -> v4, σ[3]^2 -> v5, σ[4]^2 -> v6};

(* Find distribution and characteristic function of the simplified difference 
   of the products of normally distributed random variables *)
dist2 = TransformedDistribution[x[1] x[2] - x[3] x[4],
  {x[1], x[2], x[3], x[4]} \[Distributed] dist] /.parms1;
AbsoluteTiming[cf = CharacteristicFunction[dist2, t]]

The result is messy and takes about 4 minutes to run. Now replace with the original parameters:

parms2 = parms1 /. {Rule[a_, b_] -> Rule[b, a]}
cf = cf /. parms2

That result is even messier so I won't show it here.

  • $\begingroup$ Sorry for the confusion. An earlier version used $Z_1$ and $Z_2$ but in an attempt to make things clearer, I switched those to $Z_{12}$ and $Z_{34}$ but obviously didn't make all of the changes in the text. I'm in the process of simplifying my answer and should have that done later today. $\endgroup$
    – JimB
    Commented Jan 12 at 17:29
  • $\begingroup$ What case is taking a long time? The characteristic functions cf and cf2 take 13 seconds and 0.03 seconds, respectively. $\endgroup$
    – JimB
    Commented Jan 16 at 4:26
  • $\begingroup$ If one rewrites that random variable as$ x_1 x_2-x_3 x_4$ where $(x_1,x_2,x_3,x_4)$ has multiivariate normal distribution, then the characteristic functions takes about 6 minutes to complete. Sometimes Mathematica needs a little help. (In other words, $x_1=z_5-x_3$, $x_2=z_2-z_4$, etc.) $\endgroup$
    – JimB
    Commented Jan 18 at 4:10
  • $\begingroup$ Sorry about that. Yes, $x_1=z_5-x_3$ should be $x_1=z_5-z_3$. $\endgroup$
    – JimB
    Commented Jan 18 at 16:31
  • 1
    $\begingroup$ The assignments σ[3]^2 -> v5, σ[4]^2 -> v6 in parms1 seem to be either wrong, incomplete or superficial or I do not understand its meaning. $\endgroup$ Commented Jan 24 at 23:30

For this special case the random variable $Z$ can be simplified to the product of two independent normals:

Let $Z_1=4(X_1-X_2)$ and $Z_2=X_3+X_4$. $Z_1$ and $Z_2$ are independent normals with respective means $4(\mu_1-\mu_2)$ and $\mu_3+\mu_4$ and respective variances $16(\sigma_1^2+\sigma_2^2)$ and $\sigma_3^2+\sigma_4^2$. Now

$$Z=(Z_1/4+Z_2)^2-(-Z_1/4+Z_2)^2=Z_1 Z_2$$

Finding moments of $Z$ should be straightforwared. To plot the pdf of $Z$ using Mathematica (which your profile says you have access), one can execute the following:

(* Numerical estimate of the pdf of Z *)
pdf[z_?NumericQ, μ1_?NumericQ, σ1_?NumericQ, μ2_?NumericQ, σ2_?NumericQ] := 
 NIntegrate[PDF[NormalDistribution[μ1, σ1], x] PDF[NormalDistribution[μ2, σ2], z/x]/Abs[x],
 {x, -∞, ∞}]

(* Parameters from your example *)
σz1 = Sqrt[16 (25 + 16)]
σz2 = Sqrt[9 + 4]
μz1 = 4 (2 - (-11))
μz2 = -1 + 4

(* Plot of density *)
Plot[pdf[z, μz1, σz1, μz2, σz2], {z, -400, 1000}]

PDF of Z using OP's example parameters

Expressing Z as a linear combination of two independent noncentral chisquare rv's

Define the distribution of the linear combination.


dist = TransformedDistribution[a1 y1 + a2 y2,
   {y1 \[Distributed] NoncentralChiSquareDistribution[1, m1], 
    y2 \[Distributed] NoncentralChiSquareDistribution[1, m2]}];

Define the distribution of Z expressed as the product of two independent normal rv's.

distz1z2 = TransformedDistribution[z1 z2,
  {z1 \[Distributed] NormalDistribution[μz1, σz1], 
   z2 \[Distributed] NormalDistribution[μz2, σz2]}];

If we can characterize $Z$ as a linear combination of two independent noncentral $\chi^2$ rv's which has 4 parameters, then the first 4 moments need to match with that of the characterization of $Z$ as the product of two independent normals.

sol = Solve[{Mean[distz1z2] == Mean[dist], 
     Variance[distz1z2] == Variance[dist],
     Moment[distz1z2, 3] == Moment[dist, 3], 
     Moment[distz1z2, 4] == Moment[dist, 4]},
    {a1, a2, m1, m2}][[1]] // FullSimplify

Solution to parameters of linear combination

From above we know that

parms2 = {σz1 -> Sqrt[16 (σ1^2 + σ2^2)], σz2 -> Sqrt[σ3^2 + σ4^2], 
  μz1 -> 4 (μ1 - μ2), μz2 -> μ3 + μ4};

So in terms of the original parameters the parameters of the linear combination of independent noncentral chisquare rv's are

sol /. parms2

Solution in terms of original parameters

As a minimal check, consider your example.

parms = {μ1 -> 2, μ2 -> -11, μ3 -> -1, μ4 -> 4, σ1 -> 5, σ2 -> 4, σ3 -> 3, σ4 -> 2};
zA = RandomVariate[distz1z2 /. sol /. parms2 /. parms, 1000000];
zB = RandomVariate[dist /. sol /. parms2 /. parms, 1000000];
Histogram[zA, "FreedmanDiaconis", "PDF", PlotRange -> {{-700, 1300}, {0, 0.003}}]

Histogram of density for product of two independent normals

Histogram[zB, "FreedmanDiaconis", "PDF", PlotRange -> {{-700, 1300}, {0, 0.003}}]

Histogram of density for linear combination of two noncentral chisquare random variables

  • $\begingroup$ Can the answer solely be expressed by $\chi^2$ variables as desired by the OP? $\endgroup$ Commented Dec 27, 2023 at 20:55
  • $\begingroup$ I know that's what you asked for but what I've presented is an approach to obtain the exact distribution. However, I'll take a look at if the distribution of $Z$ can be expressed as a linear combination of two independent noncentral $\chi^2$ variables. If so, that will just be equivalent to characterizing the distribution of $Z$ as the product of two independent normals. $\endgroup$
    – JimB
    Commented Dec 27, 2023 at 23:11
  • $\begingroup$ How to solve the problem without Mathematica? For this case, accidentially, after 5 minutes of calculation, Mathematica returns a result. But we do not learn anything about the steps and therefore do not know the limit of the method. $\endgroup$ Commented Dec 28, 2023 at 13:59
  • $\begingroup$ “accidentally” ? Ha, ha. You have the answer. Like is too short for me to spend time on such things. And get a faster computer. $\endgroup$
    – JimB
    Commented Dec 28, 2023 at 14:51
  • $\begingroup$ How can Mathematica solve this if we instead of the product distribution of 2 normals use the original form (eq.1) as input for the code. $\endgroup$ Commented Jan 2 at 21:15

This answer only explains why the approach fails but does not say how to get the expected result.

Set $$y_i\sim \mathcal{N}\left(\dfrac{\mu_i}{\sigma_i},1\right)$$ and write $$Z\sim(\sigma_1 y_1-\sigma_2 y_2+\sigma_3 y_3+\sigma_4 y_4)^2-(-\sigma_1 y_1+ \sigma_2 y_2+\sigma_3 y_3+\sigma_4 y_4)^2.$$ You see that the linear combinations are not orthogonal but dependent. They are independent only for the special case $\sigma^2_1+\sigma^2_2=\sigma^2_3+\sigma^2_4$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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