R中的高斯核密度估计

Question

我无法理解如何在 R 中实现以下数据集的高斯核密度估计。如果你能帮助我理解如何做的机制，我将不胜感激。我目前正在尝试获得下图底部钟形曲线的公式。如您所见，每个数据点都有一个钟形曲线。（注意图片不代表我使用的数据。）

这是我的数据：

x<-c(4.09, 4.46, 4.61, 4.30, 4.03, 5.22, 4.21, 4.07, 4.02, 4.58, 4.66, 4.05, 4.23, 5.51, 4.03, 4.72, 4.47, 4.50, 5.80, 4.30, 4.09, 4.78, 4.18, 4.45, 4.40, 5.60, 4.37, 4.42, 4.88, 4.20, 4.45, 4.10, 4.43, 4.58, 4.40, 4.38) (x 有 36 个元素)

这是核密度估计器：

（如果你看不到图片，它来自这个页面http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/tutorials/xlghtmlnode33.html）

其中 K(u)=

是高斯核函数，h=.1516是Scott选择的带宽

所以，插入我们得到 f hat (x) = 1/(36*.1516) (1/sqrt(2pi))[e^(-1/2 ((4.09-x)/.1516 )^2 + e^(-1/2 ((4.46-x)/.1516)^2 + ... + e^(-1/2 ((4.38-x)/.1516)^2]

好的。所以我们有 x 的函数。但是我们如何得到上图中每条钟形曲线的方程呢？例如，如果我们将 4.09 代入 f hat (x)，我们得到一个数字，而不是 curve/function/distribution。有人可以帮助我理解找到钟形 curve/kernel 密度估计方程的过程吗？

Answer 1

根据您的 x 值和 h 值

，此函数将 return 您的 fhat 函数

get_fhat <- function(x, h) {
  Vectorize(function(z) 1/length(x)/h*sum(dnorm((x-z)/h)))  
}

这个函数return是一个我们可以用来获取值的函数。我们 Vectorize 它因此我们可以一次将多个值传递给函数。

我们可以获得单个值或用

绘制它

fhat <- get_fhat(x, .1516)
fhat(4.09)
# [1] 0.9121099
curve(fhat, from=min(x), to=max(x))

Answer 2

图表

## Given data
x  <- c(4.09, 4.46, 4.61, 4.30, 4.03, 5.22, 4.21, 4.07, 4.02, 4.58, 4.66, 4.05, 
        4.23, 5.51, 4.03, 4.72, 4.47, 4.50, 5.80, 4.30, 4.09, 4.78, 4.18, 4.45, 
        4.40, 5.60, 4.37, 4.42, 4.88, 4.20, 4.45, 4.10, 4.43, 4.58, 4.40, 4.38)
h  <- 0.1516 

# GaussianKernel
GK <- function(u) {(1/sqrt(2*pi))*exp(-(u^2)/2)} # or dnorm(u)

这个函数给出了类似的图。

DensityGraph <- function(x, h){
  n    <- length(x)
  xi   <- seq(min(x) - sd(x), max(x) + sd(x), length.out = 512)
  # fhat without sum since we are interest in the bell shaped curves
  fhat <- sapply(x, function(y){(1/(n*h))*GK((xi - y)/h)})
  # histogram of x
  hist (x, freq = FALSE, nclass = 15, main = "Kernel density with histogram",
        xlab = paste("N = ", n, "   ", "Bandwidth = ", h))
  # add fhat with sum
  lines(xi, rowSums(fhat), lwd = 2)
  # add the bell shaped curves
  apply(fhat, 2, function(j) lines(xi, j, col = 4))
  # show data points
  rug  (x, lwd = 2, col = 2)
}


DensityGraph(x = x, h = 0.05)

蓝色钟形曲线代表 x 的每个数据点

DensityGraph(x = x, h = 0.1516)

与R中内置的密度函数比较

lines(density(x = x, bw = 0.1516), col = 3, lwd = 2)

每个x的fhat

此函数给出给定特定 x 的 fhat 的值

fhat <- function(x, h, specific_x){
  n    <- length(x)
  xi   <- seq(min(x) - sd(x), max(x) + sd(x), length.out = 512)
  f    <- rowSums(sapply(x, function(y){(1/(n*h))*GK((xi - y)/h)}))
  kde  <- data.frame(xi, fhat = f)
  indx <- which.min(abs(xi - specific_x))
  fx   <- kde[indx, "fhat"]
  list(fx = fx, kde = kde)
}

KernelDensity <- fhat(x = x, h = 0.1516, specific_x = 4.09)
KernelDensity$fx
# [1] 0.9114677
plot(KernelDensity$kde, type  = "l", lwd = 2, xlab = "")
title(xlab = paste("N = ", n, "    Bandwidth = ", h))
rug(x, lwd = 2, col = 2)

比较内置密度函数

lines(density(x, bw = 0.1516), col = 5)

R中的高斯核密度估计

Gaussian kernel density estimation in R

r

kernel-density

图表

每个x的fhat