机器算法验证 - 柯西的核密度估计 - 吾爱随笔录

柯西的核密度估计

机器算法验证 r 密度估计柯西分布

2022-03-04 04:32:42

据我了解，核密度估计不对底层密度的矩做任何假设，只需要平滑度。柯西密度函数非常平滑。尽管如此，当我尝试density()在 R 中使用 KDE 来从柯西分布中随机抽取时，我得到的答案非常不准确：

set.seed(1)
foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(density(rt(1e3, df = 1)), col = "red")

产生这个情节：

对不同的种子重复上述操作或增加样本量可能会给出进一步的不稳定估计。默认内核是高斯 in R。将内核更改为任何其他选项都不会改善输出。

问题：柯西对 KDE 违反了哪些假设？如果没有，那么为什么我们会看到 KDE 在这里失败得如此悲惨？

编辑： @cdalitz 已经确定问题出在 kde 评估密度的位置。默认值为3*bw*range(x)，对于 Cauchy 来说可能非常大。这意味着，默认情况下density会尝试512在 x 轴上稀疏分布的点处估计 KDE。

为了测试这一点，我更改了密度估计中的fromandto并查看如果我density使用两组评估点运行两次，那么密度匹配：

set.seed(1)
samp <- rt(1e4, df = 1)

bd <- 10
den1 <- density(samp, from=-bd, to=bd, n=512)
den2 <- density(samp, from =-2*bd, to = 2*bd,  n =512)

foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(den2, col = "blue", type = "b")
lines(den1, col = "red", type = "b")

这会产生以下估计：

这里的质量比以前好多了。但是，现在如果不是2*bd，我将其更改为50*bd，我得到即使在 0 左右的密度估计也有很大不同！

set.seed(1)
samp <- rt(1e4, df = 1)

bd <- 10
den1 <- density(samp, from=-bd, to=bd, n=512)
den2 <- density(samp, from =-50*bd, to = 50*bd,  n =512)

foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l', ylim = c(0,.7))
lines(den2, col = "blue", type = "b")
lines(den1, col = "red", type = "b")

稀疏点的密度评估如何改变周围的密度评估过程 $x = 0$ den1（为和选择的带宽相同den2）？任意点的 KD 估计 $x$ 是

\hat{f} (x) = \frac{1}{n h} \sum_{t = 1}^{n} K (\frac{x - x_{i}}{h}) .

$\hat{f}(x) = \dfrac{1}{nh} \sum_{t=1}^{n} K\left( \dfrac{x - x_i}{h}\right)\,.$

密度估计值不应在给定值下改变 $x = a_1$ 如果密度也在其他点进行评估。我在这里想念什么？

2个回答

问题在于 x 范围的采样density。默认情况下，density对数据范围内的 512 个值进行采样。由于 Cauchy 分布的尾部很重，因此您只能使用默认设置对密度函数进行非常稀疏的采样。

如果我们在 -20 和 +20 之间采样 512 个值，密度近似值是不错的，即使使用默认带宽选择规则：

set.seed(1)
foo <- seq(-50, 50, length = 1e3)
plot(foo, dt(foo, df = 1), type = 'l')
lines(density(rt(1e3, df = 1), from=-20, to=20, n=512), col = "red")

另请注意from，to和n仅影响评估核密度估计器的样本点，对带宽没有影响：

> set.seed(1)
> x <- rt(1e3, df = 1)
> density(x, from=-20, to=20, n=512)$bw
[1] 0.351773
> density(x)$bw
[1] 0.351773

在 R 中，该过程density给出了数据的核密度估计器 (KDE)。

为了获得非常重尾柯西分布的合适直方图，通常需要忽略样本尾部中的多个值。

然后为了显示实际的柯西密度曲线（下图黑色），需要通过除以直方图中绘制的样本比例来调整密度。

此外，要获得最佳 KDE（红色虚线），您可能需要使用参数 adjtp 调整默认带宽。（我使用默认值。）当然，如果柱状图较少，直方图会“更平滑”。KDE 完全没有参考直方图。

通过调整绘制的样本比例、KDE 的带宽和样本大小，您也许可以改进我下面的图。但是下图中的密度函数和 KDE 的一致性对于大小为 500 到 1000 的样本大致是典型的。

图的R代码：

k = diff(pt(c(-4,4),1))
set.seed(2022)
w = rt(1000, 1)        # whole sample, size 1000
y = w[x > -4 & w < 4]  # 861 plotted points
length(y)
[1] 861

hist(y, prob=T, ylim=c(0,.4), br=50, col="skyblue2")
 curve(dt(x,1)/k, add=T, lwd=2)
 lines(density(y), col="red", lwd=2, lty="dotted")

注意：对于中等大小的样本，样本的经验 CDF (ECDF) 与总体的密度函数非常匹配。检验统计量 $D$ Kolmogorov-Smirnov 拟合优度检验的最大垂直距离是两者之间的最大垂直距离。

plot(ecdf(w), lwd=3, col="red", lty="dotted")
 curve(pt(x,1), add=T)
  abline(v=0, col="green2")
  abline(h=0:1, col="green2")

ks.test(w, pt, 1)

        One-sample Kolmogorov-Smirnov test

data:  w
D = 0.018508, p-value = 0.8832
alternative hypothesis: two-sided

其它你可能感兴趣的问题

上一篇我怎么知道选择哪种参数估计方法？下一篇做什么磷( |Xn- X| ≥ε)P(|Xn−X|≥ϵ)直观地表示？