1
$\begingroup$

I need to fit many parameters to several data sets, where a few parameters are common to all data sets. I could do it in two steps, but it is better to perform the fit in a single step, since NonlinearModelFit cannot deal with the covariance matrix that would be required in the second step. A similar question was asked by April 13, 2016 – List argument for NonlinearModelFit, but somehow it didn't work for me. Boiling this problem to the minimal code, my function and example data could be

f[a_, b_][i_, x_] := a + b[[i]] x

data = Flatten[ Table[
{i,xd,RandomVariate[NormalDistribution[f[1.5, {-0.8, 1.2}][i, xd],0.5]]},
 {i,2}, {xd, 3}] , 1]

When I try

fit1 = NonlinearModelFit[data,f[a, {b1, b2}][i, x],
                         {{a, 1.1}, {b1, 0.6}, {b2, 2.2}}, {i, x}]

NonlinearModelFit complains that "the expression i cannot be used as a part specification" and other similar warnings many times, but finds the correct result, as well as

fit2 = NonlinearModelFit[data, f[a, b][i, x], {{a, 1.1}, {b, {0.6, 2.2}}}, {i, x}]

However, fit1["ParameterTable"] gives the correct table, but not fit2["ParameterTable"] – the list of complaints is varied.

I tried to define an adapter function that takes only scalar parameters that are caught in a List, which is then forwarded to function f[][] inside a Module[], but the complaints from NonlinearModelFit are very similar. Taking out the SubValues from the function definition didn't change anything.

Even if the results are correct, I'm afraid of extending this method to a large problem, truly nonlinear, without understanding what is wrong or what I'm missing.

$\endgroup$
  • 1
    $\begingroup$ To get this done in one step one needs to create dummy variables as @LouisB describes. This can be made more error-free by writing a routine to create the dummy variables rather than re-coding each time. However, one potential problem with the currently available options in NonlinearModelFit is that you have to be able (not just willing) to assume that the error variances are identical (or of a known ratio to each other) for each curve to be fit. $\endgroup$ – JimB May 19 '17 at 2:05
6
$\begingroup$

This works with the example data and function f. It could be generalized to more than 2 curves.

model = a0 + b1 (2 - i) x - b2 (1 - i) x;
nlm = NonlinearModelFit[data, model, {a0, b1, b2}, {i, x}];

{a, b} = {a0, {b1, b2}} /. nlm["BestFitParameters"];

Plot[{f[a, b][1, x], f[a, b][2, x]}, {x, 1, 3},
 Epilog -> {PointSize[1/50], Red, Point[Rest /@ data]}]

The idea is to give NonlinearModelFit the model and parameters it wants and then put the best fit parameters into the form we want.

Generalization of the Method

Here is one method of generalizing to more than 2 curves. In this example we generate 4 straight line "curves" with 10 points each. We use random slopes to generate our data curves. Next we generate the parameters and the model. Note that in this model we use Boole instead of the cruder method shown above. Fitting the data and generating the graph are as above. Here is the code and a sample of the results.

ClearAll["Global`*"]

nCurves = 4; npts = 10;

f[a_, b_][i_, x_] := a + b[[i]] x

slopes = RandomInteger[{-30, 30}, nCurves]/10.0;
data = Flatten[Table[{i, xd, RandomVariate[
      NormalDistribution[f[1.5, slopes][i, xd], 0.3]]}, 
     {i, nCurves}, {xd, npts}], 1];
{ymin, ymax} = {Min@data, Max@data};

bparams = ToExpression /@ Table["b" <> ToString[k], {k, 1, nCurves}];
model = a0 + bparams.Table[Boole[i == k], {k, nCurves}] x;
params = Append[bparams, a0];

nlm = NonlinearModelFit[data, model, params, {i, x}];

{a, b} = {a0, bparams} /. nlm["BestFitParameters"];

Plot[Evaluate@Table[f[a, b][k, x], {k, nCurves}], {x, 0, npts},
 PlotRange -> {All, {ymin, ymax}},
 Epilog -> {PointSize[1/50], Red, Point[Rest /@ data]}]

(*  Experimental`NumericalFunction::dimsl: {x} given in {i,x} should be a list of dimensions for a particular argument.  *)

We note the "error message", more like a warning, really, was generated using MMA 11.0. The plotted results show the method is properly generalized. To double check, we could re-run the code with a greater number of curves or data points. To do this does not require changing any code, except the value assigned to the variable nCurves or npts in the first line. enter image description here

Alternate Generalization

Here is an alternate generalization that uses the original "i-k" method of parameter selection instead of Boole. This method may be cruder, but it avoids the "error message". Simply replace the "model = ..." line in the above generalization with these 4 lines of code:

prod = Product[k - i, {k, nCurves}];
pr2 = Table[prod/(k - i), {k, nCurves}];
coefs = pr2/Table[pr2[[k]] /. i -> k, {k, nCurves}];
model = a0 + (bparams.coefs) x;

In the alternate model we first generate the product of all of the "$k-i$" factors. Then we divide out 1-$i$, 2-$i$, etc to get the unique product for each $b$-parameter. The actual coefficients of the $b$-parameters are the normalized products. Each coefficient is either 0 or 1, for integer $i$ from 1 to nCurves. For the test data, this alternate model gives exactly the same fit, as was intended.

$\endgroup$
  • $\begingroup$ The trouble with your solution is that I'll have to rewrite the model function for the required number of parameters in each case - I'll need a code for 3 b's, 4 b's, ... 22 b's - too much copy and paste will lead to typo errors and mistakes difficult to spot. I would like to have a single model function, so I can be sure that it is correctly coded - the true model function is not a straight line, and I do not have the same number of parameters in all cases. $\endgroup$ – Vito Vanin May 19 '17 at 1:54
  • $\begingroup$ Many thanks! I realized that I had completely missed your point after the comment by Jim Baldwin. I came to my office thinking to build the function as bparams.Join[ ConstantArray[0,i-1],{1},ConstantArray[0,nCurves-i-1] ], but your solution is very readable. I will fit your solution to my problem, check a few things, and comment afterwards. $\endgroup$ – Vito Vanin May 19 '17 at 11:49
  • $\begingroup$ I've already checked your method. It works and is easily adapted to different numbers of parameters and datasets. I have understood that one of the key issues is the definition of a "model function" that uses local variables as arguments, and it is impossible to do this thing using a function defined with SetDelayed - is this correct? Many thanks for your nice help. $\endgroup$ – Vito Vanin May 19 '17 at 12:46
  • $\begingroup$ @VitoVanin We can define, say, model2[x_,b_List] := First[b] + Last[b] x using SetDelayed and then evaluate b={a0,a1}; NonlinearModelFit[{1.1,1.8,3.1}, model2[x,b], b, x] to obtain a fit. This is very similar to your original approach. The significant difference may be that, to MMA, your 'fit[a,b][i,x]' used the continuous variable $i$ as the index, which must take discrete values. $\endgroup$ – LouisB May 19 '17 at 20:59

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.