Getting Convolve to behave in preparation for NonlinearModelFit

Question

New information added at bottom of post

This is an outgrowth of this earlier question. Let us suppose that we have an instrumental response function given here:

response[t_] = 0.2583*(1/(E^(0.25*(-2.6 + t))*(1 + 
       13.5/E^(1.67*(-2.6 + t)))^10.9) + (961*(1 - 
       E^(-0.236*(-2.6 + t))))/(E^(1.67*(-2.6 + t))*(1 + 
       13.5/E^(1.67*(-2.6 + t)))^11.9));
Plot[response[t], {t, -10, 50}, PlotRange -> Full]

And the underlying physical phenomenon I am expecting to be occurring is given by:

model[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)));

We will take some test parameters just to evaluate the shape of the curve as follows:

testparams = {Finf->4, A1->1.5, A2->2.0, k1->0.3, k2->0.06, t0->50};
Plot[model[t, Finf, A1, k1, A2, k2, t0] /. testparams, {t, 0, 200}, 
     PlotRange -> {{0, 200}, {0, 4}}]

The end goal will be to find the values of the parameters that best fit experimental data. So naturally what I want to do is find the convolution of response and model, preferably as a function although evaluated at a set of x values that match my experimental x axis could work too, and then do a NonlinearModelFit on that convolution in order to find the values of those parameters. I am having a lot of trouble with the convolution, however.

fn[x_, Finf_, A1_, k1_, A2_, k2_, t0_] = 
   Convolve[response[t], model[t, Finf, A1, k1, A2, k2, t0], t, x]

When I try it in MMA 11.3 and plot it using the same test parameters from above, I get this:

When I try it in MMA 12.0 and plot it using the same test parameters from above, I get this:

Zooming in, I see this:

In all cases I get numerous warnings along the lines of: General: [Some fraction] is too small to represent as a normalized machine number; precision may be lost.

That earlier question was (sort of) resolved by multiplying the response function by UnitStep[t]. This does not change the MMA 11.3 version's resulting plot. It does change the MMA 12.0 version's resulting plot:

So it's looking like some kind of edge effect, possibly due to the (slight) discontinuity at t=0 in the response function. So I then changed the UnitStep[t] in the response function to UnitStep[t-3] in order to get it as close as possible to matching at zero, to find:

The edge effects look to be gone (even if I have no idea where they were coming from), but look at the shape of that curve compared to the original model. It has a faster rise to the asymptote than the unconvoluted model, which is clearly wrong.

Finally, if I take the t-3 version but swap the order of two functions in the convolution, we see this:

And again, throughout, I am getting those number too small warnings. Trying to swap the order of the function in the convolution in MMA 11.3, however, seems to hang the kernel.

If I plug in the parameters first before convolution in MMA 12.0 (even without the UnitStep in the response function), I get a result that looks right:

Although with all of the other problems, I have a hard time trusting this result. Plus, having to recalculate the convolution at each set of parameters would be computationally intractable during a NonlinearModelFit. And in MMA 11.3, the results are still wrong.

Both of my functions look well-behaved enough that convolution should be possible. My response function is not all that bizarre. It is an extremely common problem in science to use an instrument response function to help model data, so it is hard to believe that this problem hasn't been come across before. Am I missing something, or doing something wrong?

As a side note, I saw a couple of other posted questions where the response was to use PiecewiseExpand to replace the UnitStep with a piecewise function in order to speed things up, but unfortunately the convolution doesn't seem to execute for me when I try that with my functions.

Edited to add a new test I ran: Ok, this has me completely baffled.

You can visualize a convolution as an integral where you swap the x axis of one of the functions and then look at different x-offsets. So I set it up with a Manipulate and then looked at the integral:

