I got a multipeak histogram from my data and want to fit it using a Gaussian Distribution function.
(f[x_, amplitude_, centroid_, sigma_] :=
amplitude Exp[-((x - centroid)^2/sigma^2)])
However, there are two peaks. How do I fit the histogram and get the two sets of values {amplitude1, centroid1,sigma1} and {amplitude2, centroid2,sigma2}?
Since nobody who supposedly thought this was a good question seems actually to have wanted to write an answer, here is one from me.
Let us define a bimodal distribution:
dist = MixtureDistribution[
{0.6, 0.4},
{NormalDistribution[-0.8, 1.3], NormalDistribution[2.7, 0.4]}
];
pdfplot = Plot[PDF[dist, x], {x, -5, 5}]
To simulate your data we draw some random variates from this distribution. HistogramList is used to convert these into a set of bin and count specifications. (For your real data, you may find it useful to try different binning methods. The various possibilities are described in the documentation.)
{bins, counts} = HistogramList@RandomVariate[dist, 10000];
For fitting we need a list of bin centres and probabilities, rather than boundaries and counts:
data = Transpose[{
Mean /@ Partition[bins, 2, 1],
counts Length[bins]/(Total@Differences[bins] Total[counts])
}];
Show[{pdfplot, ListLinePlot[data, PlotStyle -> Red]}]
Now we get to the substance of your question, i.e. actually doing the fit:
f[x_, a_, μ_, σ_] := a Exp[-(x - μ)^2/(2 σ^2)]/(Sqrt[2 Pi] σ);
model = f[x, a1, μ1, σ1] + f[x, a2, μ2, σ2];
result = FindFit[data, model,
(* reasonable initial values are important *)
{{a1, 0.5}, {μ1, -1}, {σ1, 1}, {a2, 0.5}, {μ2, 3}, {σ2, 0.5}},
x
]
(* -> {a1 -> 0.618350, μ1 -> -0.798249, σ1 -> 1.314220,
a2 -> 0.404240, μ2 -> 2.693320, σ2 -> 0.400108} *)
It matches the original:
Show[{pdfplot, Plot[model /. result, {x, -5, 5}, PlotStyle -> Red]}]
When one knows that the underlying distribution is a mixture of a known number of normal distributions, then what @SvyatoslavKorneev gives works great. However, with any reasonable sample size, there is no reason to have to rely on the assumption that the underlying bumps are the result of a number of normal distributions with the number of distributions matching what you think you see in a histogram.
A nonparametric kernel density estimate is the way to go. Mathematica has the SmoothKernelDistribution function to do so:
(* Define mixture distribution *)
mixtureDist =
MixtureDistribution[{0.4, 0.3, 0.3}, {NormalDistribution[4, 1],
NormalDistribution[10, 2], NormalDistribution[20, 3]}];
(* Sample size *)
n = 10^5;
(* Generate a random sample *)
tbl = RandomVariate[mixtureDist, n];
(* Get adaptive kernel estimate of density function *)
skd = SmoothKernelDistribution[
tbl, {"Adaptive", Automatic, Automatic}];
(* Plot true density and estimated density *)
Plot[{PDF[mixtureDist, x], PDF[skd, x]}, {x, Min[tbl], Max[tbl]},
PlotStyle -> {{LightGray, Thickness[0.02]}, Red}]
The thick gray line is the true density and the red line is the estimated density.
Yes, a kernel density estimate does not result in a nice compact formula but it also relies on fewer potentially unjustified assumptions (and likely fits better). (And, yes, I'm not answering the actual question asked.)
Note that "larger" sample sizes are needed if one can't assume a known number of normal distributions as that assumption essentially provides a whole lot of information. Here's an example with the same underlying mixture distribution but with just 100 samples:
Update
I found a comment by @AndyRoss that allows for the finding of peaks which is something else that a kernel density estimate can provide (besides means, variances, semi-interquartile ranges, etc.). That comment suggests to replace SmoothKernelDistribution with KernelMixtureDistribution. The output object allows one to estimate both the first and second derivative of the density function.
Below is the updated code that finds peaks (if any) that the data seems to support:
(* Get adaptive kernel estimate of density function *)
kmd = KernelMixtureDistribution[
tbl, {"Adaptive", Automatic, Automatic}];
(* Plot true density and estimated density *)
Plot[{PDF[mixtureDist, x], PDF[kmd, x]}, {x, Min[tbl], Max[tbl]},
PlotStyle -> {{LightGray, Thickness[0.02]}, Red}, PlotRange -> All]
(* Find peaks *)
(* First and second derivatives *)
d1[x_] := D[PDF[kmd, y], y] /. y -> x
d2[x_] := D[d1[y], y] /. y -> x
(* Find where first derivative is zero *)
xmin = Min[tbl];
xmax = Max[tbl];
sol = NSolve[{d1[x] == 0, xmin < x < xmax}, x]
(* Horizontal and vertical axes values corresponding to peaks *)
xPeaks = Pick[x /. sol, Negative[d2[x /. sol]]]
yPeaks = PDF[kmd, xPeaks]
You can write Expectation-Maximization algorithm for Gaussian Mixture Models. It is very easy. Let's consider the case with 3-modal Gaussian distribution, the extension to the arbitrary number is straight forward.
Define the list with 3-modal Gaussian random numbers
n = 10^5;
tbl = RandomVariate[
MixtureDistribution[{0.4, 0.3, 0.3}, {NormalDistribution[4, 1],
NormalDistribution[10, 2], NormalDistribution[20, 3]}], n];
Plot the PDF of the list
SmoothHistogram[tbl, PlotRange -> All]
Define Gaussian distribution
p[x_, μ_, σ_] :=
1/(σ Sqrt[2 π]) Exp[-((x - μ)^2/(2 σ^2))];
Initiate parameters for the algorithm
α1 = 0.3;
α2 = 0.3;
α3 = 0.4;
μ1 = 2;
μ2 = 8;
μ3 = 14;
σ1 = 1;
σ2 = 1;
σ3 = 1;
Loop until convergence, but I use just a finite number of iterations
Do[
denom = p[tbl, μ1, σ1]*α1 + p[tbl, μ2, σ2]*α2 + p[tbl, μ3, σ3]*α3;
ω1 = p[tbl, μ1, σ1]*α1/denom;
ω2 = p[tbl, μ2, σ2]*α2/denom;
ω3 = p[tbl, μ3, σ3]*α3/denom;
n1 = Total[ω1];
n2 = Total[ω2];
n3 = Total[ω3];
α1 = n1/n;
α2 = n2/n;
α3 = n3/n;
μ1 = 1/n1 Total[ω1*tbl];
μ2 = 1/n2 Total[ω2*tbl];
μ3 = 1/n3 Total[ω3*tbl];
σ1 = Sqrt[1/n1 Total[ω1*(tbl - μ1)^2]];
σ2 = Sqrt[1/n2 Total[ω2*(tbl - μ2)^2]];
σ3 = Sqrt[1/n3 Total[ω3*(tbl - μ3)^2]];
Print[μ1, " ", μ2, " ", μ3], {j, 1, 30}]
Plot the result
Show[SmoothHistogram[tbl],
Plot[α1 p[x, μ1, σ1] + α2 p[
x, μ2, σ2] + α3 p[
x, μ3, σ3], {x, -10, 30}, PlotStyle -> Red,
PlotRange -> All], PlotRange -> All]
FindFitand/orNonlinearModelFit. $\endgroup$