Deriving a PDF for an annulus distribution

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:

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

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


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


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


13.499060

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


3.012873

• Nice! Thanks for your answer. Best, Ben – ben18785 Apr 29 '16 at 23:47
• 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}$$ – 2012rcampion Apr 30 '16 at 5:20
• An alternative representation: TransformedDistribution[(r + d) Normalize[{x, y}], {d \[Distributed] NormalDistribution[0, σ], {x, y} \[Distributed] BinormalDistribution[0]}] – J. M.'s technical difficulties May 13 '16 at 9:32

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:

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...

• Increasing the PlotPoints (say PlotPoints -> 50) will give you a smoother plot. – Bob Hanlon Apr 30 '16 at 4:14

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] *)


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


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


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.

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