机器算法验证 - R中非负变量密度图的好方法？ - 吾爱随笔录

R中非负变量密度图的好方法？

机器算法验证 r 密度函数伽马分布内核平滑

2022-01-25 14:46:09

plot(density(rexp(100))

显然，零左侧的所有密度都表示偏差。

我希望为非统计学家总结一些数据，并且我想避免关于为什么非负数据的密度在零左侧的问题。这些图用于随机化检查；我想按治疗组和对照组显示变量的分布。分布通常是指数级的。由于各种原因，直方图很棘手。

一个快速的谷歌搜索给了我统计学家关于非负内核的工作，例如： this。

但是有没有在 R 中实现过呢？在已实施的方法中，它们中的任何一个在描述性统计方面是否以某种方式“最好”？

编辑：即使该from命令可以解决我当前的问题，也很高兴知道是否有人根据非负密度估计的文献实现了内核

4个回答

另一种方法是 Kooperberg 及其同事的方法，基于使用样条估计密度来近似数据的对数密度。我将展示一个使用来自@whuber 答案的数据的示例，这将允许比较方法。

set.seed(17)
x <- rexp(1000)

您需要为此安装logspline包；如果不是，请安装它：

install.packages("logspline")

logspline()加载包并使用函数估计密度：

require("logspline")
m <- logspline(x)

在下文中，我假设d来自@whuber 答案的对象存在于工作区中。

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

结果图如下所示，对数样条密度由红线表示

默认、截断和对数样条密度

此外，可以通过参数lbound和指定对密度的支持ubound。如果我们希望假设 0 左侧的密度为 0 并且在 0 处存在不连续性，我们可以lbound = 0在调用中使用logspline()，例如

m2 <- logspline(x, lbound = 0)

产生以下密度估计（此处显示原始m对数样条拟合，因为上一个图已经很忙了）。

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

结果图如下所示

比较有和没有支持下限的对数样条密度估计

在这种情况下，利用知识x处不趋于0 ，但类似于其他地方的标准对数样条拟合 $x = 0$ x

从空间统计的边缘加权方法中借用的一种解决方案是将左侧的密度截断为零，但对最接近于零的数据进行加权。这个想法是每个值为中心的单位总面积的内核中；内核的任何会溢出到负区域的部分都被删除，内核被重新规范化为单位面积。 $x$ $x$

例如，对于高斯核，重整化权重为 $K_h(y,x) = \exp(-\frac{1}{2}((y-x)/h)^2) / \sqrt{2\pi}$

w (x) = 1 / \int_{0}^{\infty} K (y, x) d y = \frac{1}{1 - Φ_{x, h} (0)}

$w(x) = 1 / \int_0^\infty K(y,x) dy = \frac{1}{1 - \Phi_{x, h}(0)}$

其中是均值和标准差的正态变量的累积分布函数。可比较的公式可用于其他内核。 $\Phi$ $x$ $h$

的带宽更简单——计算速度也快得多。无论如何，很难准确地规定应该如何在附近更改带宽。尽管如此，这种方法也是特设附近仍然会有一些偏差。它似乎比默认的密度估计更好。这是使用较大数据集的比较： $0$ $0$ $0$

蓝色显示默认密度，而红色显示在处为边缘调整的密度。真实的基础分布以虚线形式跟踪以供参考。 $0$

R代码

中的density函数R会抱怨权重之和不是统一的，因为它希望所有实数的积分是统一的，而这种方法使正数的积分等于统一。作为检查，后一个积分被估计为黎曼和。

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)

要按组比较分布（您在其中一个评论中说这是目标），为什么不做一些更简单的事情呢？如果 N 很大，平行箱线图效果很好；如果 N 很小（并且两者都很好地显示异常值，您说这是您的数据中的一个问题），则平行条形图起作用。

正如 Stéphane 评论的那样，您可以使用from = 0，此外，您可以在密度曲线下表示您的值rug (x)

其它你可能感兴趣的问题

上一篇如何使用传统编程语言从具有已知均值和方差的正态分布中采样？下一篇在 R 中查找四分位数