# Defining a Gaussian mixture distribution for $N$ number of Gaussian distributions

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$$.

• 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. Nov 6, 2020 at 18:31
• 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. Nov 6, 2020 at 18:35
• Also, @flinty why not provide what you supplied in pastebin as an answer. It might be helpful for others. Nov 6, 2020 at 18:37
• 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. Nov 6, 2020 at 18:45

(* 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,
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,