去除 QQ 图中心附近的无关点

机器算法验证 r 数据可视化 QQ图
2022-02-10 07:30:27

我正在尝试在 R 中绘制一个包含两个大约 120 万个点的数据集的 QQ 图(使用 qqplot,并将数据输入 ggplot2)。计算很简单,但生成的图表加载速度非常慢,因为点太多了。我尝试使用线性近似将点数减少到 10000(如果您的数据集之一大于另一个数据集,这就是 qqplot 函数所做的),但是您会丢失尾部的很多细节。

大多数朝向中心的数据点基本上是无用的——它们重叠太多,以至于每个像素可能有大约 100 个。是否有任何简单的方法可以删除过于靠近的数据,而不会将更稀疏的数据丢失到尾部?

4个回答

除了尾部之外,QQ 图具有难以置信的自相关性。在审查它们时,人们关注情节的整体形状和尾部行为。 因此,您可以通过在分布中心进行粗略二次抽样并包括足够数量的尾部来做得很好。

这是说明如何在整个数据集上采样以及如何取极值的代码。

quant.subsample <- function(y, m=100, e=1) {
  # m: size of a systematic sample
  # e: number of extreme values at either end to use
  x <- sort(y)
  n <- length(x)
  quants <- (1 + sin(1:m / (m+1) * pi - pi/2))/2
  sort(c(x[1:e], quantile(x, probs=quants), x[(n+1-e):n]))
  # Returns m + 2*e sorted values from the EDF of y
}

为了说明这一点,这个模拟数据集显示了两个大约 120 万个值的数据集之间的结构差异,以及其中一个数据集中的极少量“污染”。此外,为了使这个测试更严格,从其中一个数据集中完全排除了一个值区间:QQ 图需要显示这些值的中断。

set.seed(17)
n.x <- 1.21 * 10^6
n.y <- 1.20 * 10^6
k <- floor(0.0001*n.x)
x <- c(rnorm(n.x-k), rnorm(k, mean=2, sd=2))
x <- x[x <= -3 | x >= -2.5]
y <- rbeta(n.y, 10,13)

我们可以对每个数据集的 0.1% 进行二次抽样,并包括另外 0.1% 的极端值,从而绘制 2420 个点。总经过时间小于 0.5 秒:

m <- .001 * max(n.x, n.y)
e <- floor(0.0005 * max(n.x, n.y))

system.time(
  plot(quant.subsample(x, m, e), 
       quant.subsample(y, m, e), 
       pch=".", cex=4,
       xlab="x", ylab="y", main="QQ Plot")
  )

不会丢失任何信息:

QQ情节

在这个线程的其他地方,我提出了一个简单但有点特别的对点进行二次采样的解决方案。它速度很快,但需要一些实验才能产生出色的情节。即将描述的解决方案要慢一个数量级(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情节

编辑

我已经修改了原始代码,qq以便将第三列索引返回到原始两个数组中最长(或最短,如指定)xy,对应于所选点。这些索引指向数据的“有趣”值,因此可用于进一步分析。

我还删除了一个重复值的错误x(导致beta未定义)。

删除中间的一些数据点会改变经验分布,从而改变 qqplot。话虽如此,您可以执行以下操作并直接绘制经验分布的分位数与理论分布的分位数:

x <- rnorm(1200000)
mean.x <- mean(x)
sd.x <- sd(x)
quantiles.x <- quantile(x, probs = seq(0,1,b=0.000001))
quantiles.empirical <- qnorm(seq(0,1,by=0.000001),mean.x,sd.x)
plot(quantiles.x~quantiles.empirical) 

您必须根据要进入尾部的深度来调整 seq。如果你想变得聪明,你也可以在中间细化那个序列来加快情节。例如使用

plogis(seq(-17,17,by=.1))

是一种可能。

你可以做一个hexbin情节。

x <- rnorm(1200000)
mean.x <- mean(x)
sd.x <- sd(x)
quantiles.x <- quantile(x, probs = seq(0,1,b=0.000001))
quantiles.empirical <- qnorm(seq(0,1,by=0.000001),mean.x,sd.x)

library(hexbin)
bin <- hexbin(quantiles.empirical[-c(1,length(quantiles.empirical))],quantiles.x[-c(1,length(quantiles.x))],xbins=100)
plot(bin)