set.seed(123)

N <- 1e5
mu <- c(2,4,6)
sigma <- rep(0.5, 3)
lambda <- (1:3)/6

dat <- c(
  rnorm(N*lambda[1], mu[1], sigma[1]),
  rnorm(N*lambda[2], mu[2], sigma[2]),
  rnorm(N*lambda[3], mu[3], sigma[3])
)

xx <- seq(0, 8, by = 0.01)
hist(dat, 100, freq = FALSE, col = "lightgray", border = 1, xlab = "", main = "")
lines(xx, lambda[1]*dnorm(xx, mu[1], sigma[1]), col = "red", lwd = 2)
lines(xx, lambda[2]*dnorm(xx, mu[2], sigma[2]), col = "blue", lwd = 2)
lines(xx, lambda[3]*dnorm(xx, mu[3], sigma[3]), col = "orange", lwd = 2)


fit <- mixtools::normalmixEM(dat, k = 3)
summary(fit)


f <- function(x, mu, sigma, lambda) {
  lambda[1]*dnorm(x, mu[1], sigma[1]) + lambda[2]*dnorm(x, mu[2], sigma[2]) + lambda[3]*dnorm(x, mu[3], sigma[3])
}

f(1:10,mu,sigma,lambda)

y <- seq(0, 8, length.out = 25)
points(y, f(y,mu,sigma,lambda), pch = 16)

plot(y, f(y,mu,sigma,lambda), pch = 16, type = "b")

fy <- f(y,mu,sigma,lambda)

optf <- function(par) {
  mu <- par[1:3]
  sigma <- par[4:6]
  lambda <- par[7:9]
  sum((fy - f(y, mu, sigma, lambda))^2)
}

tmp <- optim(c(1,1,1,1,1,1,0.33,0.33,0.33),
      optf,
      method = "L-BFGS-B",
      lower = c(-Inf, -Inf, -Inf, 0, 0, 0, 0, 0, 0),
      upper = c(Inf, Inf, Inf, Inf, Inf, Inf, 1, 1, 1))

plot(y, f(y,mu,sigma,lambda), pch = 16, type = "b", xlab="y", ylab="fy", axes =FALSE)
axis(1)
axis(2)

yy <- sample(y, 5000, prob = fy, replace = TRUE)
summary(mixtools::normalmixEM(yy, k = 3))