4
$\begingroup$

As an example, for a normal distribution everyone knows that 68% will fall within one standard deviation away from the center and so on.

The following code will give me a 2D Normal distribution and plot it.

a = Table[RandomVariate[NormalDistribution[0, 1]], {i, 22000}, {j, 2}];

{SmoothDensityHistogram[a, PlotLegends -> Automatic, Mesh -> 3], 
 SmoothHistogram3D[a]}

2D Normal Distribution

The question I wish to answer is "how can I calculate the probability contained above a specific mesh line?". Of course for this particular case I could add up the area under the contour curve, for example I should be able to recover the one standard deviation percentage by adding up everything under the contour curve of the SmoothHistogram3D plot inside the area $x^2 + y^2 < 1^2$, and I could find the next by adding the area $ 1^2<x^2 + y^2 < 2^2$

So now onto my actual more difficult problem, I'm plotting the Smooth Density Histogram of a double pendulum (particularly the location of the bottom pendulum).

The following code:

deqns = {Subscript[m, 1] x1''[
  t] == (\[Lambda]1[t]/Subscript[l, 1]) x1[
   t] - (\[Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]), 
  Subscript[m, 1] y1''[
  t] == (\[Lambda]1[t]/Subscript[l, 1]) y1[
   t] - (\[Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) - 
  Subscript[m, 1] g, 
  Subscript[m, 2] x2''[
  t] == (\[Lambda]2[t]/Subscript[l, 2]) (x2[t] - x1[t]), 
  Subscript[m, 2] y2''[
  t] == (\[Lambda]2[t]/Subscript[l, 2]) (y2[t] - y1[t]) - 
  Subscript[m, 2] g};

aeqns = {x1[t]^2 + y1[t]^2 == 
  Subscript[l, 1]^2, (x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 == 
  Subscript[l, 2]^2};

ics = {x1[0] == 1, y1[0] == 0, y1'[0] == 0, x2[0] == 1, y2[0] == -1, 
  y2'[0] == 0};

params = {g -> 9.81, Subscript[m, 1] -> 1, Subscript[m, 2] -> 1, 
  Subscript[l, 1] -> 1, Subscript[l, 2] -> 1};

soldp = First[
  NDSolve[{deqns, aeqns, ics} /. params, {x1, y1, x2, 
  y2, \[Lambda]1, \[Lambda]2}, {t, 0, 15000}, 
  Method -> {"IndexReduction" -> {"Pantelides", 
    "ConstraintMethod" -> "Projection"}}]];

Will give the solution to a double pendulum (where each pendulum is of length one). ($t$ is really large so if you want to run it yourself feel free to bring it down an order of magnitude).

Here is the smooth histogram,

SmoothDensityHistogram[
  Map[Function[Evaluate[{x2[#], y2[#]} /. soldp]], 
   Range[0, 15000, 0.025]], Mesh -> 5, 
  PlotRange -> {{-2, 2}, {-2, 0.1}}]

Density Probability Double Pendulum

How can I label (and calculate) a specific region between mesh lines and so that 25% (or whatever it is) of the probability is between these two mesh lines (and so forth similar to the Normal distribution example above)?

$\endgroup$
  • $\begingroup$ something like this? $\endgroup$ – kglr Feb 19 at 19:27
  • $\begingroup$ Related. $\endgroup$ – corey979 Feb 19 at 20:12
5
$\begingroup$
soldp = First[NDSolve[{deqns, aeqns, ics} /. params, {x1, y1, x2, y2, λ1, λ2}, 
   {t, 0, 15000/2}, 
   Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}]]; 
dt = Map[Function[Evaluate[{x2[#], y2[#]} /. soldp]], Range[0, 15000/2, 0.025]]; 
sdh1 = SmoothDensityHistogram[dt, Mesh -> 5, PlotRange -> {{-2, 2}, {-2, 0.1}}]

enter image description here

Using the approach from this answer to a related q/a:

We define multivariate "quantiles" based on the height of the kernel density function. The function volume[z] gives the total probability of the set of points where density exceeds z. We find the density threshold levels corresponding to desired probability coverages (this part is very slow).

As an example, each of the two regions between the contour lines corresponding to 95%, 70% and 45% probabilities have total density 25%:

skdPDF = PDF[SmoothKernelDistribution[dt]];
volume[z_?NumericQ] := Quiet@NIntegrate[skdPDF[{s, t}] Boole[skdPDF[{s, t}] >= z], 
  {s, -∞, ∞}, {t, -∞, ∞}]

{t95, t70, t45} = Quiet[FindRoot[volume[z] - # == 0, {z, 0, 1}]] & /@ {.95, .70, .45}

{{z -> 0.0705034}, {z -> 0.128013}, {z -> 0.16069}}

mesh = {{{z /. t95, Red}, {z /. t70, Green}, {z /. t45, Purple}}};
sdh2 = SmoothDensityHistogram[dt, 
  MeshFunctions -> {skdPDF[{#, #2}] &}, 
  Mesh -> mesh, MeshStyle -> Thick, PlotRange -> {{-2, 2}, {-2, 0.1}}]

enter image description here

The regions colored Yellow and Orange each have 25% probability coverage:

SmoothDensityHistogram[dt, MeshFunctions -> {skdPDF[{#, #2}] &}, 
 Mesh -> mesh, MeshStyle -> Thick, PlotRange -> {{-2, 2}, {-2, 0.1}}, 
 MeshShading -> {None, Yellow, Orange, None}]

enter image description here

SmoothHistogram3D[dt, BoxRatios -> 1, Mesh -> mesh, 
 MeshStyle -> Thick, MeshShading -> {White, Yellow, Orange, White}, 
 Lighting -> "Neutral"]

enter image description here

Show[sdh1, Graphics[Cases[Normal[sdh2], {dir__, _Line}, All]]]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks, this is exactly what I wanted. I did a search earlier and did not see the solution you referenced at the beginning. $\endgroup$ – Josh Feb 19 at 20:53
  • $\begingroup$ @Josh,my pleasure. Thank you for the accept. $\endgroup$ – kglr Feb 19 at 20:58

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.