Удаление посторонних точек возле центра QQ-графика

Я пытаюсь построить QQ-график с двумя наборами данных около 1,2 млн. Точек в R (используя qqplot и передавая данные в ggplot2). Вычисление достаточно простое, но полученный график очень медленно загружается, потому что точек очень много. Я пробовал линейное приближение, чтобы уменьшить количество точек до 10000 (это то, что делает функция qqplot в любом случае, если один из ваших наборов данных больше другого), но тогда вы теряете много деталей в хвостах.

Большинство точек данных к центру в основном бесполезны - они перекрывают друг друга настолько, что, вероятно, около 100 на пиксель. Есть ли какой-нибудь простой способ удаления данных, которые расположены слишком близко друг к другу, не теряя более разреженные данные по направлению к хвостам?

Сюжеты QQ невероятно автокоррелированы, за исключением хвостов. При их рассмотрении основное внимание уделяется общей форме сюжета и поведению хвоста. Ergo , у вас все будет хорошо, если вы будете грубо выполнять субсэмплинг в центрах распределений и включать достаточное количество хвостов.

Вот код, иллюстрирующий, как делать выборку по всему набору данных, а также как принимать экстремальные значения.

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
}

Чтобы проиллюстрировать, этот смоделированный набор данных показывает структурное различие между двумя наборами данных приблизительно в 1,2 миллиона значений, а также очень небольшое количество «загрязнения» в одном из них. Кроме того, чтобы сделать этот тест строгим, интервал значений исключается из одного из наборов данных в целом: график 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")
  )

Никакая информация не теряется вообще:

В другом месте этой темы я предложил простое, но несколько нерегулярное решение для подбора точек. Это быстро, но требует некоторых экспериментов, чтобы получить отличные сюжеты. Решение, которое должно быть описано, на порядок медленнее (до 1,2 секунды для 10 миллионов точек), но оно адаптивное и автоматическое. Для больших наборов данных он должен давать хорошие результаты с первого раза и делать это достаточно быстро.

$D_n$

$(x,y)$$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)

Давайте построим несколько версий, используя все меньшие и меньшие значения порога. При значении .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стала неопределенной).

Удаление некоторых точек данных в середине изменило бы эмпирическое распределение и, следовательно, 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)

Другой альтернативой является параллельный боксплот; вы сказали, что у вас есть два набора данных, так что-то вроде:

y <- rnorm(1200000)
x <- rnorm(1200000)
grpx <- cut(y,20)
boxplot(y~grpx)

и вы можете настроить различные параметры, чтобы сделать его лучше с вашими данными.