mixDE <- function(dat, mu0, sd0, Np, Fr=0.8, Crp=0.8, Nit=100) {
   
   # Application of the Differential Evolution to the construction of the ML estimate
   # of the equally-weighted mixture of normal distributions
   # Note: this is just a visual demo of the DE algorithm
   #       for this problem, the EM algorithm could be better 
   # dat ... the vector of data (will be used to fit the mixture distribution)
   # mu0 ... vector of the initial estimates of means of the mixture components
   # sig0 ... vector of the initial estimates of standard devistions of the mixture components
   # Np ... the population size for the DE algorithm
   # Fr ... the differential weight of the DE algorithm
   # CRp ... the crossover probability of the DE algorithm 
   # Nit ... the number of iterations
   
   n <- length(dat) # the number of observations
   k <- length(mu0) # should be equal to length(sig0)
   m <- 2*k # the length of the feasible solutions
   
   f <- function(x) {
      # The objective function to be minimized
      # x ... a vector of size 2*k
      #       the first compnents k are the means
      #       the last k components are the standard deviations
      # Note: a negative standard deviation is replaced with a small eps>0
      
      d.vec <- rep(0, n)
      sd.vec <- pmax(x[(k + 1):m], rep(0.001, k))
      for (i.dat in 1:n)
         d.vec[i.dat] <- mean(dnorm(rep(dat[i.dat], k), 
                                mean = x[1:k], sd = sd.vec))
      - sum(log(d.vec))
   }
   
   # Create the initial population
   X <- matrix(0, nrow = Np, ncol = m)
   for (j in 1:k) {
      X[, j] <- rnorm(Np, mean = mu0[j])
      X[, j + k] <- rexp(Np, rate = 1/sd0[j])
   }

   xx <- seq(-4, 4, by = 0.1) # Just for the graphs
   
   for (it in 1:Nit) {
      print(it)
      for (i in 1:Np) {
         # Attempt to improve X[i,], i.e., the i-th member of the population
         
         # Create the donor vector y
         abc <- sample(setdiff(1:Np, i), 3, replace = FALSE)
         y <- X[abc[1],] + Fr*(X[abc[2],] - X[abc[3],])

         # Create the challenger z by the uniform recombination 
         # of X[i,] and the donor vector y 
         sel <- sample(0:1, m, replace = TRUE) # A small bug; should take Crp into account
         sel[sample(1:m, 1)] <- 1
         z <- (1 - sel)*X[i,] + sel*y
         
         # If X[i,] is not better than z, replace X[i,] by z
         if (f(X[i,]) >= f(z)) X[i,] <- z
      }
      
      # Plot the current pupulation as the corresponding densities
      if (it %% 5 == 0) {
         plot(xx, rep(0, length(xx)), type = "n", ylim = c(0, 1),
              xlab = "x", ylab = "density value", main = it)
         hist(dat, add = TRUE, freq = FALSE, col = "pink", breaks = 20)
         grid()
         points(dat, rep(0, length(dat)), pch = 19, cex = 0.1, col = "red")
         Dens <- matrix(0, ncol = length(xx), nrow = Np)
         for (i in 1:Np) {
            for (j in 1:k) {
               Dens[i,] <- Dens[i,] + dnorm(xx, mean = X[i, j],
                                            sd = max(X[i, j + k], 0.01))
            }
            lines(xx, Dens[i,]/k, col = "brown", lwd = 2)
         }
      }
   }
      
   return(X)
}

# Generate and scramble some artificial data
dat <- c(rnorm(50, mean = -1, sd = 0.2),
         rnorm(50, mean = 0, sd = 0.3),
         rnorm(50, mean = 1, sd = 0.4))
dat <- dat[sample(1:length(dat))]

# Run the DE algorithm with a very imprecise initial guess
mixDE(dat, c(-1.5, 0.5, 3), c(0.5, 0.9, 1.5), 60, Nit = 500)