8
$\begingroup$

I am using SmoothDensityHistogram on a data set of the form {{x1, y1}, {x2, y2}, …, {x_n, y_n}}, and I would like to also show the contour lines that enclose 68%, 95% and 99% of the points.

With the option MeshFunctions -> {#3 &}, Mesh -> 3 I can have 3 contour lines, but how can I set the probability at which the contours lines are?

enter image description here

As this image show, the distribution of points does not necessarily follow a binormal distribution, so I need something more general than confidence ellipses calculated with Mean and Covariance.

It seems like a common enough plot that an easy solution should exist but I can't figure out how.

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ I think you are trying to draw the confidence ellipse corresponding to those points. If that's the case, take a look at this question: How to draw confidence ellipse from a covariance matrix?. You can easily calculate statistical descriptors of your dataset using Mean and Covariance. $\endgroup$
    – MarcoB
    Commented Jun 4, 2015 at 16:55
  • 1
    $\begingroup$ Thanks @MarcoB. But my distribution is not binomial, and the confidence region cannot be described by an ellipse. I've edited my question to clarify it. $\endgroup$
    – sguillot
    Commented Jun 4, 2015 at 18:24

3 Answers 3

Reset to default
8
$\begingroup$

Here is a brute-force method (and I'm sure there are many more efficient approaches):

data = RandomVariate[BinormalDistribution[.75], 100];

(* Calculate a nonparametric density estimate *)
d = SmoothKernelDistribution[data];

(* Evaluate the estimated density function over a grid of points and sort by the density values from high to low *)
pdf = Reverse[
   Sort[Flatten[
     Table[PDF[d, {x, y}], {x, -3, 3, 0.05}, {y, -3, 3, 0.05}]]]];

(* Create a table of cumulative pdf values that correspond to the volumes of interest *)
cdf = Accumulate[pdf]/Total[pdf];
contours = 
 pdf[[Flatten[
    Table[FirstPosition[cdf, 
      p_ /; p >= alpha], {alpha, {0.68, 0.95, 0.99}}]]]];

(* Create a series of figures and then overlay them all *)
sdh = SmoothDensityHistogram[data];
cp = ContourPlot[PDF[d, {x, y}], {x, -3, 3}, {y, -3, 3}, 
   Contours -> contours, ContourShading -> None];
lp = ListPlot[data];
Show[{sdh, cp, lp}]

nonparametric density contours

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ Thanks for the solution Jim Baldwin. It works great! I would just add that a "PlotRange -> All" might be needed for the ContourPlot because depending on the contour values, some might not show. $\endgroup$
    – sguillot
    Commented Jun 4, 2015 at 19:28
  • $\begingroup$ @sguillot Agreed. Thanks. A more efficient method is needed (and maybe more accurate as I didn't use interpolation to find the contour values - just the first instance where the "cumulative" sum first equaled or exceeded the desired coverage probability) because of the my need to model animal home ranges with nonparametric density estimates. One possibility is to call the adehabitat package in R from Mathematica. (Determining the area of the associated contours is also of interest.) $\endgroup$
    – JimB
    Commented Jun 4, 2015 at 22:25
  • $\begingroup$ I find your answer very helpful, I want to ask if you can calculate the average of the point inside a certain contour. kind regards $\endgroup$
    – Renato QE
    Commented Oct 10, 2023 at 4:27
6
$\begingroup$

Here's a solution like Jim Baldwin's, but a little less brutal. I don't see a need for the mesh when you can get the density estimate for each data point:

(* two dimensional data *)
data = RandomVariate[BinormalDistribution[.5], 200];

(* nonparametric density estimate *)
d = SmoothKernelDistribution[data];

(* logged density estimate at data *)
p = PDF[d, data];

(* quantile for countour height *)
q = 1 - {0.68, 0.95, 0.99};
c = Quantile[p, q];

Show[
    DensityPlot[PDF[d, {x, y}], ##],
    ContourPlot[PDF[d, {x, y}], ##, Contours -> c, 
    ContourShading -> None],
    ListPlot[data, PlotStyle -> {Red, PointSize[0.01]}]
]&[{x, -4, 4}, {y, -4, 4}, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
7
  • $\begingroup$ Your solution is much faster but finds contours that contain a specified proportion of data points rather than a specified volume under the nonparametric density surface. I think both solutions have use. $\endgroup$
    – JimB
    Commented Jun 5, 2015 at 5:00
  • $\begingroup$ @Jim Oh, absolutely! Thanks for pointing out the difference. OP asked about points, but highest probability density regions are probably more widely used. $\endgroup$
    – Ian
    Commented Jun 5, 2015 at 11:48
  • $\begingroup$ Thanks a lot @Ian. I will try also you solution which might be more efficient, especially since I will need to do many of these contours plots. Just as a side note, I'm guessing that your solution and Jim Baldwin's should converge to the same regions in the limit of large number of points (and for Gaussian distributions), no? $\endgroup$
    – sguillot
    Commented Jun 6, 2015 at 21:47
  • $\begingroup$ I would guess the same. $\endgroup$
    – Ian
    Commented Jun 7, 2015 at 2:52
  • $\begingroup$ @Ian The quickness to which the two approaches match in terms of the volume (probability) and area covered (both amount of area and what area actually in the contour) also depends on the proportion of points enclosed. As you might imagine, there will typically be less agreement for 95% vs. 50%. On a separate issue a logical argument might be that as one has gone to the trouble of estimating a bivariate probability distribution, why not choose contours that are associated with the enclosed probability rather than the enclosed number of points. $\endgroup$
    – JimB
    Commented Jun 9, 2015 at 15:08
4
$\begingroup$ 
data = RandomVariate[BinormalDistribution[.5], 200];
skdPDF = PDF[SmoothKernelDistribution[data]];

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:

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

Find the density threshold levels corresponding to desired probability coverages (this is very slow):

{t99, t95, t68} = Quiet[FindRoot[volume[z] - # == 0, {z, 0, 1}]]& /@ {.99, .95, .68}

{{z -> 0.002508}, {z -> 0.008514}, {z -> 0.045498}}

SmoothDensityHistogram[data, MeshFunctions -> {skdPDF[{#, #2}] &}, 
  Mesh -> {{{z /. t99, Red}, {z /. t95, Green}, {z /. t68, Purple}}},
  MeshStyle -> Thick, PlotRange -> {{-4, 4}, {-4, 4}},
  Epilog -> {Black, PointSize[Medium], Point @ data}]

enter image description here

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.