如何找到经验累积密度函数 (ECDF) 的分位数
How to find quantiles of an empirical cumulative density function (ECDF)
我正在使用 ecdf() 函数从一些随机样本中计算经验累积密度函数 (ECDF)：
set.seed(0)
X = rnorm(100)
P = ecdf(X)
现在 P 给出了 ECDF，我们可以绘制它：
plot(P)
abline(h = 0.6, lty = 3)
我的问题是：如何找到样本值x，使得P(x) = 0.6，即ECDF的0.6分位数，或者ECDF 与 h = 0.6?
交点的 x 坐标
下面我就不使用ecdf()了，经验累积密度函数（ECDF）很容易自己求得
首先，我们对样本X进行升序排序：
X <- sort(X)
这些样本的 ECDF 取函数值：
e_cdf <- 1:length(X) / length(X)
然后我们可以通过以下方式绘制 ECDF：
plot(X, e_cdf, type = "s")
abline(h = 0.6, lty = 3)
现在，我们正在寻找 X 的第一个值，这样 P(X) >= 0.6。这只是：
X[which(e_cdf >= 0.6)[1]]
# [1] 0.2290196
由于我们的数据是从标准正态分布中抽样的，所以理论分位数是
qnorm(0.6)
# [1] 0.2533471
所以我们的结果非常接近。
分机
因为CDF的反函数是分位数函数（比如pnorm()的反函数是qnorm()），所以可以猜到ECDF的反函数作为样本分位数，即 ecdf() 的逆是 quantile()。这不是真的！
ECDF 是阶梯/阶跃函数，它没有反函数。如果我们围绕 y = x 旋转 ECDF，得到的曲线不是数学函数。 所以样本分位数与ECDF无关.
对于n排序的样本，样本分位数函数实际上是(x, y)的线性插值函数，其中：
- x 值为
seq(0, 1, length = n)；
- y 值正在排序样本中。
我们可以通过以下方式定义我们自己版本的样本分位数函数：
my_quantile <- function(x, prob) {
if (is.unsorted(x)) x <- sort(x)
n <- length(x)
approx(seq(0, 1, length = n), x, prob)$y
}
我们来做个测试：
my_quantile(X, 0.6)
# [1] 0.2343171
quantile(X, prob = 0.6, names = FALSE)
# [1] 0.2343171
请注意，结果与我们从 X[which(e_cdf >= 0.6)[1]] 获得的结果不同。
正是出于这个原因，我拒绝在我的回答中使用 quantile()。
我正在使用 ecdf() 函数从一些随机样本中计算经验累积密度函数 (ECDF)：
set.seed(0)
X = rnorm(100)
P = ecdf(X)
现在 P 给出了 ECDF，我们可以绘制它：
plot(P)
abline(h = 0.6, lty = 3)
我的问题是：如何找到样本值x，使得P(x) = 0.6，即ECDF的0.6分位数，或者ECDF 与 h = 0.6?
下面我就不使用ecdf()了，经验累积密度函数（ECDF）很容易自己求得
首先，我们对样本X进行升序排序：
X <- sort(X)
这些样本的 ECDF 取函数值：
e_cdf <- 1:length(X) / length(X)
然后我们可以通过以下方式绘制 ECDF：
plot(X, e_cdf, type = "s")
abline(h = 0.6, lty = 3)
现在，我们正在寻找 X 的第一个值，这样 P(X) >= 0.6。这只是：
X[which(e_cdf >= 0.6)[1]]
# [1] 0.2290196
由于我们的数据是从标准正态分布中抽样的，所以理论分位数是
qnorm(0.6)
# [1] 0.2533471
所以我们的结果非常接近。
分机
因为CDF的反函数是分位数函数（比如pnorm()的反函数是qnorm()），所以可以猜到ECDF的反函数作为样本分位数，即 ecdf() 的逆是 quantile()。这不是真的！
ECDF 是阶梯/阶跃函数，它没有反函数。如果我们围绕 y = x 旋转 ECDF，得到的曲线不是数学函数。 所以样本分位数与ECDF无关.
对于n排序的样本，样本分位数函数实际上是(x, y)的线性插值函数，其中：
- x 值为
seq(0, 1, length = n)； - y 值正在排序样本中。
我们可以通过以下方式定义我们自己版本的样本分位数函数：
my_quantile <- function(x, prob) {
if (is.unsorted(x)) x <- sort(x)
n <- length(x)
approx(seq(0, 1, length = n), x, prob)$y
}
我们来做个测试：
my_quantile(X, 0.6)
# [1] 0.2343171
quantile(X, prob = 0.6, names = FALSE)
# [1] 0.2343171
请注意，结果与我们从 X[which(e_cdf >= 0.6)[1]] 获得的结果不同。
正是出于这个原因，我拒绝在我的回答中使用 quantile()。