0
$\begingroup$

I am experimenting a bit with mixed-Gaussian distributions. If I want to define a mixture of two or three distributions, I do this with:

GaussianMixture2[w1_, w2_, u1_, u2_, s1_, s2_] := Evaluate[MixtureDistribution[{w1, w2}, {NormalDistribution[u1, s1], NormalDistribution[u2, s2]}]]
GaussianMixture3[w1_, w2_, w3_, u1_, u2_, u3_, s1_, s2_, s3_] := Evaluate[MixtureDistribution[{w1, w2, w3}, {NormalDistribution[u1, s1], NormalDistribution[u2, s2], NormalDistribution[u3, s3]}]]

However if I want to do this for $N = 10$ distributions, it becomes to get clumsy and silly. Does anyone know how I can define a function similar to those above for $N$ number of distributions.

I want to be able to use the function in FindDistributionParameters, which I can with the above functions, and get a value for each individual $u_{i}$ and $s_{i}$. In addition, It would be nice to ensure that the sum of all weights, $w_{i}$, must be equal to $1$.

$\endgroup$
4
  • 1
    $\begingroup$ You can easily put your parameters in lists like this pastebin.com/wEFvEgyS . If you want undefined parameters, replace the random parameters with Array[w,n],Array[m,n],Array[s,n]. I'm of the view it's not realistic to expect FindDistributionParameters to work with 30 parameters like this because it's far too slow. $\endgroup$
    – flinty
    Nov 6, 2020 at 18:31
  • $\begingroup$ Thanks for this! It seems to work fine for me, it takes maybe 40 mins of computation when I set precision and accuracy goals to say 1 or 2. Using NMaximise helped as well. $\endgroup$
    – user27119
    Nov 6, 2020 at 18:35
  • $\begingroup$ Also, @flinty why not provide what you supplied in pastebin as an answer. It might be helpful for others. $\endgroup$
    – user27119
    Nov 6, 2020 at 18:37
  • 1
    $\begingroup$ Done, also the weights summing to 1 doesn't really matter. You can feed them in un-normalized I think and renormalize the resulting weights without much trouble. $\endgroup$
    – flinty
    Nov 6, 2020 at 18:45

1 Answer 1

2
$\begingroup$
(* generate some data *)
SeedRandom[1];
n = 10;
means = RandomReal[{-5, 5}, n];
stddevs = RandomReal[{0.1, 2}, n];
weights = Normalize[RandomReal[1, n], Total];
mixture = MixtureDistribution[weights, 
   MapThread[NormalDistribution[#1, #2] &, {means, stddevs}]];
data = RandomVariate[mixture, 10000];
Histogram[data]

(* create the parameters *)
paraMeans = Array[m, n];
paraStddevs = Array[s, n];
paraWeights = Array[w, n];
(* random initial values *)
weightsInit = Transpose[{paraWeights, Normalize[RandomReal[1, n], Norm]}];
meansInit = Transpose[{paraMeans, ConstantArray[Mean[data], n]}];
stddevsInit = Transpose[{paraStddevs, RandomReal[1, n]}];

(* create the mixture and fit *)
mixturemodel = MixtureDistribution[paraWeights, 
   MapThread[NormalDistribution[#1, #2] &, {paraMeans, paraStddevs}]];
guessedParams = 
 FindDistributionParameters[data, mixturemodel, 
  Flatten[{weightsInit, meansInit, stddevsInit}, 1], AccuracyGoal -> 2]

(* try it out and compare *)
Histogram[RandomVariate[mixturemodel /. guessedParams, 10000]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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