# 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. –  Oleksandr R. 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. –  Oleksandr R. 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.) –  Mr.Wizard 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. –  Oleksandr R. Jul 11 '13 at 20:51
–  Sjoerd C. de Vries 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. –  Kai Tian Jul 12 '13 at 3:32