Нарисуйте плотность в верхней части гистограммы на основе установленного GMM в R

Question

Я хочу просто нарисовать плотность поверх графика гистограммы, используя средние значения и дисперсию, оцененную с использованием GMM. Я пытался сделать это, но я не смог нарисовать плотности. Ось Y всегда различна.

Это будет игрушечный пример:

Данные x поступают из двух нормальных распределений:

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

Затем я установил GMM, используя пакет mclust:

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

Два средних и (равных) дисперсии для двух гауссианцев:

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

Я могу нарисовать гистограмму с разными цветами для двух нормалей на основе значения classification, выведенного gmm. Но как я могу просто добавить две плотности для каждого гауссиана поверх этого графика?

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

Здесь есть несколько предположений, но я попробую.Прежде всего, я не думаю, что вы можете легко сделать это с помощью стандартного hist, и для него, вероятно, потребуется 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))