Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 http://mathematica.stackexchange.com/a/16833/12672). 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 http://stackoverflow.com/a/16263971 and http://www.mathworks.com/matlabcentral/answers/100210.

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 = gmdistribution.fit(data,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 = gmdistribution.fit(data,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
add comment

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.