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.
ren
( x , y)ty
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).
approx()funkcja działa w tejqqplot()funkcji.