function MixtureModelExample()
%
% This script generates some data from two different Gaussians and then
% combines the data into one big vector. It then fits a mixture model of
% two Gaussians to the data to try to recover the original Gaussians that
% generated the data (it uses the matlab function mle() to get the maximum
% likelihood mixture).
%
% Timothy F. Brady, tfbrady@mit.edu, Oct. 2008
%

    % Generate some data drawn from two Gaussians
    data = [0.4+randn(100,1).*0.15; 1+ randn(200,1).*0.25]';
    data(data < 0.05) = 0.05;
    [n,x] = hist(data);
    bar(x,n);
    
    % Make the mixture model pdf
    mixtureGauss = ...
        @(x,m1,s1,m2,s2,theta) (theta*normpdf(x,m1,s1) + (1-theta)*normpdf(x,m2,s2));
		
    % Set up parameters for the MLE function
    options = statset('mlecustom');
    options.MaxIter     = 20000;
    options.MaxFunEvals = 20000;

    % Get max likilihood parameters for our mixture model (start with some
    %   reasonable guesses about the parameters)
    p = mle(data, 'pdf', mixtureGauss, 'start', [0.5 0.1 0.5 0.1 0.5], ...
        'lowerbound', [-Inf 0 -Inf 0 0], 'upperbound', [Inf Inf Inf Inf 1], ...
        'options', options);

    % Plot and print information
    hold on;
    x = linspace(min(data),max(data),100);
    plot(x, mixtureGauss(x,p(1),p(2),p(3),p(4),p(5))*max(n), 'r', 'LineWidth', 2);
    fprintf('Gauss 1: %0.2f (+/- %0.2f)\n', p(1), p(2));
    fprintf('Gauss 2: %0.2f (+/- %0.2f)\n', p(3), p(4));
    fprintf('Mix: %0.2f proportion first gaussian\n', p(5));