Q-Q-plots er utroligt autokorreleret undtagen i halerne. Når man gennemgår dem, fokuserer man på plotens overordnede form og på halenes adfærd. Ergo , du klarer dig fint ved groft undersampling i centrene for distributionerne og med en tilstrækkelig mængde haler.

Her er kode, der illustrerer, hvordan at prøve på tværs af et helt datasæt samt hvordan man tager ekstreme værdier.

  quant.subsample <- function (y, m = 100, e = 1) {# m: størrelse på en systematisk prøve nr. e: antal ekstreme værdier i hver ende for at bruge x <- sort (y) n <- længde (x) kvantiteter <- (1 + sin (1: m / (m + 1) * pi - pi / 2 )) / 2 sortering (c (x [1: e], kvantil (x, probs = kvanter), x [(n + 1-e): n])) # Returnerer m + 2 * e sorterede værdier fra EDF af y}  

For at illustrere viser dette simulerede datasæt en strukturel forskel mellem to datasæt på ca. 1,2 millioner værdier samt en meget lille mængde "forurening" i et af dem. For at gøre denne test streng, er et interval af værdier udelukket fra et af datasættene helt og holdent: QQ-plottet skal vise et brud for disse værdier.

  set.seed (17) nx <- 1,21 * 10 ^ 6n.y <- 1,20 * 10 ^ 6k <- etage (0,0001 * nx) x <- c (rnorm (nx-k), rnorm (k, middel = 2, sd = 2)) x <- x [x < = -3 | x > = -2.5] y <- rbeta (ny, 10,13)  

Vi kan underprøve 0,1% af hvert datasæt og inkludere yderligere 0,1% af deres ekstremer, hvilket giver 2420 point til plot . Samlet forløbet tid er mindre end 0,5 sekunder:

  m <- .001 * max (nx, ny) e <- floor (0,0005 * max (nx, ny)) system.time (plot) (quant.subsample (x, m, e), quant.subsample (y, m, e), pch = ".", cex = 4, xlab = "x", ylab = "y", main = "QQ Plot "))  

Ingen information overhovedet går tabt:

Andetsteds i denne tråd foreslog jeg en enkel, men noget ad hoc løsning til undersampling af punkterne. Det er hurtigt, men kræver nogle eksperimenter for at producere store plot. Løsningen, der skal beskrives, er en størrelsesorden langsommere (tager op til 10 sekunder for 1,2 millioner point), men er adaptiv og automatisk. For store datasæt burde det give gode resultater første gang og gøre det rimeligt hurtigt.

Ideen er, at Douglas-Peucker polyline-forenklingsalgoritmen er tilpasset egenskaberne af et QQ-plot. Den relevante statistik for et sådant plot er Kolmogorov-Smirnov-statistikken $ D_n $, den maksimale lodrette afvigelse fra en monteret linje. Følgelig er algoritmen denne:

Find den maksimale lodrette afvigelse mellem linjen, der forbinder ekstremen af ​​$ (x, y) $ -parret og deres QQ-plot. Hvis dette er inden for en acceptabel brøkdel $ t $ af hele spektret af $ y $, skal du udskifte plottet med denne linje. Ellers skal du opdele dataene i dem, der går forud for punktet med den maksimale lodrette afvigelse, og dem efter dem, og anvende algoritmen rekursivt til de to stykker.

Der er nogle detaljer at tage sig af, især til klare datasæt af forskellig længde. Jeg gør dette ved at erstatte den kortere med de kvantiler, der svarer til den længere: faktisk anvendes en stykkevis lineær tilnærmelse af EDF for den kortere i stedet for dens faktiske dataværdier. ("Kortere" og "længere" kan vendes ved at indstille use.shortest=TRUE.)

