Draw density on top of histogram based on fitted GMM in R

Question

I want to simply draw the density on top of my histogram plot, using the means, and variance estimated using a GMM. I've been trying to do it, but I've been unable to draw the densities. The y-axis are always different.

This would be a toy exameple:

Data x coming from two normal distributions:

setseed(0)    
x1 <- rnorm(100,5,1)
x2 <- rnorm(100,10,1)
x <- c(x1,x2)
hist(x)

I then fit a GMM using the mclust package:

require(mclust)
gmm <- Mclust(x)
summary(gmm)

The two means, and (equal) variance for the two gaussians are:

gmm$parameters$mean ## 5.001579 and 9.931690 
gmm$parameters$variance$sigmasq ## 0.8516606

I can draw a histogram with different colors for the two normals based on the classification value outputted by the gmm. But how can I simply add two densities for each gaussian on top of this plot?

hist(x,breaks = seq(1,15,by=1),col="grey")
hist(x[gmm$classification==1],breaks = seq(1,15,by=1),col="red",add=T)
hist(x[gmm$classification==2],breaks = seq(1,15,by=1),col="blue",add=T)

Answer 1

There's a few assumptions in here, but I'll give it a try. First of all, I don't think you can easily do this with the standard hist and it likely needs ggplot2.

#libraries
library(ggplot2)
library(mclust)

#Creating your sample data
setseed(0)    
x1 <- rnorm(100,5,1)
x2 <- rnorm(100,10,1)
x <- c(x1,x2)
#Putting it in a dataframe for ggplot
df <- as.data.frame(x)

gmm <- Mclust(x)

gmm$parameters$mean ## 5.001579 and 9.931690 
gmm$parameters$variance$sigmasq ## 0.8516606

#Calculating the breaks hist() would use
brx <- pretty(range(df$x), 
              n = nclass.Sturges(df$x),min.n = 1)

#Adding the classification to the dataframe for the colors.
df$classification <- as.factor(x[gmm$classification])

#Plotting the histograms, adding the density (scaled * 80) and adding a 2nd y-axis to show that scale
ggplot(df, aes(x, fill= classification)) + 
  geom_histogram(col="grey", breaks=brx, alpha = 0.5) +
  geom_density(aes(y = 80 * ..density.. , col=classification, fill = NULL), size = 1) +
  scale_y_continuous(sec.axis = sec_axis(~./80))