I'm trying to use TransformedDistribution to find the distribution of a function of two variables that have Normal distributions. The function is $h(x, y)=x^2+y^2$.

I tried to do this two ways, first by taking intermediate steps, and then, by doing it all at once:

  1.  q1 = TransformedDistribution[x^2, x \[Distributed] NormalDistribution[0, Sqrt[Vxx]]];
     q2 = TransformedDistribution[y^2, y \[Distributed] NormalDistribution[0, Sqrt[Vyy]]];
     W = TransformedDistribution[x + y, {x \[Distributed] q1, y \[Distributed] q2}];
  2. PDF[TransformedDistribution[x^2 + y^2, {x \[Distributed] NormalDistribution[0, Sqrt[Vxx]], y \[Distributed] NormalDistribution[0, Sqrt[Vyy]]}]]

Both have the same type of results:

    out1: PDF[TransformedDistribution[x + y, {x \[Distributed] TransformedDistribution[\[FormalX]^2, \[FormalX] \[Distributed] NormalDistribution[0, Sqrt[Vxx]]], y \[Distributed] TransformedDistribution[\[FormalX]^2, \[FormalX] \[Distributed] NormalDistribution[0, Sqrt[Vyy]]]}]]
    out2: PDF[TransformedDistribution[\[FormalX]1^2 + \[FormalX]2^2, \{\[FormalX]1 \[Distributed] NormalDistribution[0, Sqrt[Vxx]], \[FormalX]2 \[Distributed] NormalDistribution[0, Sqrt[Vyy]]}]]

Is there a way to get Mathematica to return a function?


Sometimes you just need to do a brute force approach to get the pdf.

pdfx = Simplify[PDF[TransformedDistribution[x^2, x \[Distributed] NormalDistribution[0, σ]], x],
   Assumptions -> {σ > 0, x > 0}];
pdf = Integrate[(pdfx /. {x -> x1, σ -> σ1}) (pdfx /. {x -> z - x1, σ -> σ2}), {x1, 0, z}, 
  Assumptions -> {z > 0, σ1 > 0, σ2 > 0}]
(* (E^(-(1/4) z (1/σ1^2 + 1/σ2^2)) BesselI[0, 1/4 z (-(1/σ1^2) + 1/σ2^2)])/(2 σ1 σ2) *)

As a check one should do some simulations:

n = 10000;
σ10 = 1;
σ20 = 3;
xx1 = RandomVariate[NormalDistribution[0, σ10], n];
xx2 = RandomVariate[NormalDistribution[0, σ20], n];
sum = xx1^2 + xx2^2;
skd = SmoothKernelDistribution[sum, Automatic, {"Bounded", {0, ∞}, "Gaussian"}];
Plot[{PDF[skd, z], pdf /. {σ1 -> σ10, σ2 -> σ20}}, {z, 0, 40}, PlotRange -> All, 
  PlotLegends -> Placed[{"Estimated density", "True density"}, {Right, Center}]]

Estimated pdf and true pdf

I don't think the cdf has a nice closed-form (but I didn't spend much time on that).

For moments the functions MomentGeneratingFunction and CentralMomentGeneratingFunction work fine in this case:

dist = TransformedDistribution[x1^2 + x2^2, 
 {x1 \[Distributed] NormalDistribution[0, σ1], 
  x2 \[Distributed] NormalDistribution[0, σ2]}];

(* σ1^2 + σ2^2 *)

(* 2 (σ1^4 + σ2^4) *)

mgf = MomentGeneratingFunction[dist, t]
(* 1/(Sqrt[1 - 2 t σ1^2] Sqrt[1 - 2 t σ2^2]) *)
D[mgf, {t, 2}] /. t -> 0
(* 3 σ1^4 + 2 σ1^2 σ2^2 + 3 σ2^4 *)

cmgf = CentralMomentGeneratingFunction[dist, t]
(* E^(-t (σ1^2 + σ2^2))/(Sqrt[1 - 2 t σ1^2] Sqrt[1 - 2 t σ2^2]) *)
D[cmgf, {t, 2}] /. t -> 0 // Expand
(* 2 σ1^4 + 2 σ2^4 *)
  • $\begingroup$ The function that you wrote at the bottom of the first paragraph works perfectly. But, if I try to get the same result by running the code above, the result is Integrate[...], the same problem as with PDF[...]. Why doesn't it return a function? $\endgroup$
    – Anthill
    Mar 8 at 17:31
  • $\begingroup$ I was using version Windows 10. What version and operating system are you using? $\endgroup$
    – JimB
    Mar 8 at 19:36
  • $\begingroup$ I'm using version, on Windows 10 $\endgroup$
    – Anthill
    Mar 9 at 15:26

This is not an answer but a possible partial workaround.

First note that

     W = TransformedDistribution[
   x^2 + y^2, {x \[Distributed] NormalDistribution[0, \[Sigma]], 
    y \[Distributed] NormalDistribution[0, \[Sigma]]}];
PDF[W, x]

enter image description here

While the PDF is not returned by Mathematica you can get the generating function:

W = TransformedDistribution[
   x^2 + y^2, {x \[Distributed] NormalDistribution[\[Mu]1, \[Sigma]1],
     y \[Distributed] NormalDistribution[\[Mu]2, \[Sigma]2]}];


\[Phi][t_] = Expectation[Exp[I t  x], x \[Distributed] W]


enter image description here

Interestingly with $\mu1=\mu2=0$ as OP, this allows us to compute the nth moment analytically so to speak:

D[\[Phi][t], {t, n}]/I^n // FunctionExpand

enter image description here

Note that the $(\mu1,\mu2)\neq(0,0)$ also yields a formal solution for arbitrary $n$ through a 4th order recursion relation.

enter image description here

  • $\begingroup$ Mathematica returns a function for the PDFs of q1 and q2, the problem is it doesn't for W. When I try your suggestion for W, I get the same result: PDF[W, x], and not a function $\endgroup$
    – Anthill
    Mar 8 at 13:06
  • $\begingroup$ Yes, this works for $Vxx = Vyy$, but it doesn't when they have different values $\endgroup$
    – Anthill
    Mar 8 at 13:07

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.