density {stats}  R Documentation 
Kernel Density Estimation
Description
The (S3) generic function density
computes kernel density
estimates. Its default method does so with the given kernel and
bandwidth for univariate observations.
Usage
density(x, ...)
## Default S3 method:
density(x, bw = "nrd0", adjust = 1,
kernel = c("gaussian", "epanechnikov", "rectangular",
"triangular", "biweight",
"cosine", "optcosine"),
weights = NULL, window = kernel, width,
give.Rkern = FALSE, subdensity = FALSE,
warnWbw = var(weights) > 0,
n = 512, from, to, cut = 3, ext = 4,
old.coords = FALSE,
na.rm = FALSE, ...)
Arguments
x 
the data from which the estimate is to be computed. For the default method a numeric vector: long vectors are not supported. 
bw 
the smoothing bandwidth to be used. The kernels are scaled such that this is the standard deviation of the smoothing kernel. (Note this differs from the reference books cited below.)
The specified (or computed) value of 
adjust 
the bandwidth used is actually 
kernel , window 
a character string giving the smoothing kernel
to be used. This must partially match one of

weights 
numeric vector of nonnegative observation weights,
hence of same length as Note that weights are not taken into account for automatic
bandwidth rules, i.e., when 
width 
this exists for compatibility with S; if given, and

give.Rkern 
logical; if true, no density is estimated, and
the ‘canonical bandwidth’ of the chosen 
subdensity 
used only when 
warnWbw 

n 
the number of equally spaced points at which the density is
to be estimated. When 
from , to 
the left and rightmost points of the grid at which the
density is to be estimated; the defaults are 
cut 
by default, the values of 
ext 
a positive extension factor, 
old.coords 

