I'm trying to wrap my head around mixture modelling, and I've come across a small matlab script that seems relevant. In order to familiarize myself with pymix, I've decided to try rewriting the matlab script in python. I appear to have successfully fitted/decomposed the data with the expectation maximization algorithm, but I'm not sure I understand the nature of the method's output.
My question is as follows:
- How can I replicate the second figure from the matlab code (below) with the output of
- Can someone explain the output of the above method?
I should mention that I'm horribly under-educated in statistics and mathematics in general (I'm working on it!), so please assume I know nothing =)
In any case, here's the matlab script:
% 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));
And here's what I've done in python so far. Note that I'm running this code in iPython with the
--pylab=inline option, thereby importing pyplot into my main workspace:
import numpy as np import mixture # Generate some data drawn from two Gaussians data = np.concatenate((0.4 + np.random.randn(100) * 0.15, 1 + np.random.randn(200) * 0.25)) data[np.nonzero(data < .05)] = .05 print type(data), len(data) plt = Figure() hist(data, bins=50) show() # Create DataSet object mixdat = mixture.DataSet() mixdat.fromArray(data) # reasonable-guess mixture, akin to random starting point in K-Means n1 = mixture.NormalDistribution(-2, 0.4) n2 = mixture.NormalDistribution(2, 0.6) # Mixture model and EM clustering mix = mixture.MixtureModel(2, [.5, .5], [n1, n2]) postmat, _ = mix.EM(mixdat, 40, 0.1) fig = Figure() hist(data, bins=50) x = np.linspace(np.min(data), np.max(data), 100) # Now what? show()
Any other comments, criticisms, or explanations are more than welcome. Thanks for putting up with my ignorance ;-)