library('MASS')
library('mixtools')
library('RColorBrewer')
library('dplyr')
library('ellipse')
palette <- brewer.pal(3, 'Set2')

######################
## Galaxies example ##
######################

data("galaxies")

## draw histogram of galaxies data
hist(galaxies, breaks=20, freq=FALSE)

## fit mixture model
fit <- normalmixEM(x=galaxies, k=3)

## plot the mixture component densities
x_plot <- seq(min(galaxies), max(galaxies), length.out=300)  ## x for plotting
d_plot <- 0*x_plot ## placeholder for mixture density

for(i in 1:length(fit$lambda)){

  ## plot component density
  d0_plot <- dnorm(x_plot, fit$mu[i], fit$sigma[i])
  lines(x_plot, fit$lambda[i]*d0_plot)
  
  ## add to mixture density
  d_plot <- d_plot + fit$lambda[i]*d0_plot
}

## plot mixture density
lines(x_plot, d_plot, lty=2)
legend('topright', c('component', 'mixture'), lty=c(1,2), bty='n')

## get cluster assignments
z <- apply(fit$posterior, 1, which.max)
points(x=galaxies, 
       y=rep(0,length(galaxies)), 
       col=palette[z],
       pch=20)

##################
## Iris example ##
##################

## plot iris data
data("iris")
plot(iris$Sepal.Length,
     iris$Sepal.Width,
     xlab='Sepal Length',
     ylab='Sepal Width',
     main='Iris (flower) Morphology',
     pch=20)

## fit mixture model
iris.x <- iris %>% 
  select(Sepal.Length, Sepal.Width) %>%
  as.matrix()
fit <- mvnormalmixEM(x=iris.x, k=2)

## get cluster assignments
z <- apply(fit$posterior, 1, which.max)

points(
  x=iris$Sepal.Length,
  y=iris$Sepal.Width,
  col=palette[z],
  pch=20)

for(i in 1:length(fit$lambda))
  for(a in seq(0.05, 0.95, 0.2))
    ellipse(mu=fit$mu[[i]], 
            sigma=fit$sigma[[i]], 
            alpha=a,
            col=palette[i])