函数 gammamixEM（包 mixtools）的密度图形输出

Question

我正在使用包 mixtools 中的函数 gammamixEM。我如何 return 函数 normalmixEM 中密度的图形输出（即 plot(...,which=2) 中的第二个图）？

更新:

这里是 reproducible example 函数 gammamixEM:

x <- c(rgamma(200, shape = 0.2, scale = 14), rgamma(200, 
     shape = 32, scale = 10), rgamma(200, shape = 5, scale = 6))
out <- gammamixEM(x, lambda = c(1, 1, 1)/3, verb = TRUE)

这里是 reproducible example 函数 normalmixEM:

data(faithful)
attach(faithful)
out <- normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)
plot(out, which=2)

我想从函数 gammamixEM.

中获得密度的图形输出

Answer 1

给你。

out <- normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)

x          <- out 
whichplots <- 2
density = 2 %in% whichplots
loglik  = 1 %in% whichplots

def.par    <- par(ask=(loglik + density > 1), "mar") # only ask and mar are changed
mix.object <- x


k    <- ncol(mix.object$posterior)
x    <- sort(mix.object$x)
a    <- hist(x, plot = FALSE)
maxy <- max(max(a$density), .3989*mix.object$lambda/mix.object$sigma)

我不得不深入研究 plot.mixEM

的源代码

所以，现在用 gammamixEM 来做到这一点：

x <- c(rgamma(200, shape = 0.2, scale = 14), rgamma(200, 
                                                    shape = 32, scale = 10), rgamma(200, shape = 5, scale = 6))
gammamixEM.out <- gammamixEM(x, lambda = c(1, 1, 1)/3, verb = TRUE)



mix.object <- gammamixEM.out

k    <- ncol(mix.object$posterior)
x    <- sort(mix.object$x)
a    <- hist(x, plot = FALSE)
maxy <- max(max(a$density), .3989*mix.object$lambda/mix.object$sigma)

main2 <- "Density Curves"
xlab2 <- "Data" 
col2  <- 2:(k+1) 

hist(x, prob = TRUE, main = main2, xlab = xlab2, 
     ylim = c(0,maxy))


for (i in 1:k) {
  lines(x, mix.object$lambda[i] * 
          dnorm(x, 
                sd = sd(x)))
}

我相信如果你想添加标签、平滑线等，继续这个例子应该很简单。这是 plot.mixEM 函数的 source。

函数 gammamixEM（包 mixtools）的密度图形输出

Graphical output of density for the function gammamixEM (package mixtools)

r

gamma-distribution