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

Improve this question
$\endgroup$

3 Answers 3

Reset to default
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

Improve this answer
$\endgroup$
3
  • $\begingroup$ Nice! Thanks for your answer. Best, Ben $\endgroup$ Commented 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$ Commented 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$ Commented 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...

Improve this answer
$\endgroup$
1
  • $\begingroup$ Increasing the PlotPoints (say PlotPoints -> 50) will give you a smoother plot. $\endgroup$ Commented 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.

Improve this answer
$\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$ Commented 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$ Commented Jan 21, 2020 at 20:45

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.