Любые предложения о том, как я могу построить данные типа mixEM с помощью ggplot2

У меня есть образец записей 1m, полученных из моих исходных данных. (Для справки, вы можете использовать эти фиктивные данные, которые могут генерировать приблизительно подобное распределение 

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

) Я полагал, что гистограмма представляет собой смесь двух логарифмически нормальных распределений, и попытался подогнать суммированные распределения с помощью алгоритма EM, используя следующий код: 

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)

Первый график является логарифмически вероятностным, а второй (если вы нажмете return еще раз), дает примерно следующую плотность кривые:

Кривые плотности модели смеси

Как я уже упоминал, c1 имеет тип mixEM [], и функция plot () может это учесть. Я хочу заполнить кривые плотности цветами. Это легко сделать с помощью ggplot2 (), но ggplot2() не поддерживает данные типа mixEM[] и выдает следующее сообщение:

"ggplot не знает, как обращаться с данными класса mixEM" есть ли какой-либо другой подход, который я могу принять для этой проблемы? Любые предложения очень ценятся!!

Спасибо!

2 6
ggplot2 density-plot

2 ответа:

Посмотрите на структуру возвращаемого объекта (это должно быть задокументировано в справке): 

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

А теперь что: 

> 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
Лямбда -, му-и Сигма-компоненты определяют возвращаемые нормальные плотности. Вы можете построить их в ggplot, используя qplot и stat_function. Но Сначала сделайте функцию, которая возвращает масштабированные нормальные плотности: 
sdnorm =
function(x, mean=0, sd=1, lambda=1){lambda*dnorm(x, mean=mean, sd=sd)}

Затем: 

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

Введите описание изображения здесь

Или какими бы то ни было ggplot навыками вы ни обладали. Прозрачные цвета на плотностях могут быть милый.

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

Производство:

Введите описание изображения здесь

8

Вот немного другой подход, который использует geom_ploygon(...) вместо нескольких вызовов stat_function(...). Одна из проблем с stat_function(...) заключается в том, что вторичные аргументы (mu, sigma и lambda в этом примере), которые передаются с помощью параметра args=list(...), не могут быть включены в эстетическое отображение, поэтому вам нужно иметь несколько вызовов stat_function(...), как решение @Spacedman.

Этот подход строит PDF-файлы вне ggplot и использует один вызов geom_polygon(...). В результате, он работает без изменений для произвольное число распределений в смеси. 
# 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)

5