na.rm 
logical; if 
... 
further arguments for (nondefault) methods. 
Details
The algorithm used in density.default
disperses the mass of the
empirical distribution function over a regular grid of at least 512
points and then uses the fast Fourier transform to convolve this
approximation with a discretized version of the kernel and then uses
linear approximation to evaluate the density at the specified points.
The statistical properties of a kernel are determined by
\sigma^2_K = \int t^2 K(t) dt
which is always = 1
for our kernels (and hence the bandwidth
bw
is the standard deviation of the kernel) and
R(K) = \int K^2(t) dt
.
MSEequivalent bandwidths (for different kernels) are proportional to
\sigma_K R(K)
which is scale invariant and for our
kernels equal to R(K)
. This value is returned when
give.Rkern = TRUE
. See the examples for using exact equivalent
bandwidths.
Infinite values in x
are assumed to correspond to a point mass at
+/Inf
and the density estimate is of the subdensity on
(Inf, +Inf)
.
Value
If give.Rkern
is true, the number R(K)
, otherwise
an object with class "density"
whose
underlying structure is a list containing the following components.
x 
the 
y 
the estimated density values. These will be nonnegative, but can be zero. 
bw 
the bandwidth used. 
n 
the sample size after elimination of missing values. 
call 
the call which produced the result. 
data.name 
the deparsed name of the 
has.na 
logical, for compatibility (always 
The print
method reports summary
values on the
x
and y
components.
References
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988). The New S Language. Wadsworth & Brooks/Cole (for S version).
Scott, D. W. (1992). Multivariate Density Estimation. Theory, Practice and Visualization. New York: Wiley.
Sheather, S. J. and Jones, M. C. (1991). A reliable databased bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society Series B, 53, 683–690. doi:10.1111/j.25176161.1991.tb01857.x.
Silverman, B. W. (1986). Density Estimation. London: Chapman and Hall.
Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S. New York: Springer.
See Also
bw.nrd
,
plot.density
, hist
;
fft
and convolve
for the computational short
cut used.
Examples
require(graphics)
plot(density(c(20, rep(0,98), 20)), xlim = c(4, 4)) # IQR = 0
# The Old Faithful geyser data
d < density(faithful$eruptions, bw = "sj")
d
plot(d)
plot(d, type = "n")
polygon(d, col = "wheat")
## Missing values:
x < xx < faithful$eruptions
x[i.out < sample(length(x), 10)] < NA
doR < density(x, bw = 0.15, na.rm = TRUE)
lines(doR, col = "blue")
points(xx[i.out], rep(0.01, 10))
## Weighted observations:
fe < sort(faithful$eruptions) # has quite a few nonunique values
## use 'counts / n' as weights:
dw < density(unique(fe), weights = table(fe)/length(fe), bw = d$bw)
utils::str(dw) ## smaller n: only 126, but identical estimate:
stopifnot(all.equal(d[1:3], dw[1:3]))
## simulation from a density() fit:
# a kernel density fit is an equallyweighted mixture.
fit < density(xx)
N < 1e6
x.new < rnorm(N, sample(xx, size = N, replace = TRUE), fit$bw)
plot(fit)
lines(density(x.new), col = "blue")
## The available kernels:
(kernels < eval(formals(density.default)$kernel))
## show the kernels in the R parametrization
plot (density(0, bw = 1), xlab = "",
main = "R's density() kernels with bw = 1")
for(i in 2:length(kernels))
lines(density(0, bw = 1, kernel = kernels[i]), col = i)
legend(1.5,.4, legend = kernels, col = seq(kernels),
lty = 1, cex = .8, y.intersp = 1)
## show the kernels in the S parametrization
plot(density(0, from = 1.2, to = 1.2, width = 2, kernel = "gaussian"),
type = "l", ylim = c(0, 1), xlab = "",
main = "R's density() kernels with width = 1")
for(i in 2:length(kernels))
lines(density(0, width = 2, kernel = kernels[i]), col = i)
legend(0.6, 1.0, legend = kernels, col = seq(kernels), lty = 1)
## Semiadvanced theoretic from here on 
## Explore the old.coords TRUE > FALSE change:
set.seed(7); x < runif(2^12) # N = 4096
den < density(x) # > grid of n = 512 points
den0 < density(x, old.coords = TRUE)
summary(den0$y / den$y) # 1.001 ... 1.011
summary( den0$y / den$y  1) # ~= 1/(2n2)
summary(1/ (den0$y / den$y  1))# ~= 2n2 = 1022
corr0 < 1  1/(2*5122) # 1  1/(2n2)
all.equal(den$y, den0$y * corr0)# ~ 0.0001
plot(den$x, (den0$y  den$y)/den$y, type='o', cex=1/4)
title("relative error of density(runif(2^12), old.coords=TRUE)")
abline(h = 1/1022, v = range(x), lty=2); axis(2, at=1/1022, "1/(2n2)", las=1)
## The R[K] for our kernels:
(RKs < cbind(sapply(kernels,
function(k) density(kernel = k, give.Rkern = TRUE))))
100*round(RKs["epanechnikov",]/RKs, 4) ## Efficiencies
bw < bw.SJ(precip) ## sensible automatic choice
plot(density(precip, bw = bw),
main = "same sd bandwidths, 7 different kernels")
for(i in 2:length(kernels))
lines(density(precip, bw = bw, kernel = kernels[i]), col = i)
## Bandwidth Adjustment for "Exactly Equivalent Kernels"
h.f < sapply(kernels, function(k)density(kernel = k, give.Rkern = TRUE))
(h.f < (h.f["gaussian"] / h.f)^ .2)
## > 1, 1.01, .995, 1.007,... close to 1 => adjustment barely visible..
plot(density(precip, bw = bw),
main = "equivalent bandwidths, 7 different kernels")
for(i in 2:length(kernels))
lines(density(precip, bw = bw, adjust = h.f[i], kernel = kernels[i]),
col = i)
legend(55, 0.035, legend = kernels, col = seq(kernels), lty = 1)