5
$\begingroup$

Short Version: I have a parameterized function (the model) that returns a list of values. I have data that I want to fit this model to, optimizing those parameters, where the list position in the data corresponds to the list position in the model result that is returned by the function. How do I do this using NonlinearModelFit so that I can use the various statistical analyses present in the resulting FittedModel?

Full Version: I have experimental data for a uniform, discrete set of times t. I want to fit this data to a nonlinear model including convolution with a known instrument function. As I determined in a prior question, Convolve fails (Wolfram has indicated that there are machine underflow errors; they provided a workaround, but it is quite slow for my system), but ListConvolve on the same range of t works, and in fact is quite fast. Let's define the model, the instrument function, and some faux data so that we have something to work with:

F[t_, Finf_, A1_, k1_, A2_, k2_, t0_] = Finf - A1 - A2 + 
      UnitStep[t - t0] 
      (A1 + A2 - A1 E^(-k1 (t - t0)) - A2 E^(-k2 (t - t0)));
dn = {0.336025, 0.441503, 0.11445, 0.0549757, 0.0270152, 0.0132802,
      0.00652836, 0.00320924, 0.00157762, 0.000775533, 0.00038124, 
      0.000187412, 0.000092129};
tlist = Range[0, 600, 3];
data = ListConvolve[dn, F[tlist, 4, 2, 0.3, 1.5, 0.03, 50], {1, 1}, 
       0.5] + RandomVariate[NormalDistribution[0, 0.02], Length[tlist]];

Here is the instrument function (dn) vs. point number:

enter image description here

Here is the data:

