15
$\begingroup$

I have data that I want to fit with a constant plus a combination of Gaussians, with their actual number being one parameter of the fit (besides the average and the variance of every one of them).

How can I carry out and test such procedure in Mathematica?

If fitting discrete values is not possible, how can implement a series of fits for different number of Gaussians and a selection criteria the best fit?

A representative example of the data:

dat = Uncompress@Import["http://pastebin.com/raw.php?i=VHZ7XJAi"];
ListPlot[dat, PlotRange -> All, PlotTheme -> "Scientific"]

Mathematica graphics

Improve this question
$\endgroup$
14
  • 2
    $\begingroup$ Please describe your data. $\endgroup$
    – Sjoerd C. de Vries
    Commented Sep 7, 2015 at 14:57
  • 2
    $\begingroup$ If the number of Gaussians is a free parameter, wouldn't the search algorithm will most likely make it diverge? $\endgroup$
    – rhermans
    Commented Sep 7, 2015 at 15:04
  • 3
    $\begingroup$ @rhermans: I don't think the number of Gaussians would diverge in the strictest sense; but since each Gaussian has 2 free parameters, a sum of $(n-1)/2$ Gaussians plus a constant will be able to fit a generic set of $n$ data points perfectly. $\endgroup$
    – Michael Seifert
    Commented Sep 7, 2015 at 15:17
  • 5
    $\begingroup$ I don't think that there's any support for fitting to discrete parameters in Mathematica (though I could be wrong.) You might consider doing a series of fits to 1, 2, 3, ... Gaussians (which Mathematica can do) and then applying the Aikaike information criterion or a similar criterion to pick the best one. $\endgroup$
    – Michael Seifert
    Commented Sep 7, 2015 at 15:21
  • 4
    $\begingroup$ @rhermans: Not if you run out of data points to fit. For example, you could fit $n$ data points $(x_i, y_i)$ perfectly with $n$ Gaussians by letting them have infinitesimally small width, centering the $i$th Gaussian on $x_i$, and setting its height to be $y_i$. Adding another Gaussian on top of this wouldn't increase the goodness of the fit. Allowing the widths of the Gaussians to vary would presumably reduce the number required (though I'm not 100% about the estimate I gave above.) $\endgroup$
    – Michael Seifert
    Commented Sep 7, 2015 at 15:31

1 Answer 1

Reset to default
35
$\begingroup$

Solution

Following @MichaelSeifert's advice, here is a working solution.

First define a compact expression for a Gaussian

g[x_, xo_, σ_, a_] :=  a Exp[-((x - xo)^2/(2 σ^2))] /(σ Sqrt[2 π])

To get

g[x, xo, σ, Ao]

Mathematica graphics

Creates a list of variables for the $k^{\text{th}}$ term.

kvar[k_Integer] := 
 ToExpression@
  Map[StringJoin[#, ToString[k]] &, {"x", "σ", "a"}]

Mathematica graphics

Generate a model equation to fit $n$ Gaussians using Sequence to place the output of kvar as arguments for g

gmodel[n_Integer] := Sum[
  g[x, Sequence @@ kvar[i]]
  , {i, 1, n}
  ]

Generate a list of parameters for the $n$ Gaussian model.

gpars[n_Integer] := Flatten@Array[kvar, n]

When evaluated these functions look like this:

Mathematica graphics

Find fit with minimal Akaike Information Criterion ( AIC ) by calculating a series of fits with 1, 2, ...maxn Gaussian, and selecting the fit with the smallest "AIC" as defined in the NonlinearModelFit documentation.

fitg[data_, maxn_Integer] := MinimalBy[
   Table[
    {#, #["AIC"]}& @ NonlinearModelFit[data, gmodel[n], gpars[n], x]
    , {n, maxn}
    ], Last][[1, 1]]

Usage: fitg[data, maxn] where maxn is the maximum allowed number of Gaussians.

Possibly fitg[data, Length[data]]

Test

Example of data

dat = Uncompress@Import["http://pastebin.com/raw.php?i=VHZ7XJAi"]

Find fit with up to 6 Gaussians and plot

Show[
 ListPlot[dat, PlotStyle -> Red],
 Plot[Evaluate[Normal[fitg[dat, 6]]], {x, -1, 2}]
 ]

Mathematica graphics

fitg[dat, 6]["ParameterTable"]

Mathematica graphics

In this case, only three Gaussians were necessary.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ This may need more work if you need to get clever with initial guesses, but as you have not shared your data.... $\endgroup$
    – rhermans
    Commented Sep 7, 2015 at 16:16
  • 2
    $\begingroup$ Very nice. I'm not sure I could have done it quite so well myself. $\endgroup$
    – Michael Seifert
    Commented Sep 7, 2015 at 19:28
  • 1
    $\begingroup$ @rhermans thank you for this instructive answer :) $\endgroup$
    – ubpdqn
    Commented Sep 8, 2015 at 21:29
  • 1
    $\begingroup$ Please notice that there are things to improve upon here. See this question about why to use indexed variables instead of ToExpression . I should have done kvar[k_Integer] := Through[{xo, σ, a}[k]] instead. $\endgroup$
    – rhermans
    Commented Sep 9, 2015 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.