4
$\begingroup$

When plotting two kernel desitiy estimates into the same frame

plot(density(all4[,3],kernel="epanechnikov"),col='red',ylim=c(0,0.4),xlim=c(-50,50))
lines(density(all2[,3],kernel="epanechnikov"),col="blue",ylim=c(0,0.4),xlim=c(-50,50))

the result looks like this:

enter image description here

which is quite strange since the area under both plots should be unity. The data related to the red curve has an outlier (its maximum is 702.6).

> summary(all4[,3])
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-22.3100  -4.0940  -0.4151  -0.1533   3.4880 702.6000 
> summary(all2[,3])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-23.520  -4.041  -0.367  -0.109   3.533  67.780 

After removing the outlier the curves match pretty well

a4 = sort(all4[,3])
a2 = sort(all2[,3])
a4 = a4[6:999995]
a2 = a2[6:999995]
plot(density(a4,kernel="epanechnikov"),col='red',ylim=c(0,0.2),xlim=c(-50,50),lwd=3)
lines(density(a2,kernel="epanechnikov"),col="blue",ylim=c(0,0.2),xlim=c(-50,50))

> summary(a4)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-21.0200  -4.0940  -0.4151  -0.1541   3.4880  39.6700 
> summary(a2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-21.9200  -4.0410  -0.3670  -0.1091   3.5330  32.2900 

enter image description here

I am confused since the outlier should STEAL probability mass from the body of the distribution. But the red curve in the first picture seems to be inflated. Moreover, since both density curves are on the same scales, the area under both curves should be the same, also in the first plot. Where is my mistake?

EDIT: Following Momos recommendations I integrated over the densities. Indeed, the result for the weird one is far beyond unity!

library(sfsmisc)
dd1 = density(all4[,3],kernel="epanechnikov")
dd2 = density(all2[,3],kernel="epanechnikov")
integrate.xy(dd1$x,dd1$y)
[1] 1.490016
integrate.xy(dd2$x,dd2$y)
[1] 0.9977962
$\endgroup$
  • 1
    $\begingroup$ You could start with numerically integrating the two densities and see whether your assertion that they do not integrate to 1 holds. Its easily done with integrate. $\endgroup$ – Momo Jun 20 '14 at 18:56
  • 1
    $\begingroup$ Good idea, thanks! However integrate seems not to work since density(some data) is not a function. I used integrate.xy instead and there seem to be problems indeed! $\endgroup$ – Joz Jun 20 '14 at 19:17
4
$\begingroup$

The problem is discretization error.

density divides the range of data by default into n=512 intervals. The outlier causes the intervals to be at least an order of magnitude greater than otherwise. With very large amounts of data (it seems to take somewhere more than $10^5$ for this problem to become really apparent with the default value of n and this order-of-magnitude outlier), the cumulative discretization error takes over.

The general lesson--which is by no means confined to just density estimation nor to R--is that when routine tools are pushed towards their limits, we need to understand how they do their calculations and we should remember to pay attention to all the parameters that control their accuracy.


I haven't reverse-engineered the source code to pin this down definitively, but I will supply some empirical evidence using random data sets.

trial <- function(n, mean=0, sd=5, outlier=700, kernel="epanechnikov", width=512) {
  set.seed(17) # Makes everything reproducible
  x <- c(outlier, rnorm(n, mean, sd))
  #
  # Plot the densities for visual comparison.
  #
  plot(dd1 <- density(x, kernel=kernel, n=width), col="Red", xlim=c(-sd,sd)*5)
  lines(dd2 <- density(x[-1], kernel=kernel, n=width), col="Blue", xlim=c(-sd,sd)*5)
  #
  # Print crude integral estimates (trapezoidal rule) to memorialize the results.
  # 
  print (c(outlier=sum(diff(dd1$x)*(dd1$y[-1] + dd1$y[-length(dd1$y)])/2),
         non.outlier=sum(diff(dd2$x)*(dd2$y[-1] + dd2$y[-length(dd2$y)])/2)))
}
trial(10^6)                       # Looks like the error in the question
trial(10^6, kernel="rectangular") # Has similar results--problem is not the kernel
trial(10^6, width=2^13)           # Both integrals are essentially 1 now
trial(10^3)                       # Agreement is ok
trial(10^3, width=2^6)            # The problem is reproduced
trial(10^3, outlier=10^4)         # The problem is enormous due to the wide range
$\endgroup$
3
$\begingroup$

The only think I can think of is that the bandwidth selected with the outliers left in is very small and your "red" density is just a bunch of peaks at the data points. However, since it is only evaluated at n=512 points and the plot function just joins the function values it appears smoothed and with higher values everywhere.

Try this to see the effect

a <- rnorm(5000000)

plot(density(a))

lines(density(a, bw=0.001, kernel="epanechnikov"), col="red")

lines(density(a, bw=0.001, kernel="epanechnikov", n=500000), col="blue")
$\endgroup$
1
$\begingroup$

Probably the red line has a thinner tail (which extends infinitely in either direction as does the blue) and the "missing" area is captured there, which is hard to see in your graph, which focuses on the median/mode and is of too gross a scale to capture the tail difference.

$\endgroup$
  • $\begingroup$ Thanks for your suggestion, Avraham! I just edited the density-function to use an epanechnikov kernel, which has finite support. The problem remains. And more importantly, the deleting of the more extreme values in the red sample compared to the blue sample should make the tails thinner, thereby adding probability mass to the body of the distribution. But as can be seen from the plots there is a loss of probability mass in the body after removing the extreme values in the tails. I am still confused. $\endgroup$ – Joz Jun 20 '14 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.