Deriving a PDF for an annulus distribution

Question

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

Answer 1

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

Answer 2

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

Answer 3

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.