How to avoid smoothing to 0 at edges of R density plot

Question

When using the density function in R, it includes smooth transitions down to 0 at both ends of the data. Is there a way to prevent this? As a trivial example, suppose I am calculating the density function for 1000 uniformly spaced points between 0 and 1. What I'd like is a straight line from 0 to 1. Instead, I get a straight line from about 0.2 to 0.8 with smooth transitions down to 0. (Numbers are a little different for non-default kernels but the same general idea.) Of course, I know the density in this example, but not in real examples. I can shrink these end caps using bw or adjust, but that causes an undesirable reduction in smoothing within the interval.

Is there an option to truncate the kernel at the boundaries of the region, so that it estimates the density using only points inside the region? I have worked around this by mirroring all of my points about the two ends and then keeping only the part in the middle, but that seems like a crazy hack for something that should be simple. Is there a simpler way to do this?

Answer 1

In general, you need a routine which is told of any hard minimum and/or maximum for a variable and does somehow the right thing near those boundaries. The right thing could be mapping to a different space; estimating the density on a transformed scale, and then back-transforming; or reflecting probability mass backwards at a boundary.

Such routines are not common and the problem isn't even widely mentioned, so far as I know.

Otherwise put, no routine knows about a boundary unless it is told about it. Or, the software user must take extra care to do an equivalent calculation. Not the answer you seek, but for such data I think you are better off plotting the cumulative distribution function, or equivalently the quantile function. In this example, either would be a straight line for a uniform distribution, naturally. It's a broader issue, but my own prejudice is that density estimation is somewhat oversold, whereas quantile plots are still undersold. Density estimation smooths away trivial noise, to be sure, but for many very common situations, particularly highly skewed and/or bounded distributions, density estimation defaults often work poorly.

I am not well informed about density estimation code in R and don't know what adjust does. In any case this has been migrated to CV as essentially a statistical question. It seems clear that changing the kernel type or width is no kind of solution here, unless the kernel adapts at the boundaries, which is perhaps the nub of the question.

Answer 2

If it is ok for you to use a different smoother, you could try P-splines in these cases. The method I am referring to is presented in Eilers and Marx (1991).

Edit - how does this help?

Quoting Eilers and Marx,

the P-spline density smoother is not troubled by boundary effects, as for instance kernel smoothers are.

In general, P-splines combine B-splines and finite difference penalties. The density smoothing problem is a special case of GLM. So we just need to parameterize our smoothing problem accordingly.

The R-code below reproduces an example similar to the one proposed in the original question. The same code and a short explanation of the method can be found here: Kernel density estimation and boundary bias (Edit: N=60 in analogy with the other example)

# Simulate data
set.seed(1)
N = 60
x = runif(N)

# Construct histograms
his = hist(x, breaks = 60, plot = F, prob = T)
X = his$counts
u = his$mids

# Prepare basis (I-mat) and penalty (1st difference)
B = diag(length(X))
D1 = diff(B, diff = 1)
lambda = 1e6 # fixed but can be selected (e.g. AIC)
P = lambda * t(D1) %*% D1

# Smooth
tol = 1e-8
eta = log(X + 1)
for (it in 1:20) 
{
mu = exp(eta)
z = X - mu + mu * eta
a = solve(t(B) %*% (c(mu) * B) + P, t(B) %*% z)
etnew = B %*% a
de = max(abs(etnew - eta))
cat('Crit', it, de, '\n')
if(de < tol) break
eta = etnew
}

# Plot
plot(u, exp(eta) / diff(u)[1] / sum(exp(eta)), ylim = c(0, max(his$density)), type = 'l', col = 2)
lines(u, his$density, type = 'h')

The results look like this (you will notice there is no boundary bias):

Answer 3

This problem is discussed in Venables & Ripley MASS (the book) writes:

Most density estimators will not work well when the density is non-zero at an
end of its support, such as the exponential and half-normal densities. (They are
designed for continuous densities and this is discontinuity.) One trick is to reflect
the density and sample about the endpoint, say, a. Thus we compute the density
for the sample c(x, 2a-x) , and take double its density on [a, ∞) (or (−∞, a]
for an upper endpoint). This will impose a zero derivative on the estimated density at a, but the end effect will be much less severe. For details and further tricks
see Silverman (1986, §3.10). 

The cited book is here. Then they go on to mention boundary kernels, which did not exist in R/S-plus at that time. But now they do, see the following simple example:

library(bde)
set.seed(7*11*13) # My public seed
testdata <- runif(60)
bde.estimate <- bde::bde(testdata, estimator="boundarykernel")