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