# TransformedDistribution doesn't seems to work properly in this case

Im trying to get the distribution of the random variable defined by

$$Y_n:=\left|\sum_{k=1}^n e^{ i X_k}\right|$$

for arbitrary $$n$$, where $$X_k\sim U[0,2\pi]$$ i.i.d. Then I had written the code

iid[k_] := TransformedDistribution[Abs[Sum[Exp[2 Pi I a[i]], {i, k}]],
Table[Distributed[a[i], UniformDistribution], {i, k}]]


However mathematica doesn't compute it, I mean, it just repeat it definition when executed. I tried to plot it probability distribution writing

Plot[PDF[iid[5],x],{x,0,5}]


but mathematica show me an empty graph. I cannot compute either PDF[iid[2],0.5], so the above code that defines the distribution seems not working. Then I tried this definition

iid[k_] := TransformedDistribution[Norm[Sum[{Cos[a[i]], Sin[a[i]]}, {i, k}]],
Table[Distributed[a[i], UniformDistribution], {i, k}]]


However it neither worked. There is a way to get the distribution of $$Y_n$$?

• Do you know the distribution of $Y_n$ but just want to have it found using commands in Mathematica? – JimB Jul 18 at 4:30
• @JimB I dont know the exact distribution of $Y_n$, just an approximation. Anyway I want to know how to handle these cases in mathematica, assuming that I dont had an approximation – Masacroso Jul 18 at 4:32
• The exact density for $n=2$ is $\frac{2}{\pi \sqrt{4-y^2}}$ which can be determined with Mathematica commands. If I find time tomorrow, I'll add in the approach to determine that and I'll think more about the general case. – JimB Jul 18 at 4:38
• What does the absolute value add?? Or stated differently, what component is negative? And how do you propose to plot a complex random variable (assuming your i denotes the imaginary unit)? – wolfies Jul 18 at 16:33
• @wolfies all your questions are unrelated to what Im asking. I indeed plotted and studied the distribution of $Y_n$ using Julia and R, however I want to see what I could do using Wolfram mathematica and the function TransformedDistribution. Maybe what Im trying to accomplish is too heavy for symbolic computation. – Masacroso Jul 18 at 17:01

This is only an extended comment in that it is a brute force approach for $$n=2$$.

y2 = (Cos[x1] + Cos[x2])^2 + (Sin[x1] + Sin[x2])^2;

(* Suppose we want to see the region where y^2 < 1.
That region consists of the two black polygons below. *)
ContourPlot[y2, {x1, 0, 2 π}, {x2, 0, 2 π}, Contours -> {1}, AspectRatio -> 1,
ColorFunction -> GrayLevel, FrameLabel -> (Style[#, Bold, Italic,
18] &) /@ {"\!$$\*SubscriptBox[\(x$$, $$1$$]\)",
"\!$$\*SubscriptBox[\(x$$, $$2$$]\)"},


Because x1 and x2 have uniform distributions, the probability of y < 1 is just the area of the black polygons divided by the total area.

(* The x2 coordinates where x1 = 0 are found as follows *)
sol = FullSimplify[Solve[{y^2 == (Cos[x1] + Cos[x2])^2 + (Sin[x1] + Sin[x2])^2,
0 < x1 < 2 π} /. x2 -> 0, x1], Assumptions -> 0 < y < 2]
(* {{x1\[Rule](3 π)/2+ArcSin[1/2 (-2+y^2)]},{x1\[Rule]ArcCos[-1+y^2/2]}} *)
{xLow, xHigh} = {x1 /. sol[[2]], x1 /. sol[[1]]};

(* So the lower polygon has coordinates *)
poly = {{xLow, 0}, {xHigh, 0}, {2 π, xLow}, {2 π, xHigh}, {xLow, 0}};

(* As a check look at the polygon and the y = 1 contours *)
Show[ContourPlot[(Cos[x1] + Cos[x2])^2 + (Sin[x1] + Sin[x2])^2, {x1, 0, 2 π}, {x2, 0, 2 π},
Contours -> {1}, AspectRatio -> 1, ColorFunction -> GrayLevel,
FrameLabel -> (Style[#, Bold, Italic, 18] &) /@ {"\!$$\*SubscriptBox[\(x$$, $$1$$]\)",
"\!$$\*SubscriptBox[\(x$$, $$2$$]\)"},
ListPlot[poly /. y -> 1, Joined -> True, PlotStyle -> {Thickness[0.01], Red}]]


So the probability that Y < y is twice the area of the lower polygon divided by the total area (4 π^2)

area = FullSimplify[Sum[poly[[i - 1, 1]]*poly[[i, 2]] - poly[[i - 1, 2]]*poly[[i, 1]],
{i, 2, Length[poly]}]]/(4 π^2)
(* ArcCos[1-y^2/2]/π *)


Finally the pdf is the derivative of the area with respect to y:

pdf = FullSimplify[D[area, y], Assumptions -> 0 < y < 2]
(* 2/(π Sqrt[4-y^2]) *)
Plot[pdf, {y, 0, 2}, PlotRange -> {{0, 2}, {0, 2}}, Frame -> True,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"\!$$\*SubscriptBox[\(Y$$, $$2$$]\)",
"Probability density"}]


One more check: random samples:

SeedRandom[12345];
x1 = RandomVariate[UniformDistribution[{0, 2 π}], 100000];
x2 = RandomVariate[UniformDistribution[{0, 2 π}], 100000];
y2 = Sqrt[(Cos[x1] + Cos[x2])^2 + (Sin[x1] + Sin[x2])^2];
Show[Histogram[y2, 100, "PDF", Frame -> True,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"\!$$\*SubscriptBox[\(Y$$, $$2$$]\)",
"Probability density"}],
Plot[pdf, {y, 0, 2}, PlotRange -> {{0, 2}, {0, 5}}]]


Addition: $$n=3$$

Here is an estimate of the distribution of $$Y_3$$ based on simulations:

n = 10000000;
SeedRandom[12345];
x = RandomVariate[UniformDistribution[{0, 2 π}], {n, 3}];
y3 = Sqrt[(Cos[#[[1]]] + Cos[#[[2]]] + Cos[#[[3]]])^2 + (Sin[#[[1]]] +
Sin[#[[2]]] + Sin[#[[3]]])^2] & /@ x;
skd = SmoothKernelDistribution[y3,
0.015, {"Bounded", {0, 3}, "Gaussian"}];
Show[Histogram[y3, 100, "PDF", PlotLabel -> Style["n = 3", Bold, 18]],
Plot[PDF[skd, x], {x, 0, 3}, PlotRange -> {{0, 3}, {0, 1}}]]


• thank you for your time and effort, I will take a carefully look to the code – Masacroso Jul 19 at 23:51
• Note that the distribution looks wildly different for $n=3$. While there might be a nice equation for when $n>2$, I doubt it (but I've been wrong many times before). Do you really need compact formulas or will nonparametric density estimates based on simulations work for you? – JimB Jul 19 at 23:54
• no, I dont need compact formulas. I was just trying to test the capabilities of mathematica to calculate exact densities. I didn't realized that this could be too much for this case – Masacroso Jul 20 at 0:02
• Do you have in mind a software package that can do it? – JimB Jul 20 at 0:10
• no, I dont know any software that can do it – Masacroso Jul 20 at 0:49