ListPlot[{tlist, data} // Transpose, PlotRange -> Full]

enter image description here

Now I generate the model points at the experimental times. This is a function of the variable parameters I want to optimize, and involves a ListConvolve.

modelpoints[Finf_, A1_, k1_, A2_, k2_, t0_] = 
    ListConvolve[dn, F[tlist, Finf, A1, k1, A2, k2, t0], {1, 1}, 
    Finf - A1 - A2];

Thus, modelpoints returns a list of modeled, convolved y-values at the same time points as the data. Picking parameters slightly different so that the fitting process has somewhere to go:

initGuess = {Finf -> 3.9, A1 -> 2.1, k1 -> 0.2, A2 -> 1.4, k2 -> 0.04,
             t0 -> 51};
Show[ListPlot[{tlist, data} // Transpose, PlotRange -> Full], 
     ListLinePlot[{tlist, modelpoints[Finf, A1, k1, A2, k2, t0] /. initGuess} 
             // Transpose, PlotStyle -> Red, PlotRange -> Full]]

enter image description here

I can certainly go the manual route of calculating chi-squares and minimizing that value, as follows:

chiSq[Finf_, A1_, k1_, A2_, k2_, t0_] = 
      Total[(data - modelpoints[Finf, A1, k1, A2, k2, t0])^2];
guessRange = {#, 0.9 (# /. initGuess), 
              1.1 (# /. initGuess)} & /@ {Finf, A1, k1, A2, k2, t0};
NMinimize[chiSq[Finf, A1, k1, A2, k2, t0], guessRange]

(* {0.0860335, {Finf -> 3.99996, A1 -> 2.02402, k1 -> 0.311155, 
                A2 -> 1.47918, k2 -> 0.0300756, t0 -> 50.0683}} *)

But it would be wonderful if I could use NonlinearModelFit because of all of its built-in diagnostic functionality. However, I have been unable to find a way to use NonlinearModelFit in a case where the fitting form produces a list of points that should match a list of data points. Can anyone point me to a way to do this? The general idea would be something like:

NonlinearModelFit[data,modelpoints[Finf,A1,k1,A2,k2,t0][[i]], 
                  {#,#/.initGuess}&/@{Finf, A1, k1, A2, k2, t0},i]

Here I am trying (in vain) to use i as the independent variable (of course, it needs to be restricted to positive integers). But this does not work. A comment in this question suggests that NonlinearModelFit can be used if the function is discrete, but doesn't provide any examples. Further it seems that the NominalVariables option is not available for NonlinearModelFit. I have also tried a workaround by turning the model into an InterpolatingFunction, but not only is that really slow, it also gives errors.

$\endgroup$
5
+50
$\begingroup$

You may use SparseArray with Dot to fit. SparseArray gives a warning but it fits almost immediately on my slow laptop.

With all symbols as defined in OP, except SetDelayed on F and modelPoints instead of Set, then

nlm =
 NonlinearModelFit[
  Transpose@{tlist, data},
  modelpoints[Finf, A1, k1, A2, k2, t0].SparseArray[{Floor[i/3 + 1] -> 1}, Length@data],
  {#, # /. initGuess} & /@ {Finf, A1, k1, A2, k2, t0},
  i]

gives a FittedModel object. A little extra calculation was needed due to the step size of 3 starting at zero.

SparseArray continues to complain when properties are referenced but the values are returned.

nlm["BestFitParameters"]
{Finf -> 3.99836, A1 -> 2.06751, k1 -> 0.255743, A2 -> 1.42911, k2 -> 0.0289935, t0 -> 49.7843}
nlm["AdjustedRSquared"]
0.999966

From a plot of the fit (purple) the R-squared seems justified.

Show[
 ListPlot[{tlist, data} // Transpose, PlotRange -> Full, 
  PlotStyle -> LightGray],
 ListLinePlot[{tlist, modelpoints[Finf, A1, k1, A2, k2, t0] /. initGuess} // Transpose,
  PlotStyle -> Directive[Pink, Thin],
  PlotRange -> Full],
 ListLinePlot[{tlist, nlm["Function"] /@ (3 Range[0, 200])} // Transpose,
  PlotStyle -> Purple,
  PlotRange -> Full]
 ]

Mathematica graphics

While this does fit it is inefficient because it calculates the full vector for each value in order to fit. Also, SparseArray constantly complains. I'm open to any ideas to improve these points.

Hope this helps.

$\endgroup$
  • $\begingroup$ Nice! This works beautifully! I'm a little baffled by why it is necessary to use SetDelayed for F and modelPoints, but you seem to be correct that this is necessary. I agree with your self-criticism concerning the calculation of the full vector, and would welcome any improvements that anyone can make in that regard. At any rate, thank you for this solution! $\endgroup$ – Kevin Ausman Jun 10 at 23:21
5
$\begingroup$

Another option would be to convert your ListConvolve discrete model back to a continuous model with Interpolation.

F[t_, Finf_, A1_, k1_, A2_, k2_, t0_] := 
  Finf - A1 - A2 + 
   UnitStep[
     t - t0]*(A1 + A2 - A1 E^(-k1 (t - t0)) - A2 E^(-k2 (t - t0)));
dn = {0.336025, 0.441503, 0.11445, 0.0549757, 0.0270152, 0.0132802, 
   0.00652836, 0.00320924, 0.00157762, 0.000775533, 0.00038124, 
   0.000187412, 0.000092129};
tlist = Range[0, 600, 3];
data = ListConvolve[dn, F[tlist, 4, 2, 0.3, 1.5, 0.03, 50], {1, 1}, 
    0.5] + RandomVariate[NormalDistribution[0, 0.02], Length[tlist]];
initGuess = {Finf -> 3.9, A1 -> 2.1, k1 -> 0.2, A2 -> 1.4, k2 -> 0.04,
    t0 -> 51};
tdata = Transpose@{tlist, data};
lp = ListPlot[tdata, PlotRange -> Full, PlotLegends -> {"Data"}];
(* Create Interpolation Function on ListConvolved Data *)
FI[Finf_, A1_, k1_, A2_, k2_, t0_] := 
 Interpolation[
  Transpose@{tlist, 
    ListConvolve[dn, F[tlist, Finf, A1, k1, A2, k2, t0], {1, 1}, 
     0.5]}, InterpolationOrder -> 1]
nlm = NonlinearModelFit[tdata, FI[Finf, A1, k1, A2, k2, t0][t], 
   List @@@ initGuess, t, Method -> NMinimize];
fit = nlm["BestFit"];
Show[{lp, 
  Plot[fit, {t, 0.`, 600.`}, PlotStyle -> Red, 
   PlotLegends -> {"Fitted"}, PlotRange -> Full]}]
nlm["BestFitParameters"]
(*{Finf -> 3.9973407162246475, A1 -> 1.9841090792021592, k1 -> 3.185244087627753, 
 A2 -> 1.4951069600368265, k2 -> 0.032656509010415835, t0 -> 53.24451084538496} *)

Fitted data

Update Concerning Speed-up

My belief is that specifying Method->NMininmize turns the problem into an unconstrained global optimization problem. I was able to achieve a speed-up of about 3.5x by specifying some of the constrained methods such as NelderMead or SimulatedAnnealing.

{time, nlm} = 
  AbsoluteTiming@
   NonlinearModelFit[tdata, FI[Finf, A1, k1, A2, k2, t0][t], 
    List @@@ initGuess, t, Method -> NMinimize];
{timenm, nlmnm} = 
  AbsoluteTiming@
   NonlinearModelFit[tdata, FI[Finf, A1, k1, A2, k2, t0][t], 
    List @@@ initGuess, t, 
    Method -> {NMinimize, Method -> {"NelderMead"}}];
{timesa, nlmsa} = 
  AbsoluteTiming@
   NonlinearModelFit[tdata, FI[Finf, A1, k1, A2, k2, t0][t], 
    List @@@ initGuess, t, 
    Method -> {NMinimize, Method -> {"SimulatedAnnealing"}}];
time/timenm (* 3.6941030021734855` *)
time/timesa (* 3.4563409868041393` *)

I added some options to the SimulatedAnnealing that seem to speed up the process without having much effect on the fit. It was about 7x faster (variable due to the stochastic nature of SA) and took about 5.25 seconds on my machine.

{timesa, nlmsa} = 
  AbsoluteTiming@
   NonlinearModelFit[tdata, FI[Finf, A1, k1, A2, k2, t0][t], 
    List @@@ initGuess, t, 
    Method -> {NMinimize, 
      Method -> {"SimulatedAnnealing", "PerturbationScale" -> 0.5, 
        "SearchPoints" -> 2}}];
fit = nlmsa["BestFit"];
Show[{lp, 
  Plot[fit, {t, 0.`, 600.`}, PlotStyle -> Red, 
   PlotLegends -> {"Fitted"}, PlotRange -> Full]}]
nlmsa["BestFitParameters"] 
timesa(* 5.257473681307033` *)

Simulated Annealing With Options

$\endgroup$
  • $\begingroup$ Add to your code a definition for tlist, data, initGuess. Otherwise, it looks like a non-working code. $\endgroup$ – Alex Trounev Jun 9 at 6:46
  • $\begingroup$ @AlexTrounev I consolidate the code from above for clarity. Sorry about that. $\endgroup$ – Tim Laska Jun 9 at 11:19
  • $\begingroup$ Is there any way to speed up Interpolation? I think that is what is making this solution so slow. $\endgroup$ – Edmund Jun 10 at 10:14
  • 1
    $\begingroup$ @Edmund I added an update to the post concerning your question. I was able to speed up the fit about 3.5x by specifying a constrained optimization method for NMinimize. Probably, if I had formulated constraints from the beginning, NMinimize would have selected a faster constrained method. $\endgroup$ – Tim Laska Jun 10 at 12:57
  • $\begingroup$ Huh. I had tried something similar, but clearly had made a mistake, because your approach works quite well, and I am definitely impressed by the time improvements you achieved by using SimulatedAnnealing. It doesn't appear, though, that this approach is faster than the method @Edmund used, so in practice I will probably use that approach. It is nice to see how to make this approach work, though! $\endgroup$ – Kevin Ausman Jun 10 at 23:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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