K = 3;
O = 2;
numex = 10;
data = zeros(numex, O);
tmu = zeros(O, K);
tSigma = zeros(O, O, K);
tmu(:,1) = [0 0]';
tmu(:,2) = [2 2]';
tSigma(:,:,1) = eye(2);
tSigma(:,:,2) = 2*eye(2);

data(1:5,:) = gsamp(tmu(:,1), tSigma(:,:,1), 5);
data(6:10,:) = gsamp(tmu(:,2), tSigma(:,:,2), 5);


mu = repmat(mean(data')', [1 K]);
Sigma = repmat(cov(data'), [1 1 K]);
mixcoef = normalise(rand(K,1)); % have to break symmetry somewhere!

max_iter = 5;

% Netlab
mix = gmm(O, K, 'full');
opt = foptions;
opt(14) = max_iter;
mix = gmminit(mix, data, opt);
mix0 = mix; % store the starting values
mix = gmmem(mix, data, opt);

% My code
[mu, Sigma, mixcoef] = learn_mog(data', K, 'full', max_iter, mix0.centres', mix0.covars,  mix0.priors');

assert(approxeq(mu, mix.centres'));
assert(approxeq(Sigma, mix.covars));
assert(approxeq(mixcoef, mix.priors'));

