在这个线程的其他地方,我提出了一个简单但有点特别的对点进行二次采样的解决方案。它速度很快,但需要一些实验才能产生出色的情节。即将描述的解决方案要慢一个数量级(120 万个点最多需要 10 秒),但具有自适应性和自动性。对于大型数据集,它应该在第一时间就给出好的结果,而且做得相当快。
这个想法是Douglas-Peucker折线简化算法的想法,适用于 QQ 图的特征。这种图的相关统计量是Kolmogorov-Smirnov 统计量 Dn,与拟合线的最大垂直偏差。因此,算法是这样的:
求连接极值的直线之间的最大垂直偏差( x , y)对和他们的QQ图。如果这在可接受的分数范围内t全系列的y,用这条线替换情节。否则,将数据划分为最大垂直偏差点之前和之后的数据,并将算法递归地应用于这两个部分。
有一些细节需要注意,尤其是处理不同长度的数据集。我通过用对应于较长分位数的分位数替换较短的分位数来做到这一点:实际上,使用较短的 EDF 的分段线性近似而不是其实际数据值。(“更短”和“更长”可以通过设置颠倒use.shortest=TRUE
。)
这是一个R
实现。
qq <- function(x0, y0, t.y=0.0005, use.shortest=FALSE) {
qq.int <- function(x,y, i.min,i.max) {
# x, y are sorted and of equal length
n <-length(y)
if (n==1) return(c(x=x, y=y, i=i.max))
if (n==2) return(cbind(x=x, y=y, i=c(i.min,i.max)))
beta <- ifelse( x[1]==x[n], 0, (y[n] - y[1]) / (x[n] - x[1]))
alpha <- y[1] - beta*x[1]
fit <- alpha + x * beta
i <- median(c(2, n-1, which.max(abs(y-fit))))
if (abs(y[i]-fit[i]) > thresh) {
assemble(qq.int(x[1:i], y[1:i], i.min, i.min+i-1),
qq.int(x[i:n], y[i:n], i.min+i-1, i.max))
} else {
cbind(x=c(x[1],x[n]), y=c(y[1], y[n]), i=c(i.min, i.max))
}
}
assemble <- function(xy1, xy2) {
rbind(xy1, xy2[-1,])
}
#
# Pre-process the input so that sorting is done once
# and the most detail is extracted from the data.
#
is.reversed <- length(y0) < length(x0)
if (use.shortest) is.reversed <- !is.reversed
if (is.reversed) {
y <- sort(x0)
n <- length(y)
x <- quantile(y0, prob=(1:n-1)/(n-1))
} else {
y <- sort(y0)
n <- length(y)
x <- quantile(x0, prob=(1:n-1)/(n-1))
}
#
# Convert the relative threshold t.y into an absolute.
#
thresh <- t.y * diff(range(y))
#
# Recursively obtain points on the QQ plot.
#
xy <- qq.int(x, y, 1, n)
if (is.reversed) cbind(x=xy[,2], y=xy[,1], i=xy[,3]) else xy
}
作为一个例子,我使用了我之前的答案中模拟的数据(这次抛出了一个极高的异常值并且这次y
污染更多):x
set.seed(17)
n.x <- 1.21 * 10^6
n.y <- 1.20 * 10^6
k <- floor(0.01*n.x)
x <- c(rnorm(n.x-k), rnorm(k, mean=2, sd=2))
x <- x[x <= -3 | x >= -2.5]
y <- c(rbeta(n.y, 10,13), 1)
让我们绘制几个版本,使用越来越小的阈值。在 0.0005 的值并在 1000 像素高的监视器上显示时,我们将保证图上各处的误差不超过垂直像素的二分之一。这以灰色显示(只有 522 个点,由线段连接);粗略的近似值被绘制在它上面:首先是黑色,然后是红色(红色点将是黑色点的子集并重叠绘制它们),然后是蓝色(这也是子集和重叠绘制)。时间范围从 6.5(蓝色)到 10 秒(灰色)。鉴于它们可以很好地缩放,人们不妨使用大约二分之一像素作为阈值的通用默认值(例如,1/2000 用于 1000 像素高的显示器)并完成它。
qq.1 <- qq(x,y)
plot(qq.1, type="l", lwd=1, col="Gray",
xlab="x", ylab="y", main="Adaptive QQ Plot")
points(qq.1, pch=".", cex=6, col="Gray")
points(qq(x,y, .01), pch=23, col="Black")
points(qq(x,y, .03), pch=22, col="Red")
points(qq(x,y, .1), pch=19, col="Blue")
编辑
我已经修改了原始代码,qq
以便将第三列索引返回到原始两个数组中最长(或最短,如指定)x
和y
,对应于所选点。这些索引指向数据的“有趣”值,因此可用于进一步分析。
我还删除了一个重复值的错误x
(导致beta
未定义)。