# How to perform a multi-peak fitting?

I am wondering how to implement the multi-peak detecting and fitting in Mathematica. Following is an example of fitting the data using three peaks (such that the data ~ peak1 + peak2 + peak3).

The peak model is given and fixed (all peaks are fitted by the same model), but its particular form (which will be input) can be Gaussian or Lorentzian or some other customized functions. The number of peaks is unknown and should be detected automatically, and the fitting model must also be built accordingly. Is there a Mathematica function that can simply do this? Or if anyone can give an idea of how to do the multi-peak fitting using Mathematica.

(I am aware of fitting functions like FindFit,NonlinearModelFit etc., so my question is more about how to build the model and estimate the initial parameters for input of the fitting functions.)

I am expecting something like this:

PeakFit[data_, pfun_, x_]:=...


where the data is a list of points like {{x1_,y1_}..}, x_ specifies the variable to be used, and the peak function pfun is a pure function whose first three parameters control the peak height, the peak width, and the central position, and the remaining (optional) parameters are for further control of the shape of the peak. For example a Gaussian model may be described as

pfun = Function[{x}, #1 Exp[-(x - #3)^2/(2 #2^2)]] &;


Given the data and the peak function, I wish PeakFit to return a FittedModel object containing the resulting model like pfun[A_,w_,xc_][x]+....

-
All peaks are fitted by a given peak model which will be input. For this question, let us first not consider the casing of varying peak functions. Thanks for the comments, the question is modified to make this point clear. –  Everett You Jun 3 '13 at 18:57
Studying physics, I've wondered the same thing. –  rcollyer Jun 3 '13 at 19:00
MixtureDistribution will likely be helpful –  chuy Jun 3 '13 at 19:07
related –  chris Jun 3 '13 at 20:29

It is possible to include the number of peaks (denoted $n$ below) in minimum searching.

First we create some test data:

peakfunc[A_, μ_, σ_, x_] = A^2 E^(-((x - μ)^2/(2 σ^2)));

dataconfig = {{.7, -12, 1}, {2.2, 0, 5}, {1, 9, 2}, {1, 15, 2}};
datafunc = peakfunc[##, x] & @@@ dataconfig;
data = Table[{x, Total[datafunc] + .1 RandomReal[{-1, 1}]}, {x, -20, 25, 0.1}];

Show@{
Plot[datafunc, {x, -20, 25},
PlotStyle -> ({Directive[Dashed, Thick,
ColorData["Rainbow"][#]]} & /@
Rescale[Range[Length[datafunc]]]), PlotRange -> All,
Frame -> True, Axes -> False],
Graphics[{PointSize[.003], Gray, Line@data}]}


Then we define the fit function for a fixed $n$ using Least Squares criterion:

Clear[model]
model[data_, n_] :=
Module[{dataconfig, modelfunc, objfunc, fitvar, fitres},
dataconfig = {A[#], μ[#], σ[#]} & /@ Range[n];
modelfunc = peakfunc[##, fitvar] & @@@ dataconfig // Total;
objfunc =
Total[(data[[All, 2]] - (modelfunc /. fitvar -> # &) /@
data[[All, 1]])^2];
FindMinimum[objfunc, Flatten@dataconfig]
]


And an auxiliary function to ensure $n\geq 1$:

Clear[modelvalue]
modelvalue[data_, n_] /; NumericQ[n] := If[n >= 1, model[data, n][[1]], 0]


Now we can find the $n$ which minimizes our goal:

fitres = ReleaseHold[
Hold[{Round[n], model[data, Round[n]]}] /.
FindMinimum[modelvalue[data, Round[n]], {n, 3},
Method -> "PrincipalAxis"][[2]]] // Quiet


Note:

For this example, the automatic result shown above is not that good:

resfunc =
peakfunc[A[#], μ[#], σ[#], x] & /@ Range[fitres[[1]]] /. fitres[[2, 2]]

Show@{
Plot[Evaluate[resfunc], {x, -20, 25},
PlotStyle -> ({Directive[Dashed, Thick,
ColorData["Rainbow"][#]]} & /@
Rescale[Range[Length[resfunc]]]), PlotRange -> All,
Frame -> True, Axes -> False],
Plot[Evaluate[Total@resfunc], {x, -20, 25},
PlotStyle -> Directive[Thick, Red], PlotRange -> All,
Frame -> True, Axes -> False],
Graphics[{PointSize[.003], Gray, Line@data}]}


To solve the problem, we can design a penalty function, so when increasing $n$ gains relatively little, we will prefer the smaller $n$.

Here I don't present the penalty function, but only show the phenomenon it based on. Please note after $n$ achieves $4$, which is the correct peak number, the modelvalue decreases much more slower.

{#, modelvalue[data, #]} & /@ Range[1, 7] // ListLogPlot[#, Joined -> True] & // Quiet


With[{n = 4},
resfunc = peakfunc[A[#], μ[#], σ[#], x] & /@ Range[n] /. model[data, n][[2]] ]

Show@{
Plot[Evaluate[resfunc], {x, -20, 25},
PlotStyle -> ({Directive[Dashed, Thick,
ColorData["Rainbow"][#]]} & /@
Rescale[Range[Length[resfunc]]]), PlotRange -> All,
Frame -> True, Axes -> False],
Plot[Evaluate[Total@resfunc], {x, -20, 25},
PlotStyle -> Directive[Thick, Red], PlotRange -> All,
Frame -> True, Axes -> False],
Graphics[{PointSize[.003], Gray, Line@data}]}


-
Thanks for the great answer. How to design the penalty function? Something proportional to $n$, or any other functions? –  Everett You Jun 6 '13 at 4:38
@EverettYou You're welcome. I think the penalty function design is straightforward from the phenomenon I shew. I know you're a student, I think this remaining problem would be a good exercise on scientific numerical computing. –  Silvia Jun 6 '13 at 5:33
One might consider designing a method built on top of Silvia's proposal: starting from an initial estimate of $n$, do the fitting and note the $\chi^2$ value. Increment $n$, fit again, and see if $\chi^2$ decreases. Keep doing this until you see an increase in $\chi^2$, and then keep the fit that had the least value of $\chi^2$. –  Ｊ. Ｍ. Jun 6 '13 at 14:55
@0x4A4D In the absence of numerical problems, shouldn't $\chi^2$ keep decreasing as $n$ increases? More parameters make for fits with smaller residuals? Perhaps, like an F test, $n$ should be chosen as the max value such that $\chi^2$ does not significantly decrease (whatever significantly means). –  KennyColnago Jun 7 '13 at 2:43
@Kenny, ah, that's a much better formulation. Indeed, if adding the $n+1$-th term does not change $\chi^2$ much, then one can comfortably settle for the $n$-term model. –  Ｊ. Ｍ. Jun 7 '13 at 2:58

My interpretation of your question is that you want to fit a linear combination of peaked functions with non-negative coefficients.

Beware: The minimum misfit solution with non-negative coefficients is a few isolated delta-functions. Therefore, allowing peaks widths is useless, whether for least square or least absolute error, because the minimum allowed width, most resembling a delta-function, will always be chosen.

You say your question is more about initial parameter estimates and detecting peaks...

Nonlinear methods sometimes require a guess at the number of peaks, and initial values for their positions and amplitudes. Convergence may be a problem. However, a linear inversion is possible if the horizontal coordinate is specified as a vector of values. Then the algorithm searches only for the peak amplitudes at every one of these values, a linear fit. Most amplitudes will be zero (again, because the minimum misfit solution is a few, isolated delta functions). In addition, this linear method is not biased by a specification of the number of peaks.

I've used the Mathematica implementation of the non-negative least-squares algorithm NNLS of Lawson and Hanson for decades. It was written by Michael Woodhams, and is on the MathGroup Archive 2003.

-
+1. Please provide links for the algorithm NNLS of Lawson and Hanson and its implementation by Michael Woodhams. –  Alexey Popkov Jun 3 '13 at 20:39
C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974. netlib.org/lawson-hanson/all –  KennyColnago Jun 3 '13 at 21:11

The question is not so innocent as is appears. Without a penalty on the number of peaks the "best" model is overfitting the data. The answer by Silvia demonstrates this already. And, think about it, you got what you wanted: adding more peaks will fit the data better. Always!

One may revert to adding an ad-hoc penalty function on the number of peaks. But this is often unsatisfactory; nagging doubts may remain after seeing the results. Therefore, I would like to point your attention into the direction of Bayesian model selection. Model fitting and selection are two parts of the same theory - no ad-hockeries.

The "bad" news is you have to un-learn statistics and to learn Bayesian probability theory. And, yes, learn how to transform your "state of knowledge" about the problem into prior probabilities. However, this is easier than you might think.

The "good" news is that it works. E.g. I have seen satellite spectra fitted with hundreds of peaks, while simultaneously estimating the calibration parameters of the instrument being far out of reach. A hopeless task without a systematic guidance by probability theory, in my view. However, don't underestimate the computational burden. Such models may require hours-days-weeks of CPU time. Don't be put off by this, in my experience this is worth it. The Bayesian approach delivers in real scientific life, but not for the faint at heart.

In brief, how does this work. The Likelihood p(D|M) of the data D given a model M with, say, 4 peaks is p(D|M=4). (The "given" is denoted by "|".) Maximizing the Logarithm of this Likelihood, by adjusting the positions and widths of the peaks, is exactly the same as minimizing the least square error! (See book of Bishop, below.) But the Maximum Likelihood values of p(D|M=4) < p(D|M=5) < p(D|M=6) < ... , etc. Until the number of peaks equal the number of data and the least square error becomes zero.

In Bayesian model selection, the probability p(M=4|D) of a model M having 4 peaks given the data D is a viable concept. (Note the reversal of M and D about the |.) The value of the ratio of e.g. p(M=5|D)/p(M=4|D) gives a measure whether model M=5 is better than M=4. Bayes theorem yields p(M=5|D)/p(M=4|D) = p(D|M=5)/p(D|M=4) * "Ockham factor", where we recognize the above ratio of Likelihoods, which is >1 in this example.

The "Ockham factor" includes the penalties, which typically contain a ratio Exp[4]/Exp[5] < 1 from the number of M peaks in this example. The balancing between the Likelihood ratio p(D|M=5)/p(D|M=4) and the "Ockham factor" determines the most probable model. If p(M=5|D)/p(M=4|D) < 1 then the model with fewer peaks M=4 is a better model than M=5.

Anyone interested may have a look at two excellent books. 1) Data analysis, a Bayesian tutorial, by D.S. Sivia with J. Skilling (http://amzn.to/15DnwV3), and 2) Pattern Recognition and Machine Learning by C.M. Bishop (http://amzn.to/13n67ji).

-
+1 Bayesian method is a good choice if one wants to solve the problem seriously. –  Silvia Jun 6 '13 at 16:40

Here is a simple example of fitting a series of three gaussians. I hope this is helpful.

Make some data with some added noise:

f = (7/10) PDF[NormalDistribution[-12, 2], x] + PDF[NormalDistribution[0, 1], x] + (5/10) PDF[NormalDistribution[9, 3], x]
data = Table[{x, f + RandomReal[0.01]}, {x, -20, 20, 0.1}];


Set up the model, variables and a few constraints (had to tweak a few initial guesses to get it to converge):

model = (amp1 E^(-((-mu1 + x)^2/(2 sigma1^2))))/(
Sqrt[2 \[Pi]] sigma1) + (amp2 E^(-((-mu2 + x)^2/(2 sigma2^2))))/(
Sqrt[2 \[Pi]] sigma2) + (amp3 E^(-((-mu3 + x)^2/(2 sigma3^2))))/(
Sqrt[2 \[Pi]] sigma3);
vars = {{amp1, 1}, {mu1, -10}, sigma1, amp2, mu2, sigma2,
amp3, {mu3, 10}, sigma3};
cons = And @@ Thread[{amp1, sigma1, amp2, sigma2, amp3, sigma3} > 0];


Do the fit:

fit = NonlinearModelFit[data, {model, cons}, vars, x, MaxIterations -> 100]


Plot the result:

Show[Plot[fit[x], {x, -20, 20}, PlotRange -> All, PlotStyle -> Thick],
ListPlot[data, PlotStyle -> Opacity[0.5]]]


-
Now the important question: how does one automate that process? –  rcollyer Jun 3 '13 at 19:06
Thanks for the nice answer, this should be the final part of the fitting. My question is more about how to detect the peak automatically and to estimate the initial parameters accordingly. –  Everett You Jun 3 '13 at 19:34

I would despair that it is possible to build a method that can find an arbitrary number of peaks in the same way that your eye does it. Here's one approach that can work in certain situations where the width of the peaks is roughly known. The idea is simple: find the highest point on the data and use that to initialize the location of the first peak. Once the first best fit peak is found, subtract it from the data. Then repeat. With some luck the collection of peaks can be recovered. Begin by creating some data (following the construction by sOrce).

 f = PDF[NormalDistribution[100, 10], x] + (2/3) PDF[NormalDistribution[250, 20], x];
data = Table[f + RandomReal[0.01], {x, 1, 500, 1}];


Use the function

indMax[q_] := Ordering[q][[Length[q]]];


to locate the position (index) of the maximum point in the data.

len = Length[data];
indD = indMax[data]; max = Max[data];
nlmD = NonlinearModelFit[data, b Exp[-0.002 (x - indD)^2], {{b, max/2}}, x];
curveD = Table[nlmD[x], {x, 1, len, 1}];
data2 = Clip[data - curveD, {0, 1}];
indD2 = indMax[data2]; max2 = Max[data2];
nlmD2 = NonlinearModelFit[data2, b2 Exp[-0.002 (x - indD2)^2], {{b2, max2/2}}, x];
curveD2 = Table[nlmD2[x], {x, 1, len, 1}];
data3 = Clip[data2 - curveD2, {0, 1}];


Here we have done it twice, finding the first two peaks. To see how things are going:

Show[ListLinePlot[{data, curveD, curveD2},
PlotRange -> {{1, len}, All}, PlotStyle -> {Blue, Red, Red}],
Graphics[{PointSize[Large], Orange, Point[{indD, max}], Point[{indD2, max2}]}]]


This shows the detected maxima (the orange dots) and the fitted normal exponentials along with the data. In this data, (since there really only are two peaks) applying the process again will just give more (and smaller) peaks which can be weeded out by some kind of threshold, perhaps chosen from a knowledge of the noise floor.

To recap: this iteration does not need to know how many peaks there are nor where they are located. It does assume a fixed width for the peaks, and it assumes that enough is known about the nature of the problem to be able to stop the iteration.

-