Alguna sugerencia de cómo puedo trazar datos de tipo mixEM usando ggplot2

Tengo una muestra de 1 millón de registros obtenidos de mis datos originales. (Para su referencia, puede usar estos datos ficticios que pueden generar una distribución similar

b <- data.frame(matrix(rnorm(2000000, mean=c(8,17), sd=2))) c <- b[sample(nrow(b), 1000000), ] 

) Creí que el histogtwig era una mezcla de dos distribuciones logarítmicas normales e intenté ajustar las distribuciones sumdas usando el algoritmo EM usando el siguiente código:

 install.packages("mixtools") lib(mixtools) #line below returns EM output of type mixEM[] for mixture of normal distributions c1 <- normalmixEM(c, lambda=NULL, mu=NULL, sigma=NULL) plot(c1, density=TRUE) 

El primer gráfico es un gráfico de probabilidad de log y el segundo (si se presiona return nuevamente), da similar a las siguientes curvas de densidad:

Curvas de densidad modelo de mezcla

Como mencioné c1 es de tipo mixEM [] y la función plot () puede acomodar eso. Quiero llenar las curvas de densidad con colores. Esto es fácil de hacer usando ggplot2 () pero ggplot2 () no admite datos de tipo mixEM [] y arroja este mensaje:

“ggplot no sabe cómo tratar los datos de la clase mixEM” ¿Hay algún otro enfoque que pueda tomar para este problema? ¡¡Cualquier sugerencia es bienvenida!!

¡Gracias!

Mire la estructura del objeto devuelto (esto debe documentarse en la ayuda):

 > # simple mixture of normals: > x=c(rnorm(10000,8,2),rnorm(10000,17,4)) > xMix = normalmixEM(x, lambda=NULL, mu=NULL, sigma=NULL) 

Ahora que:

 > str(xMix) List of 9 $ x : num [1:20000] 6.18 9.92 9.07 8.84 9.93 ... $ lambda : num [1:2] 0.502 0.498 $ mu : num [1:2] 7.99 17.05 $ sigma : num [1:2] 2.03 4.02 $ loglik : num -59877 

Los componentes lambda, mu y sigma definen las densidades normales devueltas. Puede trazar estos en ggplot usando qplot y stat_function . Pero primero crea una función que devuelva densidades normales escaladas:

 sdnorm = function(x, mean=0, sd=1, lambda=1){lambda*dnorm(x, mean=mean, sd=sd)} 

Entonces:

 qplot(x,geom="density") + stat_function(fun=sdnorm,arg=list(mean=xMix$mu[1],sd=xMix$sigma[1], lambda=xMix$lambda[1]),fill="blue",geom="polygon") + stat_function(fun=sdnorm,arg=list(mean=xMix$mu[2],sd=xMix$sigma[2], lambda=xMix$lambda[2]),fill="#FF0000",geom="polygon") 

enter image description here

O cualquiera que ggplot habilidad de ggplot que tengas. Los colores transparentes en las densidades pueden ser agradables.

 ggplot(data.frame(x=x)) + geom_histogram(aes(x=x,y=..density..),fill="white",color="black") + stat_function(fun=sdnorm, arg=list(mean=xMix$mu[2], sd=xMix$sigma[2], lambda=xMix$lambda[2]), fill="#FF000080",geom="polygon") + stat_function(fun=sdnorm, arg=list(mean=xMix$mu[1], sd=xMix$sigma[1], lambda=xMix$lambda[1]), fill="#00FF0080",geom="polygon") 

productor:

enter image description here

Aquí hay un enfoque ligeramente diferente que usa geom_ploygon(...) lugar de múltiples llamadas a stat_function(...) . Un problema con stat_function(...) es que los argumentos secundarios (mu, sigma y lambda en este ejemplo), que se pasan usando el parámetro args=list(...) , no se pueden incluir en un mapeo estético, por lo que tienes que tener múltiples llamadas a stat_function(...) como es la solución @ Spacedman`s.

Este enfoque construye los archivos PDF fuera de ggplot y utiliza una única llamada a geom_polygon(...) . Como resultado, funciona sin modificación para un número arbitrario de distribuciones en la mezcla.

 # ggplot mixture plot gg.mixEM <- function(EM) { require(ggplot2) x <- with(EM,seq(min(x),max(x),len=1000)) pars <- with(EM,data.frame(comp=colnames(posterior), mu, sigma,lambda)) em.df <- data.frame(x=rep(x,each=nrow(pars)),pars) em.df$y <- with(em.df,lambda*dnorm(x,mean=mu,sd=sigma)) ggplot(data.frame(x=EM$x),aes(x,y=..density..)) + geom_histogram(fill=NA,color="black")+ geom_polygon(data=em.df,aes(x,y,fill=comp),color="grey50", alpha=0.5)+ scale_fill_discrete("Component\nMeans",labels=format(em.df$mu,digits=3))+ theme_bw() } library(mixtools) # two components set.seed(1) # for reproducible example b <- rnorm(2000000, mean=c(8,17), sd=2) c <- b[sample(length(b), 1000000) ] c2 <- normalmixEM(c, lambda=NULL, mu=NULL, sigma=NULL) gg.mixEM(c2) 

 # three components set.seed(1) b <- rnorm(2000000, mean=c(8,17,30), sd=c(2,3,5)) c <- b[sample(length(b), 1000000) ] library(mixtools) c3 <- normalmixEM(c, k=3, lambda=NULL, mu=NULL, sigma=NULL) gg.mixEM(c3) 

Intereting Posts