Menghapus titik asing di dekat pusat plot QQ

Saya mencoba memplot QQ-plot dengan dua set data sekitar 1,2 juta poin, dalam R (menggunakan qqplot, dan memasukkan data ke ggplot2). Perhitungannya cukup mudah, tetapi grafik yang dihasilkan sangat lambat untuk dimuat, karena ada begitu banyak poin. Saya sudah mencoba pendekatan linier untuk mengurangi jumlah poin menjadi 10.000 (ini adalah fungsi qqplot, jika salah satu set data Anda lebih besar dari yang lain), tetapi kemudian Anda kehilangan banyak detail di bagian ekor.

Sebagian besar titik data ke pusat pada dasarnya tidak berguna - mereka tumpang tindih sehingga mungkin ada sekitar 100 per piksel. Apakah ada cara sederhana untuk menghapus data yang terlalu berdekatan, tanpa kehilangan data yang lebih jarang ke arah ekor?

Plot QQ sangat berkorelasi otomatis kecuali pada bagian ekor. Dalam mengulasnya, seseorang berfokus pada bentuk keseluruhan plot dan perilaku ekor. Ergo , kamu akan baik-baik saja dengan kasar subsampling di pusat-pusat distribusi dan termasuk jumlah yang cukup dari ekor.

Berikut ini adalah kode yang menggambarkan cara mengambil sampel di seluruh dataset serta cara mengambil nilai ekstrem.

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
}

Sebagai ilustrasi, dataset simulasi ini menunjukkan perbedaan struktural antara dua dataset dengan nilai sekitar 1,2 juta serta jumlah "kontaminasi" yang sangat kecil di salah satunya. Juga, untuk membuat pengujian ini ketat, interval nilai dikecualikan dari salah satu dataset secara bersamaan: plot QQ perlu menunjukkan jeda untuk nilai-nilai tersebut.

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)

Kami dapat melakukan subsampel 0,1% dari setiap dataset dan memasukkan 0,1% ekstrem lainnya, memberikan 2.420 poin untuk plot. Total waktu yang berlalu kurang dari 0,5 detik:

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

Tidak ada informasi yang hilang sama sekali:

Di tempat lain di utas ini saya mengusulkan solusi sederhana namun agak ad hoc untuk mengamplas poin. Ini cepat, tetapi membutuhkan beberapa percobaan untuk menghasilkan plot yang bagus. Solusi yang akan dijelaskan adalah urutan besarnya lebih lambat (mengambil hingga 10 detik untuk 1,2 juta poin) tetapi bersifat adaptif dan otomatis. Untuk kumpulan data besar, itu harus memberikan hasil yang baik pertama kali dan melakukannya dengan cepat.

$D_n$

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

Ada beberapa detail yang harus diperhatikan, terutama untuk mengatasi dataset dengan panjang yang berbeda. Saya melakukan ini dengan mengganti yang lebih pendek dengan kuantil yang sesuai dengan yang lebih panjang: pada dasarnya, pendekatan linear piecewise dari EDF yang lebih pendek digunakan sebagai pengganti nilai data aktualnya. ("Lebih pendek" dan "lebih lama" dapat dibalik dengan pengaturan use.shortest=TRUE.)

Berikut ini adalah Rimplementasinya.

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
}

Sebagai contoh, saya menggunakan data yang disimulasikan seperti pada jawaban saya sebelumnya (dengan pencilan yang sangat tinggi yang dilemparkan ke dalam ydan lebih banyak kontaminasi xsaat ini):

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)

Mari kita plot beberapa versi, menggunakan nilai ambang yang lebih kecil dan lebih kecil. Pada nilai .0005 dan ditampilkan pada monitor dengan tinggi 1000 piksel, kami akan menjamin kesalahan tidak lebih dari setengah piksel vertikal di mana-mana di plot. Ini ditampilkan dalam warna abu-abu (hanya 522 poin, bergabung dengan segmen garis); perkiraan kasar diplot di atasnya: pertama berwarna hitam, kemudian merah (titik merah akan menjadi subset dari yang hitam dan overplot), kemudian biru (yang lagi-lagi merupakan subset dan overplot). Rentang waktu mulai dari 6,5 (biru) hingga 10 detik (abu-abu). Mengingat bahwa mereka skala dengan sangat baik, orang mungkin menggunakan sekitar setengah-pixel sebagai standar universal untuk ambang ( misalnya , 1/2000 untuk monitor tinggi 1000-pixel) dan dilakukan dengan itu.

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

Edit

Saya telah memodifikasi kode asli untuk qqmengembalikan kolom indeks ketiga menjadi yang terpanjang (atau terpendek, sebagaimana ditentukan) dari dua array asli, xdan y, sesuai dengan poin yang dipilih. Indeks-indeks ini menunjuk ke nilai-nilai "menarik" dari data sehingga dapat berguna untuk analisis lebih lanjut.

Saya juga menghapus bug yang terjadi dengan nilai berulang x(yang menyebabkan betatidak terdefinisi).

Menghapus beberapa titik data di tengah akan mengubah distribusi empiris dan karenanya qqplot. Ini dikatakan, Anda dapat melakukan hal berikut dan langsung plot kuantil dari distribusi empiris vs kuantil dari distribusi teoritis:

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) 

Anda harus menyesuaikan seq tergantung pada seberapa dalam Anda ingin masuk ke dalam ekor. Jika Anda ingin menjadi pintar, Anda juga dapat mengencerkan urutan di tengah untuk mempercepat plot. Misalnya menggunakan

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

adalah suatu kemungkinan.

Anda bisa melakukan hexbinplot.

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)

Alternatif lain adalah boxplot paralel; Anda mengatakan Anda memiliki dua set data, jadi sesuatu seperti:

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

dan Anda bisa menyesuaikan berbagai opsi untuk membuatnya lebih baik dengan data Anda.