R中的高斯核密度估计

Posted 2023-02-14

【中文标题】R中的高斯核密度估计

【英文标题】：Gaussian kernel density estimation in R

【发布时间】：2021-01-21 22:08:36

【问题描述】：

我无法理解如何在 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)，我们会得到一个数字，而不是曲线/函数/分布。有人可以帮我理解找到钟形曲线/核密度估计方程的过程吗？

【问题讨论】：

如果您插入x 的一个值，您将得到一个响应。要获得曲线形状，请插入许多不同的 x 值。当你连接结果时，你会得到你的曲线。

谢谢。我假设你的意思是整体曲线？但是如何获得核密度估计器，例如第一个数据点。

单个点并没有真正的密度估计。你到底是什么意思？

底部的曲线只是以每个数据点为中心的高斯分布，sd 等于所选带宽。

是的，我指的是底部的曲线。例如，第一个点上面有一个高斯分布。

【参考方案1】：

这是一个函数，它将根据您的 x 值和 h 值返回您的 fhat 函数

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

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

我们可以得到一个单一的值或用它来绘制它

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

【讨论】：

如何确定间隔的宽度？跟带宽h有关系吗？

如果您愿意，您可以将末端延长到 fhat 低于某个阈值。但是，由于您是对正常密度求和，因此它们在技术上具有无限的支持。

【参考方案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)