5
$\begingroup$

I want to use SmoothHistogram-like plot using smoothing kernel $f(x)$ which is uniform $(a,b)$ with $a,b$ adjusted for each $x$ to include 11 points on the left, and 3 points on the right of $x$.

This kind of approach is called "Nearest Neighbor" method by Silverman.

Can anyone see an efficient way to implement this, or existing code that does this?

Here's an example using SmoothHistogram, but I want density estimation to use uniform kernel, to get visualization style from searke's answer. (However, his approach doesn't scale to 10k points)

Clear[x];
dist[k_] := 
  TransformedDistribution[x^(k/2), 
   x \[Distributed] UniformDistribution[{0, 1}]];
datasets = Table[RandomVariate[#, 10000] &@dist[k], {k, 2, 10, 2}];
SmoothHistogram[datasets, PlotLegends -> Range[1, 10]]

enter image description here

$\endgroup$
9
  • 1
    $\begingroup$ What do you mean by "unbalanced data" ? Do you mean "not having a symmetric density"? Skewed to the right? Bounded only on the left or right but not both? Also, @seake presented an interesting idea that had a reasonable visual match to the underlying density for a single dataset. That does not show superiority of that method (especially because there is a section with positive density where the density is zero). $\endgroup$
    – JimB
    Nov 25, 2023 at 1:14
  • 2
    $\begingroup$ Unbalanced data means the kind of data where where the variance of density estimate from a fixed-width kernel estimator changes significantly over the support of the density. So I want to adjust the bandwidth dynamically to include roughly the same number of observations at different points of the support $\endgroup$ Nov 25, 2023 at 2:11
  • $\begingroup$ Don't you think that happens most of the time? Consider a "nice" normal distribution: SeedRandom[12345]; x = RandomVariate[NormalDistribution[0, 1], {1000, 100}]; {StandardDeviation[PDF[SmoothKernelDistribution[#], 0] & /@ x], StandardDeviation[PDF[SmoothKernelDistribution[#], 1] & /@ x], StandardDeviation[PDF[SmoothKernelDistribution[#], 2] & /@ x], StandardDeviation[PDF[SmoothKernelDistribution[#], 3] & /@ x]}. The result is {0.0456611, 0.0365265, 0.0218769, 0.00660253}. $\endgroup$
    – JimB
    Nov 25, 2023 at 2:59
  • 1
    $\begingroup$ I think you're trying to reinvent the wheel. That's called a "nearest neighbor" approach. Such a question is better suited for CrossValidated and then back here for implementation. I say that because at some point you'll need some statistical theory to back up that approach and it appears that the current assessment of the goodness is just visual. $\endgroup$
    – JimB
    Nov 25, 2023 at 5:05
  • 2
    $\begingroup$ See Silverman's book on the section of the nearest neighbor approach: ned.ipac.caltech.edu/level5/March02/Silverman/paper.pdf. $\endgroup$
    – JimB
    Nov 25, 2023 at 5:09

2 Answers 2

6
+150
$\begingroup$

Here I'm first showing the (standard) nearest neighbor approach along with the generalized nearest neighbor approach (which is related to the nearest neighbor approach but uses kernels). The asymmetric generalized nearest neigbors approach you propose will follow when I have time although I believe the results will be so similar that the faster (standard) nearest neighbor approach should be more than adequate and easily reproducible in just about any programming language.

(* Generate some data *)
SeedRandom[12345];
dist = TransformedDistribution[x^2, x \[Distributed] UniformDistribution[{0, 1}]];
n = 10000;
data = Rationalize[RandomVariate[dist, n] // Sort, 0];

(* Set values to evaluate density estimates *)
m = 1000;
t = Range[Min[data], Max[data], (Max[data] - Min[data])/(m - 1)];

(Standard) Nearest neighbor approach

(* https://ned.ipac.caltech.edu/level5/March02/Silverman/paper.pdf *)
(* k-th nearest neighbor *)
k = Sqrt[n] // Ceiling;

(* Distance to k-th nearest neigbor *)    
d = Max[Abs[Nearest[data, #, k] - #]] & /@ t;

(* Estimated density *)
pdf = (k - 1)/(2 n d);
(* Show results *)
Show[ListPlot[Transpose[{t, pdf}], Joined -> True, PlotRange -> {{0, 1}, {0, 8}}],
 Plot[PDF[dist, x], {x, 0, 1}, PlotStyle -> Red, PlotRange -> {{0, 1}, {0, 8}}]]

Standard nearest neighbor approach with true density

Generalized nearest neighbor method using a retangular kernel

(* https://ned.ipac.caltech.edu/level5/March02/Silverman/paper.pdf *)
(* k-th nearest neighbor *)
k = Sqrt[n] // Ceiling;

(* Distance to k-th nearest neighbor *)
d = Max[Abs[Nearest[data, #, k] - #]] & /@ t;

(* Rectangular kernel *)
kernel[x_] := Piecewise[{{1/2, -1 < x < 1}}]

(* Estimated density *)
pdfg = ((1/(n d[[#]])) Sum[kernel[(t[[#]] - data[[j]])/d[[#]]], {j, n}]) & /@ Range[m];

(* Show results *)
Show[ListPlot[Transpose[{t, pdfg}], Joined -> True, PlotRange -> {{0, 1}, {0, 8}}],
 Plot[PDF[dist, x], {x, 0, 1}, PlotStyle -> Red, PlotRange -> {{0, 1}, {0, 8}}]]

Generalized nearest neighbor method with true density

This looks very similar to the simpler nearest neighbor method because with the rectangular kernel the estimates are identical (for all evaluations within the range of the data - not so for outside the range of the data). However, they may not appear exactly identical when machine precision numbers are used. (Hence, that is why I used Rationalize above. In practice, one probably just needs machine precision.)

Asymmetric nearest neighbors

Here is an approach based on your desire to mimic the (standard) nearest neighbor approach but with 11 neighbors to the left and 3 neighbors to the right of an evaluation point. (Not that there isn't any but I'm unaware of any justification for doing so.)

kLower = 11; (* Number of data points less than an evaluation point *) 
kUpper =  3; (* Number of data points greater than an evaluation point *)

(* Predict data point index from value of data point *)
f = Interpolation[Transpose[{data, Range[n]}], InterpolationOrder -> 1];

(* Construct "asymmetrical neighbors" density estimate *)
pdfan = If[f[#] >= kLower && f[#] < n - kUpper,
  (kLower + kUpper)/(n (data[[Ceiling[f[#]] + kUpper]] - 
          data[[Floor[f[#]] - kLower + 1]]))] & /@ t;

(* Plot results and true density *)
Show[ListPlot[Transpose[{t, pdfan}], Joined -> True, PlotRange -> {{0, 1}, {0, 8}}],
  Plot[PDF[dist, x], {x, 0, 1}, PlotStyle -> Red, PlotRange -> {{0, 1}, {0, 8}}]]

Asymmetric nearest neighbor density estimate and true density

When the sample size is large (as is in this example) choosing larger values for kLower and kUpper will result in less variability in the density estimates. However, with all nearest neighbor estimators, one will not (never?) get a smooth estimate of density.

Asymmetric generalized nearest neighbors

Not programmed yet but with a rectangular (i.e., uniform) kernel, I suspect that this will result in an estimate nearly identical if not identical (within the range of the data) although much slower.

$\endgroup$
1
  • $\begingroup$ I thought the "skewed" approach is needed to able to visualize log-log density better, which focuses disproportionately on left hand side. But turns out symmetric estimation works as well....and using skewed version is biased -- wolframcloud.com/obj/yaroslavvb/newton/… $\endgroup$ Feb 6 at 23:52
7
$\begingroup$

Why not use the "Adaptive" option for SmoothKernelDistribution? This is in the same spirit as having a kernel density estimate based on the nearest set of data points such that areas of higher density have smaller bandwidths and areas of lower density have larger bandwidths.

SeedRandom[12345]
vals = RandomVariate[#, 100] &@MarchenkoPasturDistribution[1];
skdFixed = SmoothKernelDistribution[vals, Automatic, "Gaussian"];
skdAdaptive = SmoothKernelDistribution[vals, {"Adaptive", Automatic, 0.5},
  {"Bounded", {0, ∞}, "Gaussian"}];
Plot[{PDF[MarchenkoPasturDistribution[1], x], PDF[skdFixed, x], 
  PDF[skdAdaptive, x]}, {x, 0, 5}, PlotRange -> {All, {0, 1}},
 PlotLegends -> {"True density", "Fixed width kernel", "Adaptive kernel"},
 PlotStyle -> {Blue, {Orange, Thickness[0.02]}, Green}]

True density and fixed width and adaptive width kernel estimation

If you're needing this to work with 10,000+ data points, consider 100,000 data points:

SeedRandom[12345];
vals = RandomVariate[#, 100000] &@MarchenkoPasturDistribution[1];
skdFixed = SmoothKernelDistribution[vals, Automatic, "Gaussian"];

AbsoluteTiming[skdAdaptive = SmoothKernelDistribution[vals, 
  {"Adaptive", Automatic, 0.5}, {"Bounded", {0, ∞}, "Gaussian"}];]
(* {0.120044, Null} *)

Plot[{PDF[MarchenkoPasturDistribution[1], x], PDF[skdFixed, x], 
  PDF[skdAdaptive, x]}, {x, 0, 5}, PlotRange -> {All, {0, 1.5}},
 PlotLegends -> {"True density", "Fixed width kernel", "Adaptive kernel"},
 PlotStyle -> {Blue, {Orange, Thickness[0.02]}, Green}]

Density estimates with 100000 data points

$\endgroup$
9
  • $\begingroup$ Well, from the diagrams it's clear Adaptive setting works. However, I'm looking for uniform kernel $\endgroup$ Nov 27, 2023 at 22:58
  • 2
    $\begingroup$ Understood. As is a habit of statisticians, we many times suggest alternatives to the specific question asked. And this forum does at times get such "question/answer deviations" from experts in other fields. I think that's a strength of this forum. (And as I know I didn't answer your question on purpose, I don't expect to get any prizes.) $\endgroup$
    – JimB
    Nov 27, 2023 at 23:23
  • $\begingroup$ Also, Equation 2.4 in Silverman's book does give a "kernel" approach to the nearest neighbor approach. $\endgroup$
    – JimB
    Nov 27, 2023 at 23:26
  • 1
    $\begingroup$ To state the obvious: there is a difference between explaining and justifying. I'm not seeing much difference in the "justifying" between choosing "a sensitivity parameter between 0 and 1" and "11 observations to the left and 3 observations on the right." $\endgroup$
    – JimB
    Nov 28, 2023 at 3:11
  • 1
    $\begingroup$ If that was a requirement, much of the wonderful things that Mathematica provides would not be publishable. Both the "Adaptive" and "Bounded" aspects are even easily programmed with "basic" code even in Excel. But the nearest neighbor approach (which is a fine approach) not so easily or efficiently programmed. (But I certainly could be wrong with that last statement.) If I have time this week, I'll produce the "basic" Mathematica and R code for a bounded and adaptive density estimator. $\endgroup$
    – JimB
    Nov 28, 2023 at 16:25

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.