1
$\begingroup$

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$?

$\endgroup$
  • $\begingroup$ Do you know the distribution of $Y_n$ but just want to have it found using commands in Mathematica? $\endgroup$ – JimB Jul 18 at 4:30
  • $\begingroup$ @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 $\endgroup$ – Masacroso Jul 18 at 4:32
  • $\begingroup$ 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. $\endgroup$ – JimB Jul 18 at 4:38
  • $\begingroup$ 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)? $\endgroup$ – wolfies Jul 18 at 16:33
  • $\begingroup$ @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. $\endgroup$ – Masacroso Jul 18 at 17:01
2
$\begingroup$

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\)]\)"},
 PlotRangePadding -> None, PlotRangeClipping -> False, ImagePadding -> 50]

Contour plot where y equals 1

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\)]\)"},
  PlotRangePadding -> None, PlotRangeClipping -> False, ImagePadding -> 50],
 ListPlot[poly /. y -> 1, Joined -> True, PlotStyle -> {Thickness[0.01], Red}]]

Contour plot with lower polygon outlined in 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"}]

pdf of y_2

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

Histogram and probability density

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

Density function for n = 3

$\endgroup$
  • $\begingroup$ thank you for your time and effort, I will take a carefully look to the code $\endgroup$ – Masacroso Jul 19 at 23:51
  • $\begingroup$ 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? $\endgroup$ – JimB Jul 19 at 23:54
  • $\begingroup$ 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 $\endgroup$ – Masacroso Jul 20 at 0:02
  • $\begingroup$ Do you have in mind a software package that can do it? $\endgroup$ – JimB Jul 20 at 0:10
  • $\begingroup$ no, I dont know any software that can do it $\endgroup$ – Masacroso Jul 20 at 0:49

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.