Her er en R implementering.

  qq <- funktion (x0, y0, ty = 0,0005, use.shortest = FALSE) {qq.int <- funktion (x, y, i.min, i.max) {# x, y er sorteret og af samme længde n <-længde (y) hvis (n == 1) returnerer (c (x = x, y = y, i = i.max)) hvis (n == 2) returnerer (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, hvilket.max (abs (y-fit)))) hvis ( abs (y [i] -fit [i]) > tærsk) {samle (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))} andet {cbind (x = c (x [1], x [n]), y = c (y [1], y [n]), i = c (i.min, i.max))}} samler <- funktion (xy1, xy2) {rbind (xy1, xy2 [-1,] )} # # Forbehandling af input, så sortering sker en gang #, og de fleste detaljer udvindes fra dataene. # is.omvendt <-længde (y0) <-længde (x0) hvis (use.shortest) er.omvendt <-! er.omvendt hvis (er.omvendt) {y <- sort (x0) n <- længde (y ) x <- kvantil (y0, prob = (1: n-1) / (n-1))} andet {y <- sort (y0) n <- længde (y) x <- kvantil (x0, prob = (1: n-1) / (n-1))} # # Konverter den relative tærskel ty til et absolut. # tærsk <- t.y * diff (interval (y)) # # Få rekursivt point på QQ-plottet. # xy <- qq.int (x, y, 1, n) hvis (er.omvendt) cbind (x = xy [, 2], y = xy [, 1], i = xy [, 3]) ellers xy }  

Som et eksempel bruger jeg data simuleret som i mit tidligere svar (med en ekstrem høj outlier smidt i y og en hel del mere forurening i x denne gang):

  set.seed (17) nx <- 1,21 * 10 ^ 6n.y <- 1,20 * 10 ^ 6k <- etage (0,01 * nx) x <- c (rnorm (nx-k), rnorm (k, middel = 2, sd = 2)) x <- x [x < = -3 | x > = -2,5] y <- c (rbeta (n.y, 10,13), 1)  

Lad os plotte flere versioner ved hjælp af mindre og mindre værdier af tærsklen. Med en værdi på .0005 og vises på en skærm med en højde på 1000 pixels, ville vi garantere en fejl på ikke mere end halvdelen af ​​en lodret pixel overalt på plottet. Dette vises i gråt (kun 522 punkter, sammenføjet med linjesegmenter); de grovere tilnærmelser er plottet oven på den: først i sort, derefter i rødt (de røde punkter vil være en delmængde af de sorte og overplotte dem), derefter i blå (som igen er en delmængde og overplot). Tidspunkterne spænder fra 6,5 ​​(blå) til 10 sekunder (grå). I betragtning af at de skalerer så godt, kan man lige så godt bruge ca. en halv pixel som en universel standard for tærsklen ( f.eks , 1/2000 til en 1000 pixel høj skærm) og gøres med det.

  qq.1 <- qq (x, y) plot (qq.1, type = "l", lwd = 1, col = "Grå", xlab = "x" , ylab = "y", main = "Adaptiv QQ-plot") punkter (qq.1, pch = ".", cex = 6, col = "Grå") punkter (qq (x, y, .01), pch = 23, col = "Sort") point (qq (x, y, .03), pch = 22, col = "Rød") point (qq (x, y, .1), pch = 19, col = " Blå ")  

Rediger

Jeg har ændret den oprindelige kode til qq for at returnere en tredje kolonne med indekser til den længste (eller korteste som specificeret) af de oprindelige to arrays, x og y , svarende til de valgte punkter. Disse indekser peger på "interessante" værdier af dataene og kan derfor være nyttige til yderligere analyse.

Jeg fjernede også en fejl, der forekommer med gentagne værdier på x (som forårsagede beta skal udefineres.

Fjernelse af nogle af datapunkterne i midten ville ændre den empiriske fordeling og derfor qqplot. Når det er sagt, kan du gøre følgende og direkte tegne kvantilerne af den empiriske fordeling vs. kvantilerne for den teoretiske fordeling:

  x <- rnorm (1200000) middelværdi. X <- middel (x) sd.x <- sd (x) kvantiler. x <- kvantil (x, probs = seq (0,1, b = 0,000001)) kvantiler. spiral <- qnorm (seq (0,1, ved = 0,000001) ), middel.x, sd.x) plot (quantiles.x ~ quantiles.empirical) 

Du bliver nødt til at justere seq afhængigt af hvor dybt du vil komme ind i halerne. Hvis du vil blive klog, kan du også tynde den rækkefølge i midten for at fremskynde plottet. For eksempel er brug af

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

en mulighed.

Du kunne lave et hexbin plot.

  x <- rnorm (1200000) middelværdi. x <- middelværdi (x) sd.x <- sd ( x) kvantiler. x <- kvantil (x, probs = seq (0,1, b = 0,000001)) kvantiler. empirisk <- qnorm (seq (0,1, ved = 0,000001), gennemsnit.x, sd.x) bibliotek (hexbin) bin <- hexbin (quantiles.empirical [-c (1, længde (quantiles.empirical))], quantiles.x [-c (1, længde (quantiles.x))], xbins = 100) plot (bin)  

Et andet alternativ er en parallel boxplot; du sagde, at du havde to datasæt, så noget som:

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

og du kan justere de forskellige muligheder for at gøre det bedre med dine data.