12
$\begingroup$

I would like to derive the PDF for an ''annulus'' type distribution, defined by the parameters $\theta\sim U(0,2 \pi)$ and $d\sim N(0,\sigma)$; where $\theta$ is the angle round a circle of radius $r$, and $d$ is the deviation normal to the circle's perimeter.

I can generate sample points easy enough, by the following:

fTorusProjector[xCentre_, yCentre_, θ_,
     d_] := {xCentre + d Sin[θ], yCentre + d Cos[θ]}

fTorusRand[r_, σ_] := 
     Module[{θ = RandomReal[{0, 2 Pi}], 
   d = RandomVariate[NormalDistribution[0, σ], {1}][[1]]}, 
   fTorusProjector[r Sin[θ], r Cos[θ], θ, d]]

ListPlot[Table[fTorusRand[5, 0.5], {i, 1, 10000}]]

Which results in a plot similar to the one below: enter image description here

My question is, is there a way to derive an analytic form for the PDF of this distribution? I can't help thinking that there's probably a way I'm missing, but am not sure how to proceed here.

Even an approximate PDF / one that mimics the above behaviour would suffice.

Edit: I have tried using the following to get Mathematica to approximate a PDF of the above sample:

 empD = EmpiricalDistribution[data];

However, when I try to draw this using:

ContourPlot[PDF[empD, {x, y}], {x, -7, 7}, {y, -7, 7}]

I just get a blank plot.

Edit 2: Since $\theta$ and $d$ are independent, I can derive the joint PDF in these coordinates:

$p(\theta,d) = \frac{1}{2\pi} \times \frac{1}{2\pi} \exp(-d^2/2\sigma^2)$

I suppose I can then use Jacobians to transform back into the $(x,y)$ frame, although am not sure how to do this?

Best,

Ben

$\endgroup$

3 Answers 3

10
$\begingroup$

An alternative, and faster, way to generate samples with the desired distribution using TransformedDistribution:

ClearAll[distF]
distF[r_, σ_] := TransformedDistribution[{(r + d ) Sin[θ], (r + d ) Cos[θ]}, 
  {Distributed[d, NormalDistribution[0, σ]], Distributed[θ,
                UniformDistribution[{0, 2 Pi}]]}]
ListPlot[RandomVariate[distF[5, .5] , 10000], PlotStyle -> PointSize[Tiny]]

Mathematica graphics

data = RandomVariate[distF[5, .5] , 1000000];
empdist = SmoothKernelDistribution[data];
ContourPlot[PDF[empdist, {x, y}], {x, -8, 8}, {y, -8, 8}, PlotPoints -> 50]

Mathematica graphics

Table[fTorusRand[5, 0.5], {i, 1, 1000000}]; // AbsoluteTiming // First

13.499060

RandomVariate[distF[5, .5] , 1000000]; // AbsoluteTiming // First

3.012873

$\endgroup$
3
  • $\begingroup$ Nice! Thanks for your answer. Best, Ben $\endgroup$
    – ben18785
    Apr 29, 2016 at 23:47
  • 7
    $\begingroup$ We can actually find a closed form for the PDF: $$\frac{\exp\left\{-\frac{1}{2}\left(\frac{r_0-r} {\sigma}\right)^2\right\} + \exp\left\{-\frac{1}{2}\left(\frac{r_0+r} {\sigma}\right)^2\right\}}{\left(2\pi\right)^{3/2}\sigma r}$$ $\endgroup$ Apr 30, 2016 at 5:20
  • $\begingroup$ An alternative representation: TransformedDistribution[(r + d) Normalize[{x, y}], {d \[Distributed] NormalDistribution[0, σ], {x, y} \[Distributed] BinormalDistribution[0]}] $\endgroup$ May 13, 2016 at 9:32
7
$\begingroup$

Sorry for the trouble - I have found a way:

data = Table[fTorusRand[5, 0.5], {i, 1, 100000}];
empD1 = SmoothKernelDistribution[data];
ContourPlot[PDF[empD1, {x, y}], {x, -8, 8}, {y, -8, 8}]

Gives me what I want:

enter image description here

Edit: So using the Jacobian here, I can get an exact PDF. Still messing around with the algebra though: $x = (r+d) \sin(\theta)$ and $y = (r+d) \cos(\theta)$. So I can implicitly differentiate these to find the Jacobians...

$\endgroup$
1
  • $\begingroup$ Increasing the PlotPoints (say PlotPoints -> 50) will give you a smoother plot. $\endgroup$
    – Bob Hanlon
    Apr 30, 2016 at 4:14
6
$\begingroup$

This is a naive approach, so I wonder if I am overlooking some complication here; I would have thought that, given the description of your distribution, you could consider the "solid of revolution" generated by rotating the PDF of a normal distribution with the required parameters around the vertical $z$ axis, so something like the following:

PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]

(* Out: E^(-(1/2) (-6+Sqrt[x^2+y^2])^2)/Sqrt[2 Pi] *)

PDF

Plot3D[
  Evaluate[PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]],
  {x, -10, 10}, {y, -10, 10}, PlotPoints -> 75
]

3D plot

ContourPlot[
  Evaluate[PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]],
  {x, -10, 10}, {y, -10, 10}, PlotPoints -> 75
]

contour plot

Of course, as it is the expression is not properly normalized, so an appropriate factor should be added to make sure that its integral is equal to 1.

$\endgroup$
2
  • $\begingroup$ How can we calculate in Mathematica the volume below this function and the XY plane? Could you mention the integration code line? $\endgroup$
    – Gouz
    Jan 21, 2020 at 0:27
  • 1
    $\begingroup$ @Gouz NIntegrate[Evaluate[PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}] would do that. $\endgroup$
    – MarcoB
    Jan 21, 2020 at 20:45

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.