0
$\begingroup$

To summarize: My problem is that I have some code that won't terminate even when I use TimeConstrained to force it to terminate after 1 second. My question is whether there is an alternative solution to forcing the evaluation to terminate. I am running Mathematica 11.2 on Windows 10 with an Intel i9 processor overclocked at 4.5Gz and with 64MB RAM.

Specifics: I am using the following code to fit a model to time series data:

system = {
CC'[t] == b M[t] - CC[t] d,
M'[t] == -M[t] a t^B,
CC[0] == data[[1, 2]], M[0] == data[[1, 2]]};

model = ParametricNDSolveValue[system,  CC, {t, First[data][[1]], Last[data] 
[[1]]}, {a, B, b, d}]

fit = FindFit[data, model[a, B, b, d][t], {{a, 0}, {B, 3}, {b, .5}, {d, 
.01}}, t]

This works fine for most of the data I have tested. However, for some data, FindFit never terminates and I have to abort kernel in order to halt the evaluation.

Here are three examples of data that cause the problem:

data = {{0, 4496.5}, {24, 37416.5}, {48, 74462.7}, {72, 57204.8}};
data = {{0, 4500}, {24, 37000}, {48, 74000}, {55, 74000}, {72, 
57000}};
data = {{0, 4500}, {12, 25000}, {16, 29000}, {20, 33000}, {24, 
37000}, {30, 50000}, {35, 65000}, {40, 70000}, {48, 74000}, {55, 
74000}, {60, 73000}, {65, 70000}, {68, 65000}, {72, 57000}};

and here are two examples of data that don't cause a problem:

data = {{0, 5117}, {24, 66183.3}, {48, 119679}, {72, 79002.2}};
data = {{0, 4500}, {24, 37000}, {48, 74000}, {55, 70000}, {72, 
57000}};

(By the way, I am aware that the model with some of these data is being overfit. These are just examples).

The REAL problem is that the problem persists even when using TimeConstrained. I tried adding TimeConstrained to FindFit as follows:

fit = TimeConstrained[FindFit[data, model[a, B, b, d][t], {{a, 0}, {B, 3}, 
{b, .5}, {d, .01}}, t, MaxIterations -> 10], 1]

But this doesn't help: the computation never terminates and the only way to abort the evaluation is to quit the kernal (note that the problem occurs even when MaxIterations is set to only 10).

Since I want to run this code in a loop on a large number of data sets, I would like a way to force Mathematica to time out if the evaluation won't terminate.

Any help is much appreciated in advance!

$\endgroup$

closed as off-topic by Daniel Lichtblau, eyorble, LCarvalho, WReach, JungHwan Min Aug 5 '18 at 1:59

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – Daniel Lichtblau, eyorble, LCarvalho, WReach, JungHwan Min
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ I cannot get the code above to run. $\endgroup$ – Daniel Lichtblau Jul 20 '18 at 15:09
  • $\begingroup$ Sorry Daniel. I had an extra square bracket in my pasted code. I just edited the question and the code should run now. $\endgroup$ – dvd8000 Jul 20 '18 at 16:14
  • $\begingroup$ I might be missing something but why would you try to estimate 5 parameters (a, b, B, d, and the error variance) with just 4 (or even 5) data points? I understand you acknowledge that these are just small examples but larger amounts of data would likely not have the same issues. $\endgroup$ – JimB Jul 20 '18 at 16:30
  • $\begingroup$ The datasets that "work" aren't converging even with MaxIterations -> 1000. If you get warnings or errors, it is essential that you mention that. $\endgroup$ – JimB Jul 20 '18 at 16:39
  • $\begingroup$ Hi JimB, I just added a dataset with 14 points that also has the problem. So, while I agree that this won't be as big an issue with larger data sets, the problem does arise even when data sets are larger. $\endgroup$ – dvd8000 Jul 20 '18 at 16:51
3
$\begingroup$

I think your model is overparameterized (and not because of the amount of data). If one sets a to some value, then FindFit (and better yet NonlinearModelFit) finds a solution. If one tries lots of different values for a and plots those values against the associated estimates of B, then one gets a near-perfect correlation. Because of the near-perfect correlation, there is no unique solution and a single iteration in FindFit just continues to wander around in hyperspace. (I know that's not a great or even accurate explanation.)

So I suggest either (1) Removing a (i.e., setting it to 1) or (2) rethinking your model. Here are the results with setting a to 1. (I've also divide t by 100 to reduce round-off error.)

data = {{0, 4500}, {12, 25000}, {16, 29000}, {20, 33000}, {24, 37000}, 
  {30, 50000}, {35, 65000}, {40, 70000}, {48, 74000}, {55, 74000},
  {60, 73000}, {65, 70000}, {68, 65000}, {72, 57000}};

system = {CC'[t] == b M[t] - CC[t] d,
   M'[t] == -M[t] (t/100)^B,
   CC[0] == data[[1, 2]],
   M[0] == data[[1, 2]]};

model = ParametricNDSolveValue[system, 
   CC, {t, First[data][[1]], Last[data][[1]]}, {b, B, d}];

nlm = NonlinearModelFit[data, 
   model[b, B, d][t], {{b, 1.4}, {B, 3}, {d, .1}}, t];

nlm["CorrelationMatrix"] // MatrixForm

$$\left( \begin{array}{ccc} 1. & 0.918946 & 0.97062 \\ 0.918946 & 1. & 0.983729 \\ 0.97062 & 0.983729 & 1. \\ \end{array} \right)$$

nlm["BestFitParameters"] 
(* {b -> 0.43894, B -> 4.08547, d -> 0.0111581} *)

Show[ListPlot[data], Plot[nlm[z], {z, 0, 72}]]

Data and fit

Update

I think the main issue is still that the model is overparameterized (and that is not to be confused with overfitting). If one re-parameterizes the model such that the coefficients to be estimated are on approximately the same scale, then the fitting algorithms are much more stable. The catch, of course, is having both a good idea as to how to re-parameterize and good starting values. Here is an approach that will fit all 4 parameters: a, b, B, and d.

data = {{0, 4500}, {12, 25000}, {16, 29000}, {20, 33000}, {24, 37000}, {30, 50000},
   {35, 65000}, {40, 70000}, {48, 74000}, {55, 74000}, {60, 73000}, {65, 70000},
   {68, 65000}, {72, 57000}};

system = {CC'[t] == b M[t] - CC[t] d/100,
   M'[t] == -M[t] Exp[loga] (t/100)^B,
   CC[0] == data[[1, 2]],
   M[0] == data[[1, 2]]};

model = ParametricNDSolveValue[system, 
   CC, {t, First[data][[1]], Last[data][[1]]}, {loga, b, B, d}];

nlm = NonlinearModelFit[data, 
   model[loga, b, B, d][t], {{loga, 14}, {b, 0.4}, {B, 22}, {d, 1}}, t];
nlm["BestFitParameters"]
(* {loga -> 13.9864, b -> 0.413836, B -> 21.9799, d -> 0.924626} *)

Show[ListPlot[data], Plot[nlm[x], {x, 0, 72}]]

Data and fit

But note that loga and B are still highly correlated (0.99942) which suggests that the model is either overparameterized or not "well parameterized".

nlm["CorrelationMatrix"] // MatrixForm

$$\left( \begin{array}{cccc} 1. & -0.066721 & 0.99942 & 0.159348 \\ -0.066721 & 1. & -0.0439451 & 0.884898 \\ 0.99942 & -0.0439451 & 1. & 0.189313 \\ 0.159348 & 0.884898 & 0.189313 & 1. \\ \end{array} \right)$$

$\endgroup$

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