Does Mathematica have support mixture models? Similar to the R language libraries mclust or mixtools. The documentation for mixtools is quite good.

The closest support I can find in Mathematica is the LearnDistribution[] function as documented here. I can get a PDF out of this, but I cant figure out how to get parameter estimates out, or to define the number of distributions it should determine, or define any start conditions.

  • $\begingroup$ MixtureDistribution + FindDistributionParameters? $\endgroup$
    – b3m2a1
    Nov 13 '20 at 15:37
  • 1
    $\begingroup$ This is all in the documention. Use Method->"GaussianMixture" and the sub-options here. For example: iris = ExampleData[{"MachineLearning", "FisherIris"}, "Data"][[All, 1, {1, 3}]]; ld = LearnDistribution[iris, Method -> {"GaussianMixture", "ComponentsNumber" -> 3}]. All the information about the parameters is in ld[[1]] if you need it - that's poking into the internals of Mathematica so this part is not documented. $\endgroup$
    – flinty
    Nov 13 '20 at 15:38
  • $\begingroup$ Well, goodness. How embarrassing. Thanks @flinty $\endgroup$
    – Q.P.
    Nov 13 '20 at 15:40
  • $\begingroup$ @Q.P. it is a little more complicated than that in fact - I will write up an answer. $\endgroup$
    – flinty
    Nov 13 '20 at 15:49
  • $\begingroup$ @flinty Yes I can see, I just ran some code. It's also quite different to some of the methods I've seen in the R libraries mclust and mixtools. $\endgroup$
    – Q.P.
    Nov 13 '20 at 16:16

Mathematica often pre-Standardize's data before learning. For example, see my answer here about learning SVM hyperplanes where I also encountered this problem. Therefore things are a bit more involved than my comment about the documentation if you want to extract the right parameters. I hope you can follow along:

(* some sample data, learn a 3 component mixture *) 
iris = ExampleData[{"MachineLearning", "FisherIris"}, "Data"][[All,1,{1,3}]];
ld = LearnDistribution[iris, 
   Method -> {"GaussianMixture", "ComponentsNumber" -> 3, 
     "CovarianceType" -> "Full"}];

(* extract the parameters from the internals *)
model = ld[[1]]["Model"];
{weights, means, choleskys} = 
  model /@ {"MixingCoefficients", "Means", "CholeskyCovariances"};

(* turn the cholesky covariances into +ve def symm matrices with U.U^T*)
covariances = #.ConjugateTranspose[#] & /@ choleskys;

(* create the mixture *)
distributions = MapThread[MultinormalDistribution, {means, covariances}];
mixture = MixtureDistribution[weights, distributions];

(* we first plot against the standardized points. Mathematica standardizes 
   all data before learning the parameters *) 
pdf = PDF[mixture, {x, y}];
DensityPlot[pdf, {x, -2, 3}, {y, -2, 3}, 
 Epilog -> {Red, Point[Standardize[iris]]}, PlotPoints -> 25]


In the above, we extracted the weights, means, and cholesky covariances, and we formed the symmetric covariances using $U.U^\top$ on each. Unfortunately these parameters are have been fit to Standardize[iris], not the original points iris. Therefore we need to apply the appropriate linear transformations to get parameters that fit the original data.

(* we need to unstandardize the parameters *)
sd = DiagonalMatrix@StandardDeviation[iris];
mu = Mean[iris];
finalmeans = (sd.# + mu) & /@ means;
finalcholeskys = (sd.#) & /@ choleskys;
finalcovariances = #.ConjugateTranspose[#] & /@ finalcholeskys;

(* create the un-standardized mixture and finally we can plot against 
  the original points *)
finaldistributions = MapThread[MultinormalDistribution, {finalmeans,finalcovariances}];
finalmixture = MixtureDistribution[weights, finaldistributions];
finalpdf = PDF[finalmixture, {x, y}];

(* Plot it. The ^0.125 is for better scaling and visualization only *)
DensityPlot[finalpdf^0.125, {x, 0, 12}, {y, 0, 12}, 
 Epilog -> {Red, Point[iris]}, PlotPoints -> 50]


It looks mostly correct by eye. I'm not sure if my transformations of the Cholesky covariances are correct. Somebody with more stats knowledge might want to verify this is the right way to un-standardize a covariance matrix.

  • $\begingroup$ I'm going to accept your answer as it answers the question. But I don't think this is the solution I am looking for, the R libraries Implement maximum-likelihood approaches, and I must confess to having no idea how LearnDistribution works! $\endgroup$
    – Q.P.
    Nov 13 '20 at 17:34
  • 1
    $\begingroup$ @Q.P. thanks. Please leave the question up though as this standardization stuff is useful for the community to know. The max likelihood approach is best done in Mathematica with FindDistributionParameters[..., ParameterEstimator -> "MaximumLikelihood"] (the "MaximumLikelihood" is the default anyway). See my other answer: mathematica.stackexchange.com/a/234210/72682 $\endgroup$
    – flinty
    Nov 13 '20 at 18:07
  • 1
    $\begingroup$ @Q.P. by the way, from looking at the ld[[1]] I can see it produces this key: "LossName" -> "MeanCrossEntropy". See here and here . Also this appears: "EMIterations" -> 1 so it's doing expectation maximization. $\endgroup$
    – flinty
    Nov 13 '20 at 18:13
  • $\begingroup$ That's a very good spot!! Can you see how to get a likelihood number out as a measure of how good the model numbers are? $\endgroup$
    – Q.P.
    Nov 13 '20 at 18:35

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.