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:

  1. How can I replicate the second figure from the matlab code (below) with the output of pymix.MixtureModel.EM?
  2. 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);

% 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)

#  Create DataSet object
mixdat = mixture.DataSet()

# 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?

Any other comments, criticisms, or explanations are more than welcome. Thanks for putting up with my ignorance ;-)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.