Qualche suggerimento su come posso tracciare dati di tipo mixEM usando ggplot2

Ho un campione di record di 1 milione ottenuto dai miei dati originali. (Per riferimento, è ansible utilizzare questi dati fittizi che possono generare una distribuzione approssimativamente simile

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

) Credevo che l’istogramma fosse una combinazione di due distribuzioni log-normali e ho provato ad adattare le distribuzioni sumte usando l’algoritmo EM usando il seguente codice:

 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) 

Il primo grafico è un grafico di probabilità di log e il secondo (se si preme di nuovo il ritorno), dà un risultato simile alle seguenti curve di densità:

Curve di densità del modello di miscela

Come ho detto c1 è di tipo mixEM [] e la funzione plot () può adattarsi a questo. Voglio riempire le curve di densità con i colors. Questo è facile da fare usando ggplot2 () ma ggplot2 () non supporta i dati di tipo mixEM [] e lancia questo messaggio:

“ggplot non sa come gestire i dati del class mixEM” C’è qualche altro approccio che posso prendere per questo problema? Qualsiasi suggerimento è molto apprezzato !!

Grazie!

Guarda la struttura dell’object restituito (questo dovrebbe essere documentato nella guida):

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

Ora cosa:

 > 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 

Le componenti lambda, mu e sigma definiscono le densità normali restituite. Puoi tracciare questi in ggplot usando qplot e stat_function . Ma prima crea una funzione che restituisca le densità normali ridimensionate:

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

Poi:

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

inserisci la descrizione dell'immagine qui

O qualunque ggplot abilità ggplot che hai. I colors trasparenti sulle densità potrebbero essere carini.

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

produzione:

inserisci la descrizione dell'immagine qui

Ecco un approccio leggermente diverso che utilizza geom_ploygon(...) invece di più chiamate a stat_function(...) . Un problema con stat_function(...) è che gli argomenti secondari (mu, sigma e lambda in questo esempio), che sono passati usando il parametro args=list(...) , non possono essere inclusi in una mapping estetica, quindi devi avere più chiamate a stat_function(...) come è la soluzione di @Spacedman.

Questo approccio costruisce i PDF al di fuori di ggplot e utilizza una singola chiamata a geom_polygon(...) . Di conseguenza, funziona senza modifiche per un numero arbitrario di distribuzioni nella miscela.

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