11
$\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$
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$
  • $\begingroup$ Nice! Thanks for your answer. Best, Ben $\endgroup$ – ben18785 Apr 29 '16 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$ – 2012rcampion Apr 30 '16 at 5:20
  • $\begingroup$ An alternative representation: TransformedDistribution[(r + d) Normalize[{x, y}], {d \[Distributed] NormalDistribution[0, σ], {x, y} \[Distributed] BinormalDistribution[0]}] $\endgroup$ – J. M. will be back soon May 13 '16 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$
  • $\begingroup$ Increasing the PlotPoints (say PlotPoints -> 50) will give you a smoother plot. $\endgroup$ – Bob Hanlon Apr 30 '16 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$

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.