# Smooth kernel density with kernels following local function values

I can generate some data, and plot the result as follows,

n = 10000;
x = RandomReal[{0, 2 Pi}, n];
y = Sin[x] + RandomVariate[NormalDistribution[0, 0.01], {n}];
PlotRange -> {{0, 2 Pi}, {-1, 1}}]


which results in a kernel density plot like this,

In this plot you'll notice that using such a narrow kernel has meant that the kernel density estimates are jagged in the x-y directions. This is because, I presume, by default Mathematica uses kernels that do not have bivariate correlation.

What I would like is a kernel that smooths in the local directions given by the function (in this case, dictated by the sine function).

Does anyone know how to do this?

Note please do not suggest I use a wider kernel. I am not looking for smoother estimates at the expense of a higher bandwidth.

I have tried the following, but this does not appear to work as desired,

SmoothDensityHistogram[Thread[{x, y}], {Automatic, PDF[NormalDistribution[0, 0.1], Sin[#1]] &},
PlotRange -> {{0, 2 Pi}, {-1, 1}}]

• My suggestion for a wider bandwidth would have probably been better stated if I emphasized a "justifiable" bandwidth. At some point you'll need to justify that choice. A wider bandwith is not an "expense". It's about being appropriate.
– JimB
Jul 20, 2017 at 15:22

This is not an answer but an extended comment.

I think there are 3 issues for your consideration:

(1) I think there are conceptual issues with what you're proposing: not the use of a bivariate kernel with a correlation but rather that somehow one might know the functional form of the bivariate distribution and use that to construct a better nonparametric estimator. In other words, if you've got such knowledge of the bivariate form, why not use a parametric density estimator? It would seem to be a very rare case to know such information.

(2) The bandwidth you use seems way too small. Is it chosen on some theoretical basis or is it just small enough to show the jaggedness/waviness?

(3) I think the jaggedness/waviness that you see in the SmoothDensityHistogram plot is simply an artifact of not enough points of interpolation and not at all because the bivariate kernels used have zero internal correlation. There is an option for InterpolationPoints in SmoothKernelDistribution which one can increase and see the waviness disappear in the contours but I don't see such an option in SmoothDensityHistogram. (However, I could certainly have not recognized the equivalent option.)

Update

For your example bivariate density there is an "exact" solution for finding the bivariate density. (And I added a bit more spread for the y variable to make the figure a bit more interesting.)

xydist = TransformedDistribution[{x, Sin[x] + z},
{x \[Distributed] UniformDistribution[{0, 2 π}],
z \[Distributed] NormalDistribution[0, 0.1]}];


Here's a plot of a random sample along with the contours that contain 80, 90, and 95 percent of the probability.

xy = RandomVariate[xydist, 1000];
Show[ListPlot[xy, PlotStyle -> Red],
ContourPlot[PDF[xydist, {x, y}], {x, 0, 2 π}, {y, -1.4, 1.4},
PlotRange -> All,
Contours -> {0.02301837725113296, 0.09301835502825319,
0.1641455150120073, 0.27931426505406887},
ContourStyle -> Gray, ContourShading -> None],
ImageSize -> Large]


2nd update

I think I might have an example that might convince you that the very uniform jaggedness is due to SmoothDensityHistogram not evaluating the density at a sufficient number of points. I've tried changing the PlotPoints option but either I didn't try a large enough number or PlotPoints has little effect. Increasing MaxRecursion to 15 resulted in crashing my PC.

I used your exact code but increased the sample size by 100 times and the figure looks pretty much the same.

n = 1000000;
x = RandomReal[{0, 2 Pi}, n];
y = Sin[x] + RandomVariate[NormalDistribution[0, 0.01], {n}];
SmoothDensityHistogram[Thread[{x, y}], 0.01, PlotRange -> {{0, 2 Pi}, {-1, 1}}]


However, I now think you might be on to something for bivariate densities that look like a "linear" feature (such as a stream or road) as opposed to a bivariate distribution with a bump or two. Maybe looking at points within a specified radius of each point and calculating a correlation coefficient and then "applying" that correlation coefficient to the kernel would result in better bivariate density estimates. (One would still need to figure out how to justify a bandwidth.)