Ajustar una curva de densidad a un histograma en R

91

¿Existe una función en R que ajuste una curva a un histograma?

Digamos que tienes el siguiente histograma

hist(c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4)))

Parece normal, pero está sesgado. Quiero ajustar una curva normal que esté sesgada para envolver este histograma.

Esta pregunta es bastante básica, pero parece que no puedo encontrar la respuesta para R en Internet.

r histogram curve-fitting r-faq usuario5243421
fuente
¿Desea encontrar mys de modo que la distribución gaussiana N (m, s) se ajuste a sus datos?
SteinNorheim
No estoy seguro de lo que eso significa ...> _>
user5243421
10
@mathee: Creo que quiere decir m = media y s = desviación estándar. La distribución gaussiana es otro nombre para la distribución normal.
Peter Mortensen

Respuestas:

154

Si entiendo tu pregunta correctamente, probablemente quieras una estimación de densidad junto con el histograma:

X <- c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4))
hist(X, prob=TRUE)            # prob=TRUE for probabilities not counts
lines(density(X))             # add a density estimate with defaults
lines(density(X, adjust=2), lty="dotted")   # add another "smoother" density

Edite mucho tiempo después:

Aquí hay una versión un poco más elegante:

X <- c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4))
hist(X, prob=TRUE, col="grey")# prob=TRUE for probabilities not counts
lines(density(X), col="blue", lwd=2) # add a density estimate with defaults
lines(density(X, adjust=2), lty="dotted", col="darkgreen", lwd=2)

junto con el gráfico que produce:

ingrese la descripción de la imagen aquí

Dirk Eddelbuettel
fuente
3
+1: ¿puede hacerlo también al revés, es decir, ajustando el gráfico de densidad para que se ajuste al histograma?
vonjd
2
Sugiero dar un parámetro adicional a lines(density(X,na.rm= TRUE)ya que el vector puede contener valores NA.
Anirudh
30

Eso es fácil con ggplot2

library(ggplot2)
dataset <- data.frame(X = c(rep(65, times=5), rep(25, times=5), 
                            rep(35, times=10), rep(45, times=4)))
ggplot(dataset, aes(x = X)) + 
  geom_histogram(aes(y = ..density..)) + 
  geom_density()

o para imitar el resultado de la solución de Dirk

ggplot(dataset, aes(x = X)) + 
  geom_histogram(aes(y = ..density..), binwidth = 5) + 
  geom_density()
Thierry
fuente
28

Así es como lo hago:

foo <- rnorm(100, mean=1, sd=2)
hist(foo, prob=TRUE)
curve(dnorm(x, mean=mean(foo), sd=sd(foo)), add=TRUE)

Un ejercicio adicional es hacer esto con el paquete ggplot2 ...

John johnson
fuente
Sin embargo, si desea algo que esté sesgado, puede hacer el ejemplo de densidad anterior, transformar sus datos (por ejemplo, foo.log & lt; - log (foo) y probar lo anterior), o intentar ajustar una distribución sesgada, como el gamma o lognormal (lognormal es equivalente a tomar el registro y ajustar una normal, por cierto).
John Johnson
2
Pero eso aún requiere estimar los parámetros de su distribución primero.
Dirk Eddelbuettel
Esto se aleja un poco de la simple discusión de R, ya que nos adentramos más en las estadísticas teóricas, pero puede probar este enlace para Gamma: en.wikipedia.org/wiki/Gamma_distribution#Parameter_estimation Para lognormal, simplemente tome el registro (asumiendo todos los datos son positivos) y funcionan con datos transformados logarítmicamente. Para algo más elegante, creo que tendría que trabajar con un libro de texto de estadística.
John Johnson
3
Creo que malinterpreta cómo tanto el póster original como todas las demás respuestas se contentan con usar estimaciones no paramétricas, como un histograma de la vieja escuela o una estimación de densidad basada en datos algo más moderna. Las estimaciones paramétricas son excelentes si tiene buenas razones para sospechar una distribución. Pero ese no fue el caso aquí.
Dirk Eddelbuettel
11

Dirk ha explicado cómo trazar la función de densidad sobre el histograma. Pero a veces es posible que desee seguir el supuesto más sólido de una distribución normal sesgada y graficar eso en lugar de la densidad. Puede estimar los parámetros de la distribución y trazarlos usando el paquete sn :

> sn.mle(y=c(rep(65, times=5), rep(25, times=5), rep(35, times=10), rep(45, times=4)))
$call
sn.mle(y = c(rep(65, times = 5), rep(25, times = 5), rep(35, 
    times = 10), rep(45, times = 4)))

$cp
    mean     s.d. skewness 
41.46228 12.47892  0.99527

Gráfico de datos distribuidos sesgados-normales

Esto probablemente funcione mejor en datos que son más sesgados-normales:

Otra trama normal sesgada

fmark
fuente
3

Tuve el mismo problema, pero la solución de Dirk no pareció funcionar. Recibía este mensaje de advertencia cada vez

"prob" is not a graphical parameter

Leí ?histy encontré sobrefreq: a logical vector set TRUE by default.

el código que funcionó para mí es

hist(x,freq=FALSE)
lines(density(x),na.rm=TRUE)
Matias Andina
fuente