R density 核密度估计

说明

(S3) 通用函数 density 计算核密度估计。它的默认方法使用给定的核和单变量观测的带宽来实现这一点。

用法

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, ...)

参数

x	用于计算估计值的数据。对于默认方法，数字向量：不支持长向量。
bw	要使用的平滑带宽。对内核进行缩放，使其成为平滑内核的标准差。 (请注意，这与下面引用的参考书以及S-PLUS不同。) bw也可以是给出选择带宽的规则的字符串。看bw.nrd. 默认情况下，"nrd0"，出于历史和兼容性原因仍然是默认值，而不是作为一般建议，例如，"SJ"更适合，另见 Venables 和 Ripley (2002)。 bw 的指定(或计算)值乘以 adjust 。
adjust	使用的带宽实际上是 adjust*bw 。这样可以轻松指定“默认带宽的一半”等值。
kernel , window	给出要使用的平滑内核的字符串。这必须部分匹配 "gaussian" 、 "rectangular" 、 "triangular" 、 "epanechnikov" 、 "biweight" 、 "cosine" 或 "optcosine" 之一，默认为 "gaussian" ，并且可以缩写为唯一的前缀(单个信)。 "cosine" 比 "optcosine" 更平滑，这是文献中常见的 ‘cosine’ 内核，几乎是MSE-efficient。然而"cosine"是S使用的版本。
weights	非负观察权重的数值向量，因此与 x 长度相同。默认 NULL 相当于 weights = rep(1/nx, nx) ，其中 nx 是 x[] (的有限条目)的长度。如果 na.rm = TRUE 和 x 中有 NA ，则在计算之前删除它们和相应的权重。在这种情况下，当原始权重总计为 1 时，它们会重新缩放以继续这样做。 请注意，自动带宽规则不考虑权重，即当 bw 是字符串时。当权重与真实计数 cn 成比例时，可以使用 density(x = rep(x, cn)) 代替 weights 。
width	这是为了与 S 兼容而存在的；如果给定，而 bw 没有，则将 bw 设置为 width(如果这是字符串)，或者设置为 width 的 kernel-dependent 倍数(如果这是数字)。
give.Rkern	逻辑性；如果为 true，则不估计密度，而是返回所选 kernel 的“规范带宽”。
subdensity	仅当指定weights且其总和不为一时才使用。如果为 true，则表示需要 “sub-density”，并且不应发出任何警告信号。默认情况下，如果为 false，则当权重之和不等于 1 时，会发出 warning 信号。
warnWbw	logical ，仅在指定 weights 且 bw 为 character 时使用，即选择自动带宽选择(默认情况下)。如果为 true(默认情况下)，则会发出 warning 信号，提醒用户自动带宽选择不会考虑权重，因此可能不是最佳选择。
n	要估计密度的等距点的数量。当 n > 512 时，在计算过程中会向上舍入到 2 的幂(因为使用了 fft )，并且最终结果由 approx 进行插值。因此，将 n 指定为 2 的幂几乎总是有意义的。
from , to	要估计密度的网格的最左边和最右边的点；默认值是 range(x) 之外的 cut * bw 。
cut	默认情况下，from 和to 的值是超出数据极值的cut 带宽。这使得估计密度在极端情况下下降到大约为零。
ext	正扩展因子，默认为4。值from和to在两侧进一步扩展为lo <- from - ext * bw和up <- to + ext * bw，然后用于构建用于FFT和插值的网格，请参见上面的n。除非您知道自己在做什么，否则不要改变！
old.coords	logical 需要 R 4.4.0 之前的行为，该行为会给出大约 (1 + 1/(2n-2)) 的过大值。
na.rm	逻辑性；如果 TRUE ，则从 x 中删除缺失值。如果 FALSE 任何缺失值都会导致错误。
...	(非默认)方法的进一步参数。

细节

density.default 中使用的算法将经验分布函数的质量分散在至少 512 个点的规则网格上，然后使用快速傅里叶变换将此近似值与核的离散版本进行卷积，然后使用线性近似来评估指定点处的密度。

内核的统计特性由下式决定：\sigma^2_K = \int t^2 K(t) dt这总是= 1对于我们的内核(因此带宽bw是核的标准差)和R(K) = \int K^2(t) dt.
MSE-equivalent 带宽(针对不同内核)与\sigma_K R(K)这是尺度不变的，对于我们的内核来说等于R(K)。当以下情况时返回该值give.Rkern = TRUE。请参阅使用精确等效带宽的示例。

假定 x 中的无限值对应于 +/-Inf 处的点质量，并且密度估计是 (-Inf, +Inf) 上的 sub-density 的。

值

如果 give.Rkern 为 true，则为数字 R(K) ，否则为类 "density" 的对象，其底层结构是包含以下组件的列表。

x	估计密度的点的 n 坐标。
y	估计的密度值。这些将为非负数，但可以为零。
bw	使用的带宽。
n	消除缺失值后的样本量。
call	产生结果的调用。
data.name	x 参数的解析名称。
has.na	逻辑，为了兼容性(始终 FALSE )。

print 方法报告 x 和 y 组件上的 summary 值。

例子

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 non-unique 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 equally-weighted 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)

##-------- Semi-advanced 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/(2n-2)
summary(1/ (den0$y / den$y - 1))# ~=    2n-2 = 1022
corr0 <- 1 - 1/(2*512-2) # 1 - 1/(2n-2)
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/(2n-2)", 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)

参考

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 data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society Series B, 53, 683-690. doi:10.1111/j.2517-6161.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.

也可以看看

bw.nrd、plot.density、hist； fft 和 convolve 用于使用的计算捷径。