Applying Bayes: Estimating a bimodal distribution

Question

I'm trying to estimate a bunch of bimodal distributions, i.e. two means and two standard deviations, based on a variable number of inputs. If no input is present, a constant value should be returned. From my anecdotal knowledge of what Bayes is about, this should be exactly it, no? Adapt a prior based on incoming evidence. My question is, though, how would I go about that in detail? Maybe I'm just missing a good tutorial somehow.

Answer 1

EM algorithm

In response to @Daniel Johnson, I want to quickly show you how you can fit the EM algorithm in R. Use the package mixtools (click for a link). Then, you can use the normalmixEM function with the option k=2 to estimate the parameters of a two-component gaussian mixture distribution. Let's make an example:

# Let's look at the old faithful data

library(mixtools)

data(faithful)
plot(density(faithful$waiting), las=1, col="steelblue", lwd=2, main="")

The distribution is clearly bimodal, so we're gonna calculate a mixture model with 2 components (k=2):

out <- normalmixEM(faithful$waiting, k=2, epsilon = 1e-03, fast=TRUE)

summary(out)    

summary of normalmixEM object:
          comp 1    comp 2
lambda  0.361283  0.638717
mu     54.628096 80.099412
sigma   5.882584  5.859425
loglik at estimate:  -1034.002 

The first normal distribution has an estimated mean of $54.6$ with a standard deviation of $5.88$, the second a mean of $80.1$ with a standard deviation of $5.86$. This corresponds nicely to the peaks in the graph above.

---

Bayesian methods

To estimate the mixture models with Bayesian methodology, use the bayesmix package for R (click here for link). You'll need to install JAGS on your computer first, though (just download the exectuble file and install it). To illustrate its use, I re-run the above example of old faithful data. I choose independent (option independence) and uninformative priors (option priorsUncertain). Further, we run 10000 MCMC samples (option n.iter=10000) and discard the first 1000 as burn-in samples (option burn.in=1000):

library(bayesmix)

data(faithful)

bayesmod <- BMMmodel(faithful$waiting, k=2, priors=list(kind = "independence", parameter = "priorsUncertain", hierarchical = NULL)) # k=2 for two components

jcontrol <- JAGScontrol(variables = c("mu", "tau", "eta", "S"), burn.in = 1000, n.iter = 10000, seed = 10)

z <- JAGSrun(faithful$waiting, model = bayesmod, control = jcontrol, tmp = FALSE, cleanup = TRUE)

zSort <- Sort(z, by = "mu")

zSort

The model output is ("mu" denotes the estimates of the means and "sigma2" the estimates of the variances):

Markov Chain Monte Carlo (MCMC) output:
Start = 1001 
End = 11000 
Thinning interval = 1 

 Empirical mean, standard deviation and 95% CI for eta 
         Mean     SD   2.5%  97.5%
eta[1] 0.3622 0.0318 0.3013 0.4251
eta[2] 0.6378 0.0318 0.5749 0.6987

 Empirical mean, standard deviation and 95% CI for mu 
       Mean     SD  2.5% 97.5%
mu[1] 54.63 0.7365 53.25 56.12
mu[2] 80.08 0.5128 79.02 81.05

 Empirical mean, standard deviation and 95% CI for sigma2 
           Mean    SD  2.5% 97.5%
sigma2[1] 36.11 7.209 24.61 52.76
sigma2[2] 35.38 5.066 26.90 46.34

The estimated means are $54.63$ and $80.08$ with standard deviations (square roots of the variances "sigma2") of $\sqrt(36.11)\approx 6.01$ and $\sqrt(35.38)\approx 5.95$. This is very close to the estimates we've calculated with normalmixEM.

Answer 2

I think "A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models" is exactly what you're looking for. Particularly the third section: Finding Maximum Likelihood Mixture Densities Parameters via EM.