0
$\begingroup$

I have been trying to put toghether a visual method to help fit experimental data with peaks. I have stumbled upon this answer which is very similar to where I got. But it is not sufficiently general for what I need. So I reworked it into a simplified version that has the properties I need but which does not work.

The idea is to be able to have arbitrary number of functions and parameters put together in a model which is defined ahead of time. So first we define the same data and model as in the answer referenced above but in a more general way with indexed variables.

ClearAll["Global`*"]
fdata = Table[{x, 
    8 x^3 - 7 x^2 - 10 x + 1 + RandomReal[{-5, 5}]}, {x, -2, 2, 
    0.1}];
n = 3;
model = Sum[a[i] x^i, {i, 0, n}]

Note that model is identical to the one used in the referenced answer for simplification. At the second part of the post there is one of the full models that I am working with to illustrate what I really need.

And now the modified DynamicModule section from the reference above:

DynamicModule[
 {sol},
 With[
  {
   manlocalmodel = model /. (a[#] -> av[#] & /@ Range[0, n]),
   fitlocalmodel = model /. (a[#] -> af[#] & /@ Range[0, n]),
   data = (av[#] & /@ Range[0, n])
   },
  Manipulate[
   If[
    computeFlag == True,
    sol =
     FindFit[
      fdata,
      fitlocalmodel,
      {af[#], data[[# + 1]]} & /@ Range[0, n],
      x
      ];
    Print[sol];
   (av[#]=(af[#]/. sol))&/@Range[0,n];
    computeFlag = False;
    ];
   Column[
    {
     Dynamic[Button["Compute", computeFlag = True]],
     Show[
      Plot[manlocalmodel, {x, -2, 2}, PlotStyle -> Black],
      ListPlot[fdata, PlotStyle -> Red],
      ImageSize -> 300, PlotRange -> {{-2.05, 2.05}, All}
      ]
     }
    ],(*Manipulate variables*)
   {
    {computeFlag, False}, ControlType -> None},
   Evaluate[
    Sequence @@ ({{av[#], 1}, -10, 10, Appearance -> "Open"} & /@ 
       Range[0, n])]
   ] (*end of Manipulate*)
  ] (* end of With *)
 ] (*end of DynamicModule*)

The problem with this code is that it does not feed back the fitted values into the Manipulate and the responsible line is:

av[#] = (af[#] /. sol)) & /@ Range[0, n];

Now, if I replace this line with the explicit assignments, it works, that is:

av[0] = af[0] /. sol;
av[1] = af[1] /. sol;
av[2] = af[2] /. sol;
av[3] = af[3] /. sol;

I have been doing some alchemy with this but to no avail. I simply do not understand the scoping and Set and indexed variables and other basics sufficiently. So the question is, how to enable the feedback to the Manipulate with general number of indexed variables. Also, if someone want to play, below is the full code with data and maybe there are improvements that could help the code.


The full problem

The full problem code is as this and requires to get some data from the wolfram cloud. The first part of the code explicitly provides some data and at the end, the initial guesses so that fit can be found without needing to deal with setting up the sliders (this is to save time now). Then there is the DynamicModule part. In the version here, it works due to the explicit assignment of the fit solution to the indexed variables but this requires manual input and is not what I need.

Additionally, I am not even sure, if I am doing this the best way possible. Any suggestions to make the code below better are welcome (some interesting questions that pop up to me, is the With construct necessary within the DynamicModule?, should model be provided as a function of t_,μ_,a_,σ_ possibly with SetDelayed?, can the creation of controls be simplified? etc)

Note: I consider 2 types of peaks, those whose position and height I have determined (I use FindPeaks et al.) and those that are there (due to too broad shoulders and some other effects visible in the data) but whose position cannot be determined. The number of determined peaks is given by nd and their positions and amplitudes are saved in μd and ad respectibely. The number of undetermined peaks are given by nu and thus the number of each set of parameters (amplitude, position and variance) is nd + nu. Also, the fit is not ideal because the peaks are really not Gaussians but are asymetric, this is irrelevant at this point but the model should be able to handle any combination of custom functions (as long as FindFit or NonlinearModelFit can deal with the model).

ClearAll["Global`*"]
cdata = CloudGet[
   "https://www.wolframcloud.com/obj/8df13783-4827-4fd1-be7b-\
65875d9b02ee"];
fdata = cdata[[3097 ;; 52097]];
μd = {5422.6501499999995`, 6450.7998`, 6711.8498500000005`, 
   6824.2998`, 8795.4004`, 9502.0996`};
ad = {0.00267136`, 0.0210852`, 0.0228746`, 0.0311949`, 0.00351935`, 
   0.00333709`};
{t0, t1} = {4690.5`, 11110.5`};
dmmax = 0.0311949`;
plotbcg = 
  ListLinePlot[fdata, PlotRange -> {{t0, t1}, {0, dmmax + 0.1*dmmax}}];
Format[σ[i_]] := Subscript[σ, i]
Format[a[i_]] := Subscript[a, i]
Format[μ[i_]] := Subscript[μ, i]
nd = Length@μd;
nu = 2;
ad2 = ConstantArray[0, nd + nu];
ad2[[;; nd]] = ad;
μd2 = ConstantArray[{t0, t0, t1}, nd + nu];
μΔ = 200;
μd2[[;; nd, 1]] = μd;
μd2[[;; nd, 2]] = μd - μΔ;
μd2[[;; nd, 3]] = μd + μΔ;
model3 = 
 Sum[a[i]*Exp[-((t - μ[i])^2/(2 σ[i]^2))], {i, 1, nd}] + 
  Sum[a[i]*Exp[-((t - μ[i])^2/(2 σ[i]^2))], {i, nd + 1, 
    nd + nu}]
{ag, μg, σg} = 
  Transpose@{{0.00115`, 5422.6501499999995`, 
     50}, {0.011300000000000001`, 6450.7998`, 120.`}, {0.01535`, 
     6711.8498500000005`, 91.`}, {0.01145`, 6824.2998`, 
     11.`}, {0.00351935`, 8795.4004`, 176.`}, {0.00333709`, 
     9502.0996`, 104.`}, {0.01065`, 6270.`, 423.`}, {0.01105`, 6800, 
     24.`}};


DynamicModule[
 {sol},
 With[
  {
   localmodel = 
    model3 /. ((({a[#] -> 
             av[#], μ[#] -> μv[#], σ[#] -> \
σv[#]}) & /@ Range[nd + nu]) // Flatten),
   fitmodel = 
    model3 /. ((({a[#] -> 
             af[#], μ[#] -> μf[#], σ[#] -> \
σf[#]}) & /@ Range[nd + nu]) // Flatten),
   data = {av[#], μv[#], σv[#]} & /@ Range[nd + nu]
   },
  Manipulate[
   If[
    computeFlag == True,
    sol = FindFit[
      fdata,
      fitmodel,
      Evaluate@
       Flatten[{{af[#], data[[#, 1]]}, {μf[#], 
            data[[#, 2]]}, {σf[#], data[[#, 3]]}} & /@ 
         Range[nu + nd], 1],
      t];
    Print[sol];
    av[1] = af[1] /. sol; av[2] = af[2] /. sol; av[3] = af[3] /. sol; 
    av[4] = af[4] /. sol; av[5] = af[5] /. sol; av[6] = af[6] /. sol; 
    av[7] = af[7] /. sol; av[8] = af[8] /. sol;
    μv[1] = μf[1] /. sol; μv[2] = μf[2] /. 
      sol; μv[3] = μf[3] /. sol; μv[4] = μf[4] /. 
      sol; μv[5] = μf[5] /. sol; μv[6] = μf[6] /. 
      sol; μv[7] = μf[7] /. sol; μv[8] = μf[8] /. 
      sol;
    σv[1] = σf[1] /. sol; σv[
      2] = σf[2] /. sol; σv[3] = σf[3] /. 
      sol; σv[4] = σf[4] /. sol; σv[
      5] = σf[5] /. sol; σv[6] = σf[6] /. 
      sol; σv[7] = σf[7] /. sol; σv[
      8] = σf[8] /. sol;
    computeFlag = False;
    ];
   Show[
    Quiet@Plot[
      localmodel, {t, t0, t1},
      PlotStyle -> {Red},
      PlotRange -> {{t0, t1}, {0, dmmax + 0.1*dmmax}},
      PlotPoints -> 50,
      MaxRecursion -> 5
      ],
    plotbcg
    ],
   Evaluate[
    Sequence @@ 
     Flatten[{{{av[#], ag[[#]], Subscript[a, #]}, 0, dmmax, 
          Appearance -> "Labeled"}, {{μv[#], μg[[#]], 
           Subscript[μ, #]}, μd2[[#, 2]], μd2[[#, 3]], 
          Appearance -> "Labeled"}, {{σv[#], σg[[#]], 
           Subscript[σ, #]}, 0.01, (t1 - t0)/10, 
          Appearance -> "Labeled"}} & /@ Range[nd + nu], 1]
    ],
   Dynamic[Button["Compute", computeFlag = True]], 
   Dynamic[Button["Print", Print[data]]]
   ]
  ]
 ]
$\endgroup$
8
  • $\begingroup$ I am a big supporter of getting good starting values: "Good starting values are your best friends." But I don't see a description of why better starting values are needed here. (1) Do naïve starting values end up with clearly wrong results? (2) Do better starting values result in quicker execution times for NonlinearModelFit? (3) Is the problem that you have so much data (~ 50,000 data points) which are close together such that maybe sampling only every other 100-th data point would get either the answer or the better starting values? In short, why the need for better starting values? $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 16:42
  • $\begingroup$ @JimB The reason is that with blind values, FindFit does not converge. Additionally, there is the issue with peaks that are not clearly determined. One visually knows where they belong but Mathematica might fall into a different minimum $\endgroup$
    – atapaka
    Commented Jul 21, 2023 at 19:33
  • $\begingroup$ What counts is whether Mathematica ends up with the smallest root mean square error (or AICc) or your visual approach. "Visual" fits don't usually follow minimizing the sum of squares. So are you suggesting that one not use a standard metric for fitting and that you'll know it when you see it? (I know that sounds sarcastic but I'm really trying to understand what you think the issue is and if there's a basis for that.) $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 19:40
  • $\begingroup$ @JimB I am not sure how to answer that. I kept encountering FindFit does not converge with some "dumb" choice of initial guesses. Also, if I remember correctly minimization techniques such as least squares sometimes suffer from getting stuck in a local minimum, not sure if this is applicable for this type of data but it can happen. Choosing proper initial values should prevent this. $\endgroup$
    – atapaka
    Commented Jul 21, 2023 at 19:45
  • $\begingroup$ No disagreement with the usual need for good starting values. (My answer to mathematica.stackexchange.com/questions/165054/… shows why some starting values fail.) But if you have 50,000 data points, that also can cause problems. $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 20:00

1 Answer 1

1
$\begingroup$

It's not clear why better starting values (over naïve starting values) are needed for your "full problem". I'm not saying they are not needed, it just doesn't seem necessary for the "full problem" you give.

Because there are so many data points (around 50,000) and that there seems a very smooth transition between sample points, just taking 1/50th of the data points is more than adequate in this case.

cdata = CloudGet["https://www.wolframcloud.com/obj/8df13783-4827-4fd1-be7b-65875d9b02ee"];
data = cdata[[3097 ;; 52097]];

(* Choose every 50-th data point *)
data50 = data[[Range[1, Length[data], 50]]];

(* Number of Gaussian kernels to use *)
n = 10;

(* Create constraints: all "standard deviations" are positive and 
   the means are in consecutive order *)
constraints = And @@ Join[Table[c[i] > 0, {i, n}], {Less @@ Table[b[i], {i, n}]}]

(* Initial values:  spread out means evenly through the data,
   put all of the heights at the mean of the response, and
   make the "standard deviation" one-fourth of the distance between means *)
{xmin, xmax} = MinMax[data[[All, 1]]]
inits = Join[Table[{a[i], Mean[data[[All, 2]]]}, {i, n}],
  Table[{b[i], xmin + (xmax - xmin) i/(n + 1)}, {i, n}],
  Table[{c[i], (xmax - xmin)/(4 (n + 1))}, {i, n}]]

(* Define basis function *)
kernel[a_, b_, c_, x_] := a Exp[-(x - b)^2/(2 c^2)]

(* Perform fit *)
nlm = NonlinearModelFit[data50,
   {Sum[kernel[a[i], b[i], c[i], x], {i, 1, n}], constraints}, inits, x];

(* Plot resulting fit and all data *)
Show[ListPlot[data, PlotRange -> All, PlotStyle -> PointSize[0.01]],
 Plot[nlm[x], {x, xmin, xmax}, PlotStyle -> Red, PlotRange -> All]]

Data and fit

The fit "looks" OK but I'd argue only if the objective is to parsimoniously approximate the data as opposed to being able to make predictions for future data with some reliable measure of precision.

Why? The residuals are highly correlated but the model fit assumes uncorrelated errors. In that case the error predictions are unreliable. And there is no replication at any predictor value so it's impossible to estimate and separate the lack-of-fit error from the measurement error.

ListPlot[Transpose[{data50[[All, 1]], nlm["FitResiduals"]}], PlotRange -> All,
  Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"x", "Residual"}] // Quiet

Plot of residuals vs predictor

Update:

As we've both mentioned in comments, maybe the objective function of minimizing the sum of squares over the range of the predictors misses small but "known" peaks. Maybe putting more weight around locations where a "visible peak" occurs might be a way to go. I just don't know the subject matter to make the suggestion more concrete than that.

However, whatever fit is obtained, I would think the plotting the contributions of the individual peaks is a useful diagnostic. Here is such a plot for $n=10$ (which has a lower AICc value than for $n=8$). The blue dots represent the data, the green line represents the overall fit, and the red curves are the individual kernels.

Data, fit, and individual contributions

We see that the "bump/peak" near 5,400 is not included even with $n=10$. If some other set of 8 or 10 kernels is used that includes a peak at 5,400, then I would think that the overall mean square error would be higher because of not fitting as well elsewhere.

In short, I still don't think that the main issue is about starting values in this case but maybe an overall minimization of square of the errors just can't find all of the peaks and maybe an ad hoc weighting of points near observed peaks will produce more realistic results.

$\endgroup$
7
  • $\begingroup$ I will surely check your approach, especially your choice of initial values. However, one big question, why are you fitting 10 kernels? One of the reasons of the visual method suggested in the question is to avoid way too many peaks. In the posted data there are very clear 6 peaks. Additionally, there is a left shoulder on the highest peak suggesting another peak that is hidden. There needs to be at least one more to get the width right between t0 and approximately 7000. The ultimate goal is to match the model after integration (data is a DTG curve). $\endgroup$
    – atapaka
    Commented Jul 21, 2023 at 19:38
  • $\begingroup$ Also note that when I try to add 9th peak to my approach, FindFit tends to complain a lot about lack of convergence unless the peak is very specifically positioned initially. $\endgroup$
    – atapaka
    Commented Jul 21, 2023 at 19:39
  • $\begingroup$ I arbitrarily set the number of kernels to 10. Standard procedure would be to use AICc to choose a reasonable number of kernels. However, because you also have serial correlation that is not accounted for, that approach probably won't be satisfactory. $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 19:42
  • $\begingroup$ I avoid FindFit because it does not have any capability of producing standard errors or any diagnostic statistics as the "fit" always needs to be examined for meeting the assumptions of the model and error structure. In practice one typically just wants to know if there are major deviations from those assumptions. $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 19:45
  • 1
    $\begingroup$ If the underlying process that generates the data is something like researchgate.net/figure/…, then standard regression is not necessarily the way to go because identification of peaks is more important the fitting the rising or falling limbs or the low values. Or (I know with hydrocarbon chromatogarms for termites) it's the area under the peak associated with a certain hydrocarbon that's critical to estimate (as opposed to just the height of the peak). Life is complicated. $\endgroup$
    – JimB
    Commented Jul 21, 2023 at 20:18

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.