Manipulate[
  Show[Plot[{model[t, Finf, A1, k1, A2, k2, t0] /. testparams, 
             response[-t + offset] 10}, {t, 0, 200}, PlotRange -> Full, 
            PlotPoints -> 200, ImageSize -> Large, 
            LabelStyle -> Directive[Background -> White]], 
       ListPlot[{{offset, 
                  NIntegrate[(model[x, Finf, A1, k1, A2, k2, t0] /. 
                  testparams) response[-x + offset], {x, 0, 200}]}}]], 
       {{offset, 30}, 0, 250}]

Strictly, you are supposed to integrate from -Infinity to Infinity. Here I am restricting the integral to the viewing window of 0-200, but the result is the same.

That makes NO sense at all!

Further information:

Yeah, NIntegrate is simply giving incoherent results. Witness:

mdl[x_] = (model[x, Finf, A1, k1, A2, k2, t0] /. testparams)
rsp[x_] = response[-x + 121]
Plot[{mdl[x], rsp[x], mdl[x] rsp[x]}, {x, 0, 200}, 
      Filling -> {3 -> 0}, Background -> White, 
      LabelStyle -> Background -> White]
NIntegrate[mdl[x] rsp[x], {x, 0, 200}, AccuracyGoal -> 10]

So first I define the functions as just being functions of x. I plot them individually plus their product (with the product being filled to the x axis, since we will be integrating). The functions are well-behaved, and the product is >= zero over the whole integration range (barring a VERY small ringing region in the response function). Yet NIntegrate gives me a result of zero. What in the actual @$!&?

Answer 1

All of the Convolve examples in the documentation appear to rational coefficients. The complexity of your kernel may make it difficult to find analytical solutions in a general way. ListConvolve, however, appears to be straightforward to use if you can make your problem cyclic (I simply reflected your functions about the end time). The following is an example.

First, define the functions:

response[t_] = 
  0.2583*(1/(E^(0.25*(-2.6 + t))*(1 + 
           13.5/E^(1.67*(-2.6 + t)))^10.9) + (961*(1 - 
          E^(-0.236*(-2.6 + t))))/(E^(1.67*(-2.6 + t))*(1 + 
           13.5/E^(1.67*(-2.6 + t)))^11.9));
model[Finf_, A1_, k1_, A2_, k2_, t0_][t_] := 
  Finf - A1 - A2 + 
   UnitStep[
     t - t0] (A1 + A2 - A1 E^(-k1 (t - t0)) - A2 E^(-k2 (t - t0)));
delta[t0_, tf_, n_] := (tf - t0)/(n - 1)
grid[t0_, tf_, n_] := N@Range[t0, tf, delta[t0, tf, n]]

Next, set up some test parameters:

(* Test Parm *)
Finf = 4;
A1 = 1.5;
A2 = 2.0;
k1 = 0.3;
k2 = 0.06;
t0 = 50;

Then, follow this procedure to ListConvolve:

