# Contour lines over SmoothDensityHistogram

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?

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.

• 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. – MarcoB Jun 4 '15 at 16:55
• 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. – sguillot Jun 4 '15 at 18:24

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}]


• 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. – sguillot Jun 4 '15 at 19:28
• @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.) – JimB Jun 4 '15 at 22:25
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}]


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,