Obteniendo densidad de probabilidad de datos

Question

Necesito analizar algunos datos sobre sesiones de internet para una línea DSL. Quería echar un vistazo a cómo se distribuyen las duraciones de la sesión. Pensé que una forma simple de hacer esto sería comenzar haciendo una gráfica de densidad de probabilidad de la duración de todas las sesiones.

He cargado los datos en R y he usado la función density(). Por lo tanto, era algo como esto

plot(density(data$duration), type = "l", col = "blue", main = "Density Plot of Duration", 
    xlab = "duration(h)", ylab = "probability density") 

Soy nuevo en R y este tipo de análisis. Esto fue lo que encontré al pasar por google. Conseguí un plan, pero me dejaron algunas preguntas. ¿Es esta la función correcta para hacer lo que estoy tratando de hacer o hay algo más?

En la gráfica encontré que la escala del eje Y era de 0 ... 1.5. No entiendo cómo puede ser 1.5, ¿no debería ser de 0 ... 1?

Además, me gustaría obtener una curva más suave. Dado que el conjunto de datos es realmente grande, las líneas son realmente irregulares. Sería mejor tenerlos suavizados cuando presente esto. ¿Cómo voy a hacer eso?

Answer 1

Debe jugar con el parámetro de ancho de banda (bw) para cambiar la suavidad de la curva. En general, R hace un buen trabajo y automáticamente le da una curva agradable y suave, pero tal vez ese no sea el caso para su conjunto de datos específico.

En cuanto a la llamada que está utilizando, sí, es correcto, type="l" no es necesario, es el valor predeterminado utilizado para trazar objetos de densidad. El área debajo de la curva (es decir, la integral de -Inf a + Inf de su función de densidad) será = 1.

Ahora, ¿es la curva de densidad lo mejor para usar en su caso? Quizás, quizás no ... realmente depende del tipo de análisis que quieras hacer. Probablemente usar hist sea suficiente, y quizás aún más informativo, ya que puede seleccionar contenedores específicos de duración (vea ?hist para obtener más información).

Answer 2

Como dijo nico, deberías echar un vistazo a hist, pero también puedes combinar los dos. Entonces podría llamar a la densidad con lines en su lugar. Ejemplo:

duration <- rpois(500, 10) # For duration data I assume Poisson distributed 
hist(duration, 
    probability = TRUE, # In stead of frequency 
    breaks = "FD",  # For more breaks than the default 
    col = "darkslategray4", border = "seashell3") 
lines(density(duration - 0.5), # Add the kernel density estimate (-.5 fix for the bins) 
    col = "firebrick2", lwd = 3) 

deben darle algo como:

Tenga en cuenta que la estimación de la densidad del núcleo asume un kernel de Gauss como predeterminado. Pero el ancho de banda es a menudo el factor más importante. Si llama density directamente reporta el ancho de banda estimado por defecto:

> density(duration) 

Call: 
     density.default(x = duration) 

Data: duration (500 obs.);  Bandwidth 'bw' = 0.7752 

     x     y    
Min. : 0.6745 Min. :1.160e-05 
1st Qu.: 7.0872 1st Qu.:1.038e-03 
Median :13.5000 Median :1.932e-02 
Mean :13.5000 Mean :3.895e-02 
3rd Qu.:19.9128 3rd Qu.:7.521e-02 
Max. :26.3255 Max. :1.164e-01 

Aquí es 0.7752. Verifique sus datos y juegue con ellos como lo sugirió nico. Es posible que desee mirar ?bw.nrd.

Answer 3

Iba a agregar esto como un comentario a la respuesta anterior, pero es demasiado grande. El sesgo aparente se debe a la forma en que se agrupan los valores en un histograma. A menudo es un error usar histogramas para datos discretos. Consulte a continuación ...

set.seed(1001) 
tmpf <- function() { 
    duration <- rpois(500, 10) # For duration data I assume Poisson distributed 
    hist(duration, 
     probability = TRUE, # In stead of frequency 
     breaks = "FD",  # For more breaks than the default 
     col = "darkslategray4", border = "seashell3", 
     main="",ann=FALSE,axes=FALSE,xlim=c(0,25),ylim=c(0,0.15)) 
    box() 
    lines(density(duration), # Add the kernel density estimate 
     col = "firebrick2", lwd = 3) 
    par(new=TRUE) 
    plot(table(factor(duration,levels=0:25))/length(duration), 
     xlim=c(0,25),ylim=c(0,0.15),col=4,ann=FALSE,axes=FALSE) 
} 

par(mfrow=c(3,3),mar=rep(0,4)) 
replicate(9,tmpf())