Gaussian Mixture Model parameters from density

Question

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

Answer 1

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)
curve(dmix,add=T,col="red")

Answer 2

If you are willing to code a bit, you can implement your own version of the EM algorithm that takes into account the density value of a grid point. So instead of using the usual likelihood function:

$$ \log L(\Theta) = \sum_{t=1}^L\left [{\sum_{i=1}^N \phi(\boldsymbol r_t|\boldsymbol \mu_i,\boldsymbol \Sigma_i)} \right] $$

where the outer sum is over the data points and the inner sum is over the GMM. You would instead want a density-weighted likelihood function:

$$ \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] $$

This would also affect the M-step - you'd have to multiply the responsibilities (that is, the probability of each data point) by the density at that point, and make sure things are still normalized correctly. Something like this was done in this paper:

http://www.ncbi.nlm.nih.gov/pubmed/18708469