Gaussian Mixture Model parameters from density

How do I estimate parameters of subpopulations in a 1D gaussian mixture model when I already have density (measured on a grid) of the mixture?

All the algorithms I can find (like the well-known EM algorithm) assume that only samples from the mixture are available. My experiment directly yields density values, not separate samples. Of course I could sample manually from that density, but I feel there must be a better way

You can use minimum squared errors in order to estimate/fit a mixture density to your data (Note that this method also inherits problems of uniqueness of the estimators, as any other approach in the context of finite mixtures).

Basically, the idea is to minimize the distances between a mixture density (with fixed number of mixture components) and the observed density values, as a function of the parameters of the mixture you want to fit. The following R code shows an example with simulated data. The simulated data is obtained by simulating from a two-component Gaussian mixture and then using the command hist() (which emulates your context where only the density values are observed). As you can see, the estimators are very accurate in this example. The accuracy depends on how "informative" is the grid where you observe the density values.

rm(list=ls())
library(mixtools)

# Simulated data
set.seed(100)
n <- 500
lambda <- rep(1, 2)/2
mu <- c(0, 5)
sigma <- rep(1, 2)
mixnorm.data <- rnormmix(n, lambda, mu, sigma)
##A histogram of the simulated data.
hist(mixnorm.data,breaks=50)

# Binning the data
x.data <- hist(mixnorm.data,breaks=50,plot=FALSE)$mids den.data <- hist(mixnorm.data,breaks=50,plot=FALSE)$density

# Sum of squared errors
ld <- function(param){
mu1 = param[1]
mu2 = param[2]
sigma1 = param[3]
sigma2 = param[4]
eps = param[5]
if(sigma1>0 & sigma2>0 & eps >0 & eps<1){
return(sum((eps*dnorm(x.data,mu1,sigma1) + (1-eps)*dnorm(x.data,mu2,sigma2) - den.data)^2))
}
else return(Inf)
}

# Optimization step
optim(c(0,5,1,2,0.5),ld)

# Estimators of the parameters
MSEPAR <- optim(c(0,5,1,2,0.5),ld)\$par

# Fitted density
dmix <- Vectorize(function(x){
mu1 = MSEPAR[1]
mu2 = MSEPAR[2]
sigma1 = MSEPAR[3]
sigma2 = MSEPAR[4]
eps = MSEPAR[5]
return(  eps*dnorm(x,mu1,sigma1) + (1-eps)*dnorm(x,mu2,sigma2))
})

hist(mixnorm.data,breaks=50,probability=T)
$$\log L(\Theta) = \sum_{t=1}^L\left [{\sum_{i=1}^N \phi(\boldsymbol r_t|\boldsymbol \mu_i,\boldsymbol \Sigma_i)} \right]$$
$$\log L(\Theta) = \sum_{t=1}^Lp(\boldsymbol r_t)\left [{\sum_{i=1}^N \phi(\boldsymbol r_t|\boldsymbol \mu_i,\boldsymbol \Sigma_i)} \right]$$