Contour lines over SmoothDensityHistogram

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


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


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