我正在寻找一种计算“阈值”的方法，其中一个经验分布的分位数等于另一个经验分布的 1-分位数。我曾认为这类似于基于估计重叠内核密度问题的区域的先前问题，或者可以通过内核密度的交集来估计。我不再相信是这样的。
这是一个玩具示例，我可以通过猜测得到一个近似值：
set.seed(11)
d <- data.frame(f=c(rep("a", 50), rep("b", 60)), x=c(rnorm(50), runif(60, 0, 3)))
ddply(d, .(f), summarize,ql=quantile(x, .1741),qu=quantile(x, .8259))
f ql qu
1 a -0.9932201 0.6501115
2 b 0.6501501 2.2738458
这就是与内核密度相比阈值所在的位置：
我想对成对比较矩阵执行此操作，例如 10 种不同的经验分布。所以我肯定需要一个比猜测更好的方法。
因为您将针对对分布执行此操作，所以您需要一种相当有效的方法。
该问题要求求解（至少近似）形式为的方程，其中是逆经验 CDF。等效地，您可以解决其中是经验 CDF。最好使用不假设函数是可微分（甚至连续）的求根方法来完成，因为这些函数是不连续的：它们在数据值处跳转。
在R，uniroot将完成这项工作。尽管它假设函数是连续的（我相信它使用布伦特方法），R但经验 CDF 的实现使它们看起来足够连续。要使此方法起作用，您需要将根括在已知边界之间，但这很容易：它必须位于两个数据集的并集范围内。
代码非常简单：给定两个数据数组x和y，创建它们的经验 CDF 函数F.x和F.y，然后调用uniroot。这就是你所需要的。
overlap <- function(x, y) {
F.x <- ecdf(x); F.y <- ecdf(y)
z <- uniroot(function(z) F.x(z) + F.y(z) - 1, interval<-c(min(c(x,y)), max(c(x,y))))
return(list(Root=z, F.x=F.x, F.y=F.y))
}
它相当快：应用于大小从到对 10 个数据集，它总共在秒内找到了答案。
或者，请注意所需的点是两个分布的相等混合的中位数。 当两个数据集大小相同时，只需获取所有数据并集的中位数即可！您可以通过计算加权中位数将其推广到不同大小的数据集。此功能可通过分位数回归（在quantreg包中）获得，该回归可容纳权重：根据常数回归数据并与数据集的大小成反比加权。
overlap.rq <- function(x, y) {
library(quantreg)
fit <- rq(c(x,y) ~ 1, data=d,
weights=c(rep(1/length(x), length(x)), rep(1/length(y), length(y))))
return(coef(fit))
}
时序测试表明，这至少比寻根方法慢三倍，并且对于更大的数据集也不能很好地扩展：在前面的对数据集的测试中，它花费了秒，慢了十倍以上。主要优点是加权中位数的这种特殊实现将在答案看起来不唯一时发出警告，而布伦特的方法倾向于在可能答案的区间中间找到唯一答案。
作为演示，这里是两个经验 CDF 的图以及显示两个解决方案的垂直线（以及标记和水平的水平线）。在这种特殊情况下，这两种方法产生相同的答案，因此只出现一条垂直线。
#
# Generate some data.
#
set.seed(17)
x <- rnorm(32, 5, 2)
y <- rgamma(10, 2)
#
# Compute the solution two ways.
#
solution <- overlap(x, y)
solution.rq <- overlap.rq(x, y)
F.x <- solution$F.x; F.y <- solution$F.y; z <- solution$Root
alpha <- c(F.x(z$root), F.y(z$root))
#
# Plot the ECDFs and the results.
#
plot(interval, 0:1, type="n", xlab="z", ylab="Probability", main="CDFs")
curve(F.x(x), add=TRUE, lwd=2, col="Red")
curve(F.y(x), add=TRUE, lwd=2, col="Blue")
abline(v=z$root, lty=2)
abline(v=solution.rq, lty=2, col="Green")
abline(h=alpha, lty=3, col="Gray")
我想到了使用经验累积分布函数。答案近似于任何所需的有效数字程度。这是我想出的：
CDF.intersect<-function(a, b){
#a and b are vectors of the same metric, intent is to find cdf
if(median(a) < median(b)){
Fn1<-ecdf(a)
Fn2<-ecdf(b)
} else{
Fn1<-ecdf(b)
Fn2<-ecdf(a)
}
x<-seq(min(c(a,b)), max(c(a,b)), length.out=100000)
for (i in 1:100000){
y<-(1-Fn1(x[i]))-Fn2(x[i])
z<-x[i]
if (y<=0.00001) break
}
out<-data.frame("Threshold"=z, "Upper Quantile of Lower Distribution"= 1-Fn1(z),
"Lower Quantile of Upper Distribution" = Fn2(z))
return(out)
}