How to perform a multi-peak fitting?

Question

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]+....

Answer 1

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

Answer 2

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

Answer 3

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.

Answer 4

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

Answer 5

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.

Answer 6

If you know your peak model a priori (say that you have the area normalized version in a vector psf) then you could also put this psf in a convolution kernel H formed by shifted copies of the psf along the diagonal so that your signal y = Hx and then use best subset regression (L0 pseudonorm penalized regression) to estimate x, taking advantage that the original signal x should have been a sparse spike train and that most coefficients should be estimated as zero. This approach was described in the paper "Deconvolution of pulse trains with the L0 penalty" by Johan de Rooi & Paul Eilers. I've got R code implementing this method here (sorry, not in Mathematica, but should be easy enough to port). In an extended method, the paper also shows how the peak model (if it's identical for all peaks) can be estimated from the data. See also here for a nice recent paper with a description of some possible approaches. There is also a new R package, l0ara, that can fit L0 penalized linear models. Though you would have to edit the code a bit to take into account the nonnegativity constraints on your coefficients that you would like in this case. If your noise is Poisson this can be taken into account by using observation weights 1/(y+1) (i.e. one over variance weights), and you can multiply your covariate matrix and your dependent variable y by sqrt(weights) (even better would be to use an L0 penalized GLM with identity link function, but l0ara only has Poisson with log link).