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.
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 (
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.
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
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 0.3622 0.0318 0.3013 0.4251 eta 0.6378 0.0318 0.5749 0.6987 Empirical mean, standard deviation and 95% CI for mu Mean SD 2.5% 97.5% mu 54.63 0.7365 53.25 56.12 mu 80.08 0.5128 79.02 81.05 Empirical mean, standard deviation and 95% CI for sigma2 Mean SD 2.5% 97.5% sigma2 36.11 7.209 24.61 52.76 sigma2 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
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.