14
$\begingroup$

As far as I understand, kernel density estimation does not make any assumptions on the moments of the underlying density, and just requires smoothness. The Cauchy density function is quite smooth. Even still, when I try to do KDE using density() in R for random draws from Cauchy distribution, I get incredibly inaccurate answers:

set.seed(1)
foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(density(rt(1e3, df = 1)), col = "red")

produces this plot: enter image description here

Repeat the above with different seeds or increasing the sample size can give further erratic estimates. The default kernel is Gaussian in R. Changing the kernel to any of the other options doesn't improve the output.

Question: What assumptions does Cauchy violate for KDEs? If it doesn't, then why do we see KDEs failing so miserably here?


Edit: @cdalitz has identified that the problem is where the kde is evaluating the density. The default is 3*bw*range(x), which for Cauchy can be quite large. Which means, by default density tries to estimate the KDE at 512 points sparsely distributed on the x-axis.

To test this, I change the from and to in the density estimation and see that if I run density twice with two sets of evaluating points, so the densities match:

set.seed(1)
samp <- rt(1e4, df = 1)

bd <- 10
den1 <- density(samp, from=-bd, to=bd, n=512)
den2 <- density(samp, from =-2*bd, to = 2*bd,  n =512)

foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(den2, col = "blue", type = "b")
lines(den1, col = "red", type = "b")

This produces the estimates below: enter image description here

The quality here is much better than before. However, now if instead of 2*bd, I change this to 50*bd, I get that the density estimate even around 0 is very different!

set.seed(1)
samp <- rt(1e4, df = 1)

bd <- 10
den1 <- density(samp, from=-bd, to=bd, n=512)
den2 <- density(samp, from =-50*bd, to = 50*bd,  n =512)

foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l', ylim = c(0,.7))
lines(den2, col = "blue", type = "b")
lines(den1, col = "red", type = "b")

enter image description here

How does evaluating the density at sparse points change the density evaluation process around $x = 0$ (the bandwidth chosen is the same for both den1 and den2)? The KD estimate at any point $x$ is $$ \hat{f}(x) = \dfrac{1}{nh} \sum_{t=1}^{n} K\left( \dfrac{x - x_i}{h}\right)\,. $$

The density estimate shouldn't change at a given value of $x = a_1$ if the density is also being evaluated at other points. What am I missing here?

$\endgroup$
15
  • 2
    $\begingroup$ The default bandwidth bw.nrd uses 1.06 times the minimum of the standard deviation and the interquartile range divided by 1.34 times the sample size to the negative one-fifth power. This quantity thus converges even in the case of a distribution with no moments. $\endgroup$
    – Xi'an
    Feb 2, 2022 at 7:19
  • 1
    $\begingroup$ This 1992 Biometrika paper by Chiu examines automatic selections of the prior for arbitrary densities, including the Cauchy as a test case, and with bandwidths that end up being close to those for the Normal test case. $\endgroup$
    – Xi'an
    Feb 2, 2022 at 8:21
  • 1
    $\begingroup$ @Greenparker: since $\hat\sigma_n$ includes the interquartile range, it does not diverge. $\endgroup$
    – Xi'an
    Feb 2, 2022 at 8:39
  • 2
    $\begingroup$ To avoid further search in the wrong direction: The default bandwidth selection works well in this case! Beware however that density chooses sample points in the range of the data at which the estimator is evaluated. Your observed problem is just an artifact created by this sample point selection. $\endgroup$
    – cdalitz
    Feb 2, 2022 at 9:05
  • 1
    $\begingroup$ No, the error is already there before the call to approx. You can verify this by copying the code of density.default and set a breakpoint (with browser()) before the call to approx. Note that you must replaxe C_BinDist with stats:::C_BinDist, because otherwise it is not found. $\endgroup$
    – cdalitz
    Feb 3, 2022 at 9:09

3 Answers 3

9
$\begingroup$

