# How to make a multi-peak histogram fitting

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}?

• No offence, but I really think this is far too basic a question for this site. If you want to learn how to fit data, just consult the documentation for FindFit and/or NonlinearModelFit. Jul 11 '13 at 19:22
• This question appears to be too basic for this site in its present form. If there are any subtleties involved, these have not been made clear enough for the question to be answerable. Jul 11 '13 at 19:24
• @OleksandrR. "We" used to answer some very basic and/or essentially mathematical questions on this site without complaint. Do you feel it is necessary to close this rather than simply ignore it if it doesn't interest you? (I ask honestly.) Jul 11 '13 at 20:26
• @Mr.Wizard by voting to close I do not demand that it is closed. Others may, as they see fit, either agree with me and vote to close, or disagree and write an answer. The question is not only very basic but also not especially well posed, and I notice that the number of questions sharing these characteristics is rapidly increasing. I voted to close because I have no desire for the site to become dominated by simple questions that the askers can surely solve by themselves with minimal (but still non-zero) effort. Such questions are a waste of time for everyone else. Jul 11 '13 at 20:51
• Jul 11 '13 at 21:20

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]}]


• Got it. Thanks. Jul 12 '13 at 3:32

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]


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[

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

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