Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm using Mathematica 9.0 to fit multidimensional observations with a gaussian mixture. To do so I first adapted these steps to ensure that the mixture parameters are well defined (see Runing time was unacceptable (I waited some hours and aborted it afterwards), so I tried the fitting in Matlab. Matlab instantly retrieved an error message from the clustering algorithm that states "ill-conditioned covariance created at iteration 3". To overcome this I tried 2 solutions concerning the covariance matrices definitions following the suggestions in and

To reproduce my problem and illustrate my experiments here I enclose some data, part of my original data, this being 7 samples each 5-dimensional, which I try to fit in a 2 gaussian mixture:

data = {{0, 2.76, 10.68, 12.4133, 9.26667}, {0, 0, 0, 0, 1.13333}, {6.06667, 8.32, 6.86667, 6.6, 4.49333}, {1.25333, 5.96, 7.49333, 5.81333, 5.05333}, {0, 0, 1.13333, 6.42667, 7.41333}, {0, 0, 0, 3.25333, 3.82667}, {0.76, 6.28, 10.28, 9.05333, 5.12}};
data // MatrixForm
d = Dimensions[data]
nDim = Last[d];
nMixtures = 2;

One solution is to use a single covariance matrix for all the mixtures. What follows is the implementation in Mathematica:

For[k = 1, k < nMixtures + 1, k++, 
    m[k] = Table[mm[k, i], {i, nDim}]];(*means*)
For[k = 1, k < nMixtures, k++, 
    w[k] = Exp[ww[k]]/((1 + Exp[ww[k]])*(nMixtures - 1))]; (*mixture weights*)
c = Table[If[i >= j, cc[i, j], 0], {i, nDim}, {j, nDim}];(*triangular matrix*)
s = c.Transpose[c];(*covariance matrix*)
M = Append[Array[w, nMixtures - 1], 1 - Total[Array[w, nMixtures - 1]]];(*all weights, which sum up to 1*)
Q = {};
For[i = 1, i < nMixtures + 1, i++, 
    Q = Insert[Q, MultinormalDistribution[m[i], s], i]];(*all mixtures, which share the same covariance matrix*)

The corresponding fitting takes 85 seconds in my PC using the following command:

edf_shared = EstimatedDistribution[data, MixtureDistribution[M, Q]]; // Timing

The same implementation in Matlab takes 0.0047 seconds:

data = [0,0,6.0667,1.2533,0,0,0.76; 2.76,0,8.32,5.96,0,0,6.28; 10.68,0,6.8667,7.4933,1.1333,0,10.28; 12.4133,0,6.6,5.8133,6.4267,3.2533,9.0533; 9.2667,1.1333,4.4933,5.0533,7.4133,3.8267,5.1200]';
edf_shared =,nMixtures,'SharedCov',true);

Another way to go around this problem is, instead of using the same covariance matrix for all gaussians, add a small epsilon to the diagonal of the covariance matrices. In Mathematica it takes 29,247 seconds (about 8 hours):

For[k = 1, k < nMixtures + 1, k++, 
    m[k] = Table[mm[k, i], {i, nDim}]];(*means*)
For[k = 1, k < nMixtures, k++, 
    w[k] = Exp[ww[k]]/((1 + Exp[ww[k]])*(nMixtures - 1))];(*mixture weigths*)
For[k = 1, k < nMixtures + 1, k++, 
    c[k] = Table[If[i >= j, cc[k, i, j], 0], {i, nDim}, {j, nDim}]];(*triangular matrices*)
epsilon = 0.00005;
epsilonMatrix = DiagonalMatrix[ConstantArray[epsilon, nDim]];(*diagonal matrix*)
For[k = 1, k < nMixtures + 1, k++, 
    s[k] = c[k].Transpose[c[k]] + epsilonMatrix];(*covariance matrices*)
M= Append[Array[w, nMixtures - 1], 1 - Total[Array[w, nMixtures - 1]]];
Q = {};
For[i = 1, i < nMixtures + 1, i++, 
    Q = Insert[Q, MultinormalDistribution[m[i], s[i]], i]];(*all mixtures, with regularized covariance matrices*)
edf_epsilon = EstimatedDistribution[data,MixtureDistribution[M, Q]]; // Timing

The corresponding implementation in Matlab takes 0.0055 seconds:

edf_epsilon =,nMixtures,'Regularize', 1e-5);

Neadless to say, when evaluating the estimated distribution function on a point, execution time is also an issue. The implementation using shared covariances in Mathematica requires 276 seconds, while in Matlab is just 0.41 seconds. The first approach involving a shared matrix, when applied to this new point takes 107 seconds in Mathematica and 0.20 seconds in Matlab.

I also tried playing with the EstimatedDistribution function parameters such as WorkingPrecision, AccuracyGoal, ParameterEstimator, etc., I didn't see any substantial improvement. At most, when I reduced the WorkingPrecision to 5 I managed to cut down the estimation time to just 40 seconds.

My question is whether there is a better way to do this, I mean to estimate the parameters of a gaussian distribution in Mathematica while coping for ill-conditioned covariance matrices. By better I mean faster while keeping a relatively good accuracy in the estimations, in the Matlab direction I enclosed here. Perhaps there is some issue with my code, or the approaches themselves are intrinsically inconvenient in Mathematica.

Any suggestion will be welcome, thanks in advance!

share|improve this question

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.