The problem lies in the sampling of the x-range by density. By default, density samples 512 values in the range of the data. As the Cauchy distribution has a heavy tail, you will only achieve a very sparse sampling of the density function with the default settings.

If we sample 512 values between -20 and +20, the density approximation is decent, even with the default bandwidth selection rule:

set.seed(1)
foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(density(rt(1e3, df = 1), from=-20, to=20, n=512), col = "red")

enter image description here

Note also that from, to and n only affect the sample points at which the kernel density estimator is evaluated and have no effect on the bandwidth:

> set.seed(1)
> x <- rt(1e3, df = 1)
> density(x, from=-20, to=20, n=512)$bw
[1] 0.351773
> density(x)$bw
[1] 0.351773
$\endgroup$
6
  • $\begingroup$ Thanks! But there seems to be something going on in the density function. Consider running the below code where density is evaluated at three points: at 0, -a, +a, where a may change. The evaluated density at 0 should be the same, but it isn't: set.seed(1) samp <- rt(1e3, df = 1) bd <- 10 dum1 <- density(samp, from=-bd, to=bd, n=3) dum2 <- density(samp, from =-500, to = 500, n =3) foo <- seq(-50, 50, length = 1e3) plot(foo, dt(foo, df = 1), type = 'l', ylim = c(0,.50)) lines(dum2, col = "blue", type = "b") lines(dum1, col = "red", type = "b") $\endgroup$ Feb 2, 2022 at 11:00
  • 1
    $\begingroup$ The code typeset is horrible, so I'll edit the question. $\endgroup$ Feb 2, 2022 at 11:03
  • 1
    $\begingroup$ @Greenparker Please read the documentation of the parameter n: It is rounded to a power of 2 (for using FFT), and later it is interpolated. This explains slight differences in the values due to interpolation (3 is not a power of 2). I wonder, however, why this is necessary since there is an alternative algorithm (FFTW) instead of FFT that does not have this restriction. $\endgroup$
    – cdalitz
    Feb 2, 2022 at 11:08
  • 1
    $\begingroup$ The rounding to the power of 2 occurs only when n > 512 (or atleast that's what the documentation says. I've also edited the question, to highlight my larger question. Thanks a ton for bringing some insight, but I still don't quite see how sparse evaluations will change the density estimate at a particular point (without interpolation -- that's why I now plot the points) $\endgroup$ Feb 2, 2022 at 11:29
  • 1
    $\begingroup$ @Greenparker You can view the source code of density by entering desity.default (without parentheses). Then you will see that n is automatically set to a minimum value of 512, even when a smaller value is provided by the user. This is later approximately reversed by interpolating at the actual values. $\endgroup$
    – cdalitz
    Feb 3, 2022 at 8:57
7
$\begingroup$

In R, the procedure density gives a kernel density estimator (KDE) of data.

In order to get a suitable histogram of the very heavy tailed Cauchy distribution, it is usually necessary to disregard more than a few values in the tails of the sample.

Then in order show the actual Cauchy density curve (black below), it is necessary to adjust the density by dividing by the proportion of the sample plotted in the histogram.

Also, to get a the best KDE (dotted red) you may need to use parameter adj tp adjust the default bandwidth. (I used the default.) Of course, the histogram would be 'smoother' if it had fewer bins; the KDE is made entirely without reference to the histogram.

By adjusting the proportion of the sample plotted, the bandwidth of the KDE, and the sample size, you may be able to improve on my plot below. But the agreement of the density function and the KDE in the the plot below is roughly typical for samples of size 500 to 1000.

enter image description here

R code for figure:

k = diff(pt(c(-4,4),1))
set.seed(2022)
w = rt(1000, 1)        # whole sample, size 1000
y = w[x > -4 & w < 4]  # 861 plotted points
length(y)
[1] 861

hist(y, prob=T, ylim=c(0,.4), br=50, col="skyblue2")
 curve(dt(x,1)/k, add=T, lwd=2)
 lines(density(y), col="red", lwd=2, lty="dotted")

