los - superponer graficas en r ggplot

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

Question

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 histograma era una mezcla de dos distribuciones logarítmicas normales e intenté ajustar las distribuciones sumadas 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:

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!

Answer 1

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)

Answer 2

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

O cualquier habilidad 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: