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

| improve this answer | |
$\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.'s discontentment 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...

| improve this answer | |
$\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.

| improve this answer | |
$\endgroup$
  • $\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 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 at 20:45

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.