Usuwanie obcych punktów w pobliżu środka wykresu QQ

Próbuję wykreślić wykres QQ z dwoma zestawami danych około 1,2 miliona punktów, w R (używając qqplot i wprowadzając dane do ggplot2). Obliczenia są dość łatwe, ale wynikowy wykres jest boleśnie powolny do ładowania, ponieważ jest tak wiele punktów. Próbowałem aproksymacji liniowej, aby zmniejszyć liczbę punktów do 10000 (to właśnie robi funkcja qqplot, jeśli jeden z twoich zestawów danych jest większy od drugiego), ale wtedy tracisz wiele szczegółów w ogonach.

Większość punktów danych w kierunku centrum jest w zasadzie bezużyteczna - pokrywają się tak bardzo, że prawdopodobnie około 100 na piksel. Czy istnieje prosty sposób usuwania danych, które są zbyt blisko siebie, bez utraty bardziej rzadkich danych w kierunku ogonów?

Wykresy QQ są niezwykle autokorelowane, z wyjątkiem ogonów. Przeglądając je, skupiamy się na ogólnym kształcie fabuły i zachowaniu ogona. Ergo , poradzisz sobie z gruboziarnistym podpróbkowaniem w środkach dystrybucji i dołączeniem wystarczającej ilości ogonów.

Oto kod ilustrujący, jak próbkować w całym zestawie danych, a także jak przyjmować ekstremalne wartości.

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
}

Aby to zilustrować, ten symulowany zestaw danych pokazuje różnicę strukturalną między dwoma zestawami danych o wartości około 1,2 miliona wartości, a także bardzo niewielką ilość „zanieczyszczenia” w jednym z nich. Ponadto, aby ten test był bardziej rygorystyczny, przedział wartości jest całkowicie wykluczony z jednego z zestawów danych: wykres QQ musi pokazywać przerwę dla tych wartości.

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)

Możemy podpróbować 0,1% każdego zestawu danych i uwzględnić kolejne 0,1% ich skrajności, co daje 2420 punktów do wykreślenia. Całkowity czas, który upłynął, wynosi mniej niż 0,5 sekundy:

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")
  )

Żadne informacje nie zostaną utracone:

W innym miejscu tego wątku zaproponowałem proste, ale nieco ad hoc rozwiązanie podpróbkowania punktów. Jest szybki, ale wymaga eksperymentów, aby stworzyć świetne fabuły. Rozwiązanie, które ma zostać opisane, jest o rząd wielkości wolniejsze (zajmuje do 10 sekund dla 1,2 miliona punktów), ale jest adaptacyjne i automatyczne. W przypadku dużych zbiorów danych za pierwszym razem powinien dawać dobre wyniki i robić to dość szybko.

$D_n$

$(x,y)$$t$$y$

Jest kilka szczegółów, którymi należy się zająć, zwłaszcza w przypadku zestawów danych o różnej długości. Robię to, zastępując krótszy kwantylem odpowiadającym dłuższemu: w efekcie zamiast rzeczywistych wartości danych stosuje się częściowe przybliżenie liniowe EDF krótszego. („Krótsze” i „dłuższe” można odwrócić poprzez ustawienie use.shortest=TRUE).

Oto Rimplementacja.

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
}

Jako przykład wykorzystuję dane symulowane jak w mojej wcześniejszej odpowiedzi (z ekstremalnie wysoką wartością odstającą yi w xtym czasie znacznie większym zanieczyszczeniem ):

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)

Narysujmy kilka wersji, używając coraz mniejszych wartości progu. Przy wartości 0,0005 i wyświetlaniu na monitorze o wysokości 1000 pikseli gwarantowalibyśmy błąd nie większy niż połowa pionowego piksela wszędzie na wykresie. Jest to pokazane na szaro (tylko 522 punkty, połączone segmentami linii); grubsze aproksymacje są wykreślone na nim: najpierw na czarno, potem na czerwono (czerwone punkty będą podzbiorem czarnych i nadplanują je), a następnie na niebiesko (które znowu są podzbiorem i nadplotem). Zakresy czasowe wynoszą od 6,5 (niebieski) do 10 sekund (szary). Biorąc pod uwagę, że skalują się tak dobrze, równie dobrze można użyć około połowy piksela jako uniwersalnej wartości domyślnej dla progu ( np. 1/2000 dla monitora o wysokości 1000 pikseli) i można to zrobić.

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")

Edytować

Zmodyfikowałem oryginalny kod, qqaby zwracał trzecią kolumnę indeksów do najdłuższej (lub najkrótszej, jak określono) z dwóch oryginalnych tablic xi yodpowiadającej wybranym punktom. Indeksy te wskazują na „interesujące” wartości danych, a zatem mogą być przydatne do dalszej analizy.

Usunąłem również błąd występujący przy powtarzających się wartościach x(które spowodowały betaniezdefiniowanie).

Usunięcie niektórych punktów danych w środku zmieniłoby rozkład empiryczny, a zatem qqplot. Biorąc to pod uwagę, możesz wykonać następujące czynności i bezpośrednio wykreślić kwantyle rozkładu empirycznego w porównaniu z kwantylami rozkładu teoretycznego:

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) 

Będziesz musiał dostosować sekwencję w zależności od tego, jak głęboko chcesz dostać się do reszki. Jeśli chcesz być sprytny, możesz również rozrzedzić tę sekwencję na środku, aby przyspieszyć fabułę. Na przykład za pomocą

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

jest możliwość.

Możesz zrobić hexbinfabułę.

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)

Inną alternatywą jest równoległy wykres pudełkowy; powiedziałeś, że masz dwa zestawy danych, więc coś takiego:

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

i możesz dostosować różne opcje, aby poprawić swoje dane.