# TransformedDistribution not returning a function

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}];
PDF[W]

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?

• The problem seems to be the symbolic variance. Is it acceptable to use numerical values. Mar 8 at 14:23
• You might want to look at stats.stackexchange.com/questions/67533/….
– JimB
Mar 8 at 17:14

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;
SeedRandom[12345];
σ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}]]

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

Mean[dist]
(* σ1^2 + σ2^2 *)

Variance[dist]
(* 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 *)
• 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? Mar 8 at 17:31
• I was using version 12.2.0.0 Windows 10. What version and operating system are you using?
– JimB
Mar 8 at 19:36
• I'm using version 11.3.0.0, on Windows 10 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]

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

Then

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

returns

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

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

• 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 Mar 8 at 13:06
• Yes, this works for $Vxx = Vyy$, but it doesn't when they have different values Mar 8 at 13:07