rk5*_*567 7 r em ggplot2 density-plot mixture-model
我有一个从原始数据中获得的1m记录样本.(供您参考,您可以使用可能产生大致相似分布的虚拟数据
b <- data.frame(matrix(rnorm(2000000, mean=c(8,17), sd=2)))
c <- b[sample(nrow(b), 1000000), ]
Run Code Online (Sandbox Code Playgroud)
我认为直方图是两个对数正态分布的混合,我试图使用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)
Run Code Online (Sandbox Code Playgroud)
第一个图是对数似然图,第二个图(如果再次点击返回),给出类似于以下密度曲线:

正如我所提到的,c1的类型为mixEM [],而plot()函数可以容纳它.我想用颜色填充密度曲线.这很容易使用ggplot2()但ggplot2()不支持mixEM []类型的数据并抛出此消息:
"ggplot不知道如何处理类mixEM的数据"我还能采取其他方法解决这个问题吗?任何建议都非常感谢!!
谢谢!
查看返回对象的结构(这应该在帮助中记录):
> # simple mixture of normals:
> x=c(rnorm(10000,8,2),rnorm(10000,17,4))
> xMix = normalmixEM(x, lambda=NULL, mu=NULL, sigma=NULL)
Run Code Online (Sandbox Code Playgroud)
怎么办:
> 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
Run Code Online (Sandbox Code Playgroud)
lambda,mu和sigma组件定义返回的正常密度.您可以使用qplot和在ggplot中绘制这些图stat_function.但首先要创建一个返回缩放正常密度的函数:
sdnorm =
function(x, mean=0, sd=1, lambda=1){lambda*dnorm(x, mean=mean, sd=sd)}
Run Code Online (Sandbox Code Playgroud)
然后:
qplot(x,geom="density") + stat_function(fun=sdnorm,args=list(mean=xMix$mu[1],sd=xMix$sigma[1], lambda=xMix$lambda[1]),fill="blue",geom="polygon") + stat_function(fun=sdnorm,args=list(mean=xMix$mu[2],sd=xMix$sigma[2], lambda=xMix$lambda[2]),fill="#FF0000",geom="polygon")
Run Code Online (Sandbox Code Playgroud)

或者ggplot你拥有的任何技能.密度上的透明颜色可能很好.
ggplot(data.frame(x=x)) +
geom_histogram(aes(x=x,y=..density..),fill="white",color="black") +
stat_function(fun=sdnorm,
args=list(mean=xMix$mu[2],
sd=xMix$sigma[2],
lambda=xMix$lambda[2]),
fill="#FF000080",geom="polygon") +
stat_function(fun=sdnorm,
args=list(mean=xMix$mu[1],
sd=xMix$sigma[1],
lambda=xMix$lambda[1]),
fill="#00FF0080",geom="polygon")
Run Code Online (Sandbox Code Playgroud)
生产:

这是一种略有不同的方法,它使用geom_ploygon(...)而不是多次调用stat_function(...).一个问题stat_function(...)是使用args=list(...)参数传递的辅助参数(在此示例中为mu,sigma和lambda)不能包含在美学映射中,因此您必须多次调用stat_function(...)@ Spacedman的解决方案.
这种方法在ggplot之外构建PDF并使用单个调用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)
Run Code Online (Sandbox Code Playgroud)

# 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)
Run Code Online (Sandbox Code Playgroud)