starttime = 0;
finaltime = 200;
npoints = 200;
(* Create the grid *)
g = grid[starttime, finaltime, npoints];
(* Digitize the response/kernel *)
rdigitized = N@response[#] &@ g;
(* Need to make cyclic *)
rdigitized = rdigitized~Join~Reverse[rdigitized];
(* Normalization *)
rdigitized = rdigitized/Total[rdigitized];
(* Digitize the model *)
mdigitized = N@(model[Finf, A1, k1, A2, k2, t0][#]) &@g;
mdigitized =  mdigitized~Join~Reverse[mdigitized];
(* Perform the convolution *)
conv = ListConvolve[rdigitized, mdigitized, 1];
discreteConv = Transpose[{g, conv~Take~npoints}];
(* Display the results *)
Show[ListPlot[discreteConv, PlotLegends -> {"Convolution"}], 
 Plot[model[Finf, A1, k1, A2, k2, t0][t], {t, 0, 200}, 
  PlotRange -> {{0, 200}, {0, 4}}, PlotStyle -> {Red, Dashed}, 
  PlotLegends -> {"Model"}]]

It runs pretty fast and seems to give a reasonable result and you will not be limited to models that have analytical convolutions. You will have to restrict your time range, but that would likely be a lesser headache.

Update

I wanted to test the interactivity of the ListConvolve approach so I wrapped the functionality in module and created a simple Manipulate. It seems be pretty fast up to 1000 digitized points.

Here is the updated code:

response[t_] = 
  0.2583*(1/(E^(0.25*(-2.6 + t))*(1 + 
           13.5/E^(1.67*(-2.6 + t)))^10.9) + (961*(1 - 
          E^(-0.236*(-2.6 + t))))/(E^(1.67*(-2.6 + t))*(1 + 
           13.5/E^(1.67*(-2.6 + t)))^11.9));
model[Finf_, A1_, k1_, A2_, k2_, t0_][t_] := 
  Finf - A1 - A2 + 
   UnitStep[
     t - t0] (A1 + A2 - A1 E^(-k1 (t - t0)) - A2 E^(-k2 (t - t0)));
delta[t0_, tf_, n_] := (tf - t0)/(n - 1)
grid[t0_, tf_, n_] := N@Range[t0, tf, delta[t0, tf, n]]
cyclic = N@(#~Join~Reverse[#] &)@(dummy /@ #) &;
discreteConvolve[starttime_, finaltime_, npoints_, r_, m_] := 
 Module[{gr, g, rdigitized, mdigitized, conv, discreteConv},
  (* Create the cyclic grid *)
  gr = grid[starttime, finaltime, npoints];
  g = cyclic@gr;
  (* Digitize the response/kernel *)
  rdigitized = r @@@ g;
  (* Normalization *)
  rdigitized = rdigitized/Total[rdigitized];
  (* Digitize the model *)
  mdigitized = m @@@ g;
  (* Perform the convolution *)
  conv = ListConvolve[rdigitized, mdigitized, 1];
  discreteConv = Transpose[{gr, conv~Take~npoints}];
  (* Display the results *)
  discreteConv]
Manipulate[
 Show[ListPlot[
   discreteConvolve[0, finaltime, npoints, response, 
    model[Finf, A1, k1, A2, k2, t0]], PlotLegends -> {"Convolution"}],
   Plot[model[Finf, A1, k1, A2, k2, t0][t], {t, 0, finaltime}, 
   PlotRange -> {{0, 200}, {0, 4}}, PlotStyle -> {Red, Dashed}, 
   PlotLegends -> {"Model"}]],
 {{Finf, 4}, 2, 6},
 {{A1, 1.5}, 1, 3},
 {{A2, 2}, 1, 3},
 {{k1, 0.3}, 0.1, 0.5},
 {{k2, 0.06}, 0.03, 0.09},
 {{t0, 50}, 0, 100},
 {{finaltime, 200}, 200, 400},
 {{npoints, 200}, 100, 1000, 100}
 ]

Here is a Manipulate animation (speed up 2x):

Update 2

It turns out that ListConvolve does not like my cyclic implementation for anything but even functions so I simply replaced it with:

cyclic = N@(# &)@(dummy /@ #) &;

I also found the maximal overhangs and minimal padding {1,-1},0 seemed to have the least amount of start-up effect or by replace conv with:

conv = ListConvolve[rdigitized, mdigitized, {1, -1}, 0];

There is a start up effect that you need to consider, but once beyond seems to change the solution little beyond a simple shift.

Answer 2

Wolfram Technical Support has weighed in. The problem is a machine underflow error, and remains under investigation. A workaround can be implemented by using Rationalize and increasing WorkingPrecision (although I had to go up to WorkingPrecision->50 to get good results using the example above in MMA 11.3):

Plot[Rationalize[fn[t, Finf, A1, k1, A2, k2, t0] /. testparams, 0] // Evaluate, 
     {t, 0, 200}, PlotRange -> {{0, 200}, {0, 4}}, WorkingPrecision -> 50] // Quiet

In practice, an answer based on @TimLaska's answer looks to be quicker (keep in mind reading his answer, though, that the padding/overlap has to be handled differently; read the comment thread to his answer for details), although I do lament the seeming inability to use the ListConvolve approach with NonlinearModelFit, as I discussed in this question.