Reputation: 31
I'm currently trying to plot the components found via EM algorithm. However, the estimated densities do not extend fully to the end. It looks like this:
My code is:
plot(EM_data, which=2, xlim= c(0, 80), xlab2= "", yaxt= "n", main2 ="", lwd2=0.8, border = "azure3")
lines(density(EM_data), lty=2, lwd=0.8)
The plot is truncated wether I specify xlim or not. xlim2 is not defined for this type of plot. Where am I going wrong?
Upvotes: 1
Views: 984
Reputation: 2047
The method to plot mixEM only draws within the range of the data, if you want to extend the densities you must build your own function.
Use something like this:
Example data:
library(mixtools)
data(faithful)
attach(faithful)
set.seed(100)
EM_data<-normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)
mixtools plot:
plot(EM_data, which=2, xlim= c(30, 110), xlab2= "", yaxt= "n", main2 ="",
lwd2=0.8, border = "azure3")
lines(density(EM_data$x), lty=2, lwd=0.8)
Adaptation by extending densities:
a <- hist(EM_data$x, plot = FALSE)
maxy <- max(max(a$density), 0.3989 * EM_data$lambda/EM_data$sigma)
hist(EM_data$x, prob = TRUE, main = "", xlab = "", xlim= c(30, 110),
ylim = c(0, maxy), yaxt= "n", border = "azure3")
for (i in 1:ncol(EM_data$posterior)) {
curve(EM_data$lambda[i] * dnorm(x, mean = EM_data$mu[i], sd = EM_data$sigma[i]),
col = 1 + i, lwd = 0.8, add = TRUE)
}
lines(density(EM_data$x), lty=2, lwd=0.8)
Upvotes: 2