Note: For moderately large samples, the empirical CDF (ECDF) of the sample matches the density function of the population very well. The test statistic $D$ of the Kolmogorov-Smirnov goodness-of-fit test, is the maximum vertical distance between the two.

enter image description here

plot(ecdf(w), lwd=3, col="red", lty="dotted")
 curve(pt(x,1), add=T)
  abline(v=0, col="green2")
  abline(h=0:1, col="green2")

ks.test(w, pt, 1)

        One-sample Kolmogorov-Smirnov test

data:  w
D = 0.018508, p-value = 0.8832
alternative hypothesis: two-sided
$\endgroup$
6
  • 4
    $\begingroup$ Strictly speaking, there are no "outliers" in a Cauchy sample, all points within the sample are from the same Cauchy distribution. $\endgroup$
    – Xi'an
    Feb 2, 2022 at 7:07
  • 1
    $\begingroup$ Yes, I changed the wording a bit before I saw your comment. $\endgroup$
    – BruceET
    Feb 2, 2022 at 7:12
  • 1
    $\begingroup$ Thanks @BruceET. Certainly, I can understand how the code may be modified to make visually acceptable KDEs. I am more curious as to where the theory of KDEs brings down in the practical implementation for Cauchy. That is, why is it "usually necessary to disregard more than a few values in the tails of the sample" $\endgroup$ Feb 2, 2022 at 7:58
  • 1
    $\begingroup$ The optimal choice of bandwidth depends on the standard deviation. I may be crawling out on a limb here, but I think other parts of the theory of kernel density estimation may also depend on the existence of the population standard deviation. For now, I'll leave that for users who know more about KDE than I to clarify. $\endgroup$
    – BruceET
    Feb 2, 2022 at 8:18
  • 2
    $\begingroup$ @Greenparker: The only characterisation of the target density $f$ that matters in the mean integrated square error (MISE) and the choice of an optimal bandwidth is$$\int (f")^2(x)\,\text dx$$which is well-defined for the Cauchy density. $\endgroup$
    – Xi'an
    Feb 2, 2022 at 16:11
2
$\begingroup$

By default, the density function in R computes the density as a set of 512 points, rather than actually computing the density function. If you would like to obtain the actual density function you can use the KDE function in the utilities package. This function computes all the probability functions for the KDE and gives these probability functions as elements of the output object. It also allows you to easily load the probability functions to the global environment (or another environment).

#Load the library
library(utilities)

#Create data from the Cauchy distribution
set.seed(1)
DATA <- rt(1e3, df = 1)

#Generate the KDE and load probability functions to global environment
MY_KDE <- KDE(DATA, df = 1, to.environment = TRUE)
MY_KDE

  Kernel Density Estimator (KDE) 
 
  Computed from 1000 data points in the input 'DATA'
  Estimated bandwidth = 0.351773  
  Input degrees-of-freedom = 1.000000  
 
  Probability functions for the KDE are the following: 
 
      Density function:                   dkde * 
      Distribution function:              pkde * 
      Quantile function:                  qkde * 
      Random generation function:         rkde * 
 
  * This function is presently loaded in the global environment

Once you create the object and set to.environment = TRUE the four probability functions for the KDE are loaded to the global environment and can be called just like the probability functions for standard distributional families. We can now use the density function dkde to compare the Cauchy distribution to the KDE generated from data from the Cauchy distribution. As you can see in the plot, in the present case, there is a bit less data in the middle of the distribution than would be expected, but otherwise the fit is quite close.

#Plot the KDE against the Cauchy distribution
xx <- seq(-50, 50, length = 1e3)
y1 <- dt(xx, df = 1)
y2 <- dkde(xx)
plot(xx, y1, type = 'l', col = 'red',
     main = 'Cauchy Distribution vs KDE from Cauchy Data',
     xlab = 'Value', ylab = 'Density')
lines(xx, y2, type = 'l', col = 'blue')

enter image description here

$\endgroup$

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.