La eliminación de extraños puntos cerca del centro de un QQ-plot

Question

Estoy tratando de trazar un QQ-plot con dos conjuntos de datos de alrededor de 1,2 millones de puntos, en R (utilizando qqplot, y la alimentación de los datos en ggplot2). El cálculo es bastante fácil, pero el gráfico resultante es muy lento para cargar, porque hay muchos puntos. Lo he intentado en aproximación lineal para reducir el número de puntos a 10000 (esto es lo que el qqplot función de todos modos, si uno de los conjuntos de datos más grande que el otro), pero luego se pierde una gran cantidad de detalle en las colas.

La mayoría de los puntos de datos hacia el centro son básicamente inútiles - que se superponen de modo que mucho de lo que hay, probablemente, unos 100 por píxel. ¿Hay alguna forma más sencilla de extracción de datos que está demasiado cerca juntos, sin perder el más datos dispersos hacia las colas?

Answer 1

Q-Q parcelas son increíblemente autocorrelated excepto en las colas. En la revisión de los mismos, uno se centra en la forma general de la trama y en la cola de comportamiento. Ergo, estará bien por el grueso de submuestreo en los centros de las distribuciones, incluyendo una cantidad suficiente de las colas.

Aquí está el código que ilustra cómo la muestra a través de todo un conjunto de datos así como la forma de tomar los valores extremos.

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
}

Para ilustrar, este conjunto de datos simulados muestra una diferencia estructural entre dos conjuntos de datos de aproximadamente 1,2 millones de valores, así como una cantidad muy pequeña de "contaminación" en uno de ellos. Además, para hacer esta prueba rigurosa, de un intervalo de valores está excluido de uno de los conjuntos de datos por completo: el QQ plot necesita mostrar una pausa para esos valores.

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)

Podemos submuestra 0.1% de cada conjunto de datos e incluir otro 0.1% de sus extremos, dando 2420 puntos de la parcela. Total de tiempo transcurrido es de menos de 0.5 segundos:

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

Ninguna información se pierde en lo absoluto:

Answer 2

En otra parte de este hilo, he propuesto una simple pero algo ad hoc solución de submuestreo de los puntos. Es rápido, pero requiere un poco de experimentación para producir grandes parcelas. La solución a punto de ser descrito es de un orden de magnitud más lento (de hasta 10 segundos de 1,2 millones de puntos), pero es adaptable y automática. Para grandes conjuntos de datos, se debería dar buenos resultados en el primer tiempo y hacerlo con bastante rapidez.

La idea es que de la de Douglas-Peucker polilínea simplificación del algoritmo, adaptada a las características de un QQ plot. La estadística relevante para esa parcela es la prueba de Kolmogorov-Smirnov estadístico $D_n$, la máxima desviación vertical de una cocina equipada línea. En consecuencia, el algoritmo es este:

Encontrar el máximo de la desviación vertical entre la línea que une los extremos de la $(x,y)$ pares y sus QQ plot. Si esto está dentro de un aceptable fracción $t$ de la gama completa de $y$, reemplace la parcela con esta línea. De lo contrario, la partición de los datos en los anteriores el punto de máxima desviación vertical y los de después y de forma recursiva aplicar el algoritmo para las dos piezas.

Hay algunos detalles que cuidar, especialmente para hacer frente con conjuntos de datos de diferente longitud. Puedo hacer esto mediante la sustitución de la menor por los cuantiles correspondientes a más de uno: en efecto, una aproximación lineal por tramos de la FED de la menor se utiliza en lugar de sus valores reales de los datos. ("Corto" y "largo" puede ser revertido por la configuración de la use.shortest=TRUE.)

Aquí está una R de ejecución.

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
}

Como ejemplo puedo usar datos simulados como en mi anterior respuesta (con una alta extrema outlier lanzado en y y un poco más de contaminación en x este momento):

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)

Vamos parcela varias versiones, utilizando más y más pequeños los valores del umbral. En un valor de .0005 y mostrar en un monitor de 1000 píxeles de alto, estaríamos garantizando un error de no más de la mitad vertical de píxeles en todas partes en la parcela. Esto se muestra en gris (sólo 522 puntos, acompañado por segmentos de línea); la burdas aproximaciones se representan en la parte superior: la primera en negro, luego en color rojo (los puntos rojos será un subconjunto de los negros y overplot ellos), luego en azul (que de nuevo son un subconjunto y overplot). Los horarios rango de 6.5 (azul) a 10 segundos (en color gris). Dado que la escala tan bien, uno podría usar alrededor de la mitad de píxeles como un universal predeterminado para el umbral (por ejemplo, 1/2000 de 1000 píxeles de alto monitor) y hacer con ella.

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

Editar

He modificado el código original de qq a devolver una tercera columna de los índices en la más larga (o más corto, según se especifica) de la original de dos matrices, x y y, correspondiente a los puntos que están seleccionados. Estos índices señalan a "interesantes" los valores de los datos y por lo tanto podrían ser útiles para su posterior análisis.

También he eliminado un bug que ocurre con la repetición de los valores de x (lo que causó beta a ser indefinido).

Answer 3

La eliminación de algunos de los puntos de datos en el medio iba a cambiar la distribución empírica y, por tanto, la qqplot. Dicho esto, usted puede hacer lo siguiente, y directamente de la trama de los cuantiles de la distribución empírica frente a los cuantiles de la distribución teórica:

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) 

Usted tendrá que ajustar la seq dependiendo de qué tan profunda que usted desea conseguir en las colas. Si usted desea conseguir inteligente también puede delgada que la secuencia de en medio a la velocidad de la trama. Por ejemplo el uso de

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

es una posibilidad.

Answer 4

Usted podría hacer una hexbin de la parcela.

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)

Answer 5

Otra alternativa es un paralelo boxplot; usted dijo que había dos conjuntos de datos, algo así como:

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

y usted puede ajustar diversas opciones para hacerlo mejor con sus datos.