¿Buenos métodos para los gráficos de densidad de variables no negativas en R?

Question

plot(density(rexp(100))

Obviamente, toda la densidad a la izquierda de cero representa un sesgo.

Quiero resumir algunos datos para los que no son estadísticos, y quiero evitar preguntas sobre por qué los datos no negativos tienen densidad a la izquierda de cero. Los gráficos son para comprobar la aleatoriedad; quiero mostrar las distribuciones de las variables por grupos de tratamiento y control. Las distribuciones suelen ser de tipo exponencial. Los histogramas son complicados por varias razones.

Una rápida búsqueda en Google me da trabajos de estadísticos sobre núcleos no negativos, por ejemplo: este .

¿Pero se ha implementado algo de esto en R? De los métodos implementados, ¿alguno de ellos es el "mejor" de alguna manera para la estadística descriptiva?

EDIT: aunque el from puede resolver mi problema actual, sería bueno saber si alguien ha implementado núcleos basados en la literatura sobre la estimación de la densidad no negativa

Answer 1

Una alternativa es el enfoque de Kooperberg y sus colegas, basado en la estimación de la densidad utilizando splines para aproximar la densidad logarítmica de los datos. Mostraré un ejemplo utilizando los datos de la respuesta de @whuber, que permitirá comparar los enfoques.

set.seed(17)
x <- rexp(1000)

Necesitarás el Logspline instalado para ello; instálelo si no lo está:

install.packages("logspline")

Cargue el paquete y calcule la densidad utilizando el logspline() función:

require("logspline")
m <- logspline(x)

En lo que sigue, asumo que el objeto d de la respuesta de @whuber está presente en el espacio de trabajo.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

El gráfico resultante se muestra a continuación, con la densidad logspline mostrada por la línea roja

Además, se puede especificar el soporte de la densidad mediante argumentos lbound y ubound . Si queremos suponer que la densidad es 0 a la izquierda de 0 y que hay una discontinuidad en 0, podríamos utilizar lbound = 0 en la llamada a logspline() por ejemplo

m2 <- logspline(x, lbound = 0)

El resultado es la siguiente estimación de la densidad (mostrada aquí con el original m logspline se ajustan como la figura anterior ya se estaba ocupando).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

El gráfico resultante se muestra a continuación

En este caso, aprovechar el conocimiento de x resulta en una estimación de la densidad que no tiende a 0 en $x = 0$ pero es similar al ajuste logspline estándar en otros lugares sobre x

Answer 2

Una solución, tomada de los enfoques de ponderación de bordes de las estadísticas espaciales, es truncar la densidad de la izquierda en cero, pero ponderar al alza los datos más cercanos a cero. La idea es que cada valor $x$ se "reparte" en un núcleo de área total unitaria centrado en $x$ Cualquier parte del núcleo que se extienda a territorio negativo se elimina y el núcleo se renormaliza a la unidad de área.

Por ejemplo, con un núcleo gaussiano $K_h(y,x) = \exp(-\frac{1}{2}((y-x)/h)^2) / \sqrt{2\pi}$ el peso de renormalización es

$$w(x) = 1 / \int_0^\infty K(y,x) dy = \frac{1}{1 - \Phi_{x, h}(0)}$$

donde $\Phi$ es la función de distribución acumulativa de una variante Normal de media $x$ y la desviación estándar $h$ . Existen fórmulas comparables para otros núcleos.

Esto es más sencillo -y mucho más rápido en el cálculo- que tratar de reducir los anchos de banda cerca de $0$ . Es difícil prescribir exactamente cómo deben cambiarse los anchos de banda cerca de $0$ De todos modos Sin embargo, este método también es ad hoc : todavía habrá algún sesgo cerca de $0$ . Parece que funciona mejor que la estimación de la densidad por defecto. Aquí hay una comparación utilizando un conjunto de datos de gran tamaño:

El azul muestra la densidad por defecto mientras que el rojo muestra la densidad ajustada para el borde en $0$ . La verdadera distribución subyacente está trazada como línea de puntos para referencia.

---

Código R

Le site density función en R se quejará de que la suma de pesos no es la unidad, porque quiere que la integral sobre todos los números reales sea la unidad, mientras que este enfoque hace que la integral sobre los números positivos sea igual a la unidad. Como comprobación, esta última integral se estima como una suma de Riemann.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)

Answer 3

Para comparar distribuciones por grupos (que dices que es el objetivo en uno de tus comentarios) ¿por qué no algo más sencillo? Los gráficos de caja paralelos funcionan bien si N es grande; los gráficos de franjas paralelas funcionan si N es pequeño (y ambos muestran bien los valores atípicos, que usted dice que es un problema en sus datos).

Answer 4

Como comenta Stéphane se puede utilizar from = 0 y, además, puede representar sus valores bajo la curva de densidad con rug (x)