List Argument for NonlinearModelFit - another version

Question

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.

Answer 1

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.

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.