r - 获得 ECDF 的导数

Question

是否可以区分 ECDF？以下面得到的为例。

set.seed(1)

a <- sort(rnorm(100))
b <- ecdf(a)

plot(b)

我想对 b 求导以获得它的概率密度函数 (PDF)。

Answer 1

n <- length(a)  ## `a` must be sorted in non-decreasing order already
plot(a, 1:n / n, type = "s")  ## "staircase" plot; not "line" plot

However I'm looking to find the derivative of b

在基于样本的统计中，估计密度(对于连续随机变量)不是通过微分从 ECDF 获得的，因为样本大小是有限的，并且 ECDF 不可微分。相反，我们直接估计密度。我想 plot(density(a)) 才是您真正想要的。

---

几天后..

警告:以下只是数值解，没有统计依据!

我把它作为一个练习来了解 R 包 scam for shape constrained additive models，Wood 教授的 mgcv 的子包早期博士生 Pya 博士。

逻辑是这样的:

• 使用 scam::scam，将单调递增的 P 样条拟合到 ECDF(您必须指定您想要的结数)； [请注意，单调性并不是唯一的理论约束。要求平滑的 ECDF 在其两个边缘上被“剪裁”:左边缘为 0，右边缘为 1。我目前正在使用 weights 来施加这样的约束，通过赋予非常大的权重在两个边缘]
• 使用 stats::splinefun，通过节点的单调插值样条和节点处的预测值重新参数化拟合样条；
• 返回插值样条函数，它也可以计算一阶、二阶和三阶导数。

为什么我希望它起作用:

随着样本量的增加，

• ECDF 收敛于 CDF；
• P 样条是一致的，因此平滑的 ECDF 对于 ECDF 将越来越无偏；
• 平滑 ECDF 的一阶导数对于 PDF 将越来越无偏。

谨慎使用:

• 您必须自己选择结数；
• 导数未归一化，因此曲线下面积为1；
• 结果可能相当不稳定，仅适用于大样本量。

函数参数:

• x:样本向量；
• n.knots:结数；
• n.cells:绘制导函数时的网格点数

您需要从 CRAN 安装 scam 包。

library(scam)

test <- function (x, n.knots, n.cells) {

  ## get ECDF
  n <- length(x)
  x <- sort(x)
  y <- 1:n / n
  dat <- data.frame(x = x, y = y)  ## make sure `scam` can find `x` and `y`

  ## fit a monotonically increasing P-spline for ECDF
  fit <- scam::scam(y ~ s(x, bs = "mpi", k = n.knots), data = dat,
                    weights = c(n, rep(1, n - 2), 10 * n))
  ## interior knots
  xk <- with(fit$smooth[[1]], knots[4:(length(knots) - 3)])
  ## spline values at interior knots
  yk <- predict(fit, newdata = data.frame(x = xk))
  ## reparametrization into a monotone interpolation spline
  f <- stats::splinefun(xk, yk, "hyman")

  par(mfrow = c(1, 2))

  plot(x, y, pch = 19, col = "gray")  ## ECDF
  lines(x, f(x), type = "l")          ## smoothed ECDF
  title(paste0("number of knots: ", n.knots,
               "\neffective degree of freedom: ", round(sum(fit$edf), 2)),
        cex.main = 0.8)

  xg <- seq(min(x), max(x), length = n.cells)
  plot(xg, f(xg, 1), type = "l")     ## density estimated by scam
  lines(stats::density(x), col = 2)  ## a proper density estimate by density

  ## return smooth ECDF function
  f
  }

---

## try large sample size
set.seed(1)
x <- rnorm(1000)
f <- test(x, n.knots = 20, n.cells = 100)

f 是 stats::splinefun 返回的函数(阅读 ?splinefun)。

一个天真的、类似的解决方案是在 ECDF 上做插值样条而不平滑。但这是一个非常糟糕的主意，因为我们没有一致性。

g <- splinefun(sort(x), 1:length(x) / length(x), method = "hyman")
curve(g(x, deriv = 1), from = -3, to = 3)

提醒:强烈建议使用 stats::density 进行直接密度估计。