NonlinearModelFit with indexed parameters

I have an instrument calibration problem for a multi-channel instrument. Some of the calibration parameters are common between channels, while others are per-channel. Toy example:

data = {{1, 0.5, 1.1}, {1, 1.2, 2.2}, {2, 0.3, 0.6}, {2, 2.2, 5.5}}

Column 1 is the (integer) channel number. Now if I try:

NonlinearModelFit[data, a[n] Exp[b x], {a[1], a[2], b}, {n, x}]

I get a NonlinearModelFit::nrlnum: error. The indexed variable confuses the machinery. The following works:

avar[vars_, n_?NumericQ] := vars[[Round[n]]]
NonlinearModelFit[data, avar[{a[1], a[2]}, n] Exp[b x], {a[1], a[2], b}, {n, x}]

But this seems just a bit clumsy. Is there a better way?

=================

1

 

It looks like you are trying to perform simultaneous fitting. Check these threads: (1), (2).
– Alexey Popkov
Jun 28 at 19:45

  

 

Rob Sewell’s answer in your thread 1 is essentially what I did that worked. I’m interested in whether there’s a simpler way, for the case of many channels.
– John Doty
Jun 28 at 19:59

=================

2 Answers
2

=================

If I understand correctly you have basically independent data per channel but the fit is “coupled” by having the same parameter b. Here is one approach.

sol = Last@
NMinimize[
Total[(a[#[[1]]] Exp[ b #[[2]] ] – #[[3]])^2 & /@ data],
{a[1], a[2], b}]

{a[1] -> 0.579786, a[2] -> 0.462069, b -> 1.12535}

GraphicsRow@Table[Show[{
ListPlot[Select[data, #[[1]] == c &][[All, 2 ;;]],
PlotStyle -> {PointSize[.1], Red}],
Plot[a[c] Exp[b x] /. sol, {x, 0, 3}]}], {c, 2}]

If that’s what you are after maybe there is a way to get there with NonlinearSolve..

Edit, here is a trick that seems to work:

NonlinearModelFit[data ,
Switch[x, 1., a[1], 2., a[2]] Exp[ b y ] , {a[1], a[2], b}, {x, y}]

and you see you get the same result as in the NMinimize approach:

%[“BestFitParameters”]

{a[1] -> 0.579786, a[2] -> 0.462069, b -> 1.12535}

something like this for multiple channels:

ch = Union@data[[All, 1]]

{1,2}

NonlinearModelFit[data ,
Switch[x, Evaluate[Sequence @@ Flatten[{N@#, a[#]} & /@ ch]]]
Exp[b y ] ,
Evaluate[{Sequence @@ (a[#] & /@ ch), b}], {x, y}]

  

 

Doing it that way loses all the handy machinery behind the FittedModel object that NonlinearModelFit yields.
– John Doty
Jun 28 at 19:50

1

 

see edit, although I think I just reinvented one of the linked answers. Maybe this question should be marked as a dup.
– george2079
Jun 28 at 20:09

  

 

The Switch or similar approaches (like the one I’m using with Round and Part) get rather clumsy when you have many channels. But maybe that’s the best way that exists.
– John Doty
Jun 28 at 20:25

  

 

see the last edit, its not so bad if you build up the expression programmatically.
– george2079
Jun 28 at 21:01

  

 

Looks good. Now the only problem is explaining Evaluate[Sequence @@ Flatten[{N@#, a[#]} & /@ ch]]] to a pack of IDL and Python people. But asking for a solution to that would be out of scope here ツ
– John Doty
Jun 28 at 21:44

Consider this toy example with 4 channels:

(* {channel, x, y} *)
data = {
{1, 1, 0.19027753120368535}, {1, 2, 0.4204740826587965}, {1, 3, 0.40095769865041475},
{1, 4, 0.23634200938719693}, {1, 5, 0.22578075443150813}, {1, 6, 0.37813496845699},
{1, 7, 0.4838646226571333}, {1, 8, 0.5372735805892475}, {1, 9, 0.3303194817904658},
{1, 10, 0.47604460078275945}, {2, 1, 0.36327171621720744}, {2, 2, 0.33367975557748897}, {2, 3, 0.27105972341189977},
{2, 4, 0.554434177667526}, {2, 5, 0.5315562384552447}, {2, 6, 0.7388343020079361},
{2, 7, 0.5590302156326521}, {2, 8, 0.42818337259989014}, {2, 9, 0.7743736078022455},
{2, 10, 0.8404533078473905}, {3, 1, 0.279037237641659}, {3, 2, 0.529158651066043},
{3, 3, 0.5780675357790563}, {3, 4, 0.7656938190688338}, {3, 5, 0.546753332962597},
{3, 6, 0.7124116571206809}, {3, 7, 0.9849439289565503}, {3, 8, 0.8980520707221048},
{3, 9, 0.980035007948905}, {3, 10, 1.004441552438515}, {4, 1, 0.9501529392948554},
{4, 2, 0.9433949139708255}, {4, 3, 1.0706241417300806}, {4, 4, 1.132668643121858},
{4, 5, 1.3889831839563211}, {4, 6, 1.3213678187002478}, {4, 7, 1.56705164788335},
{4, 8, 1.7701910964370222}, {4, 9, 2.0352693167578426}, {4, 10, 2.1176549936039266}};

If the model can be linearized (which is not the case for the question above), then the NominalVariables and IncludeConstantBasis options can be used with LinearModelFit:

lm = LinearModelFit[data, {ch, x}, {ch, x}, NominalVariables -> ch,
IncludeConstantBasis -> False];
lm // Normal
(* 0.07105395604248158 x
-0.022849825172828374 DiscreteIndicator[ch,1,{1,2,3,4}]+
0.14869088348829945 DiscreteIndicator[ch,2,{1,2,3,4}]+
0.33706272113684577 DiscreteIndicator[ch,3,{1,2,3,4}]+
1.0389391113119841 DiscreteIndicator[ch,4,{1,2,3,4}] *)
lm[“BestFitParameters”]
(* {-0.022849825172828374, 0.14869088348829945,
0.33706272113684577, 1.0389391113119841, 0.07105395604248158} *)

A nonlinear approach is to create dummy variables:

(* Create dummy variables for each channel *)
data2 = Table[Flatten[
{Table[DiscreteIndicator[data[[i, 1]], j, Range[1, 4]], {j, 4}],
data[[i, {2, 3}]]}, 1],
{i, Length[data[[All, 1]]]}]
(* {{1,1,0.22989224889107301},{1,2,0.07861933963986154},1,3,0.4089512564643703},
{1,4,0.26494490181995184},{1,5,0.23084921113958554}, {1,6,0.40770619997939406},
{1,7,0.5187790893182533},{1,8,0.5772050391967247},{1,9,0.3977081849079005},
{1,10,0.6704837525540908},{2,1,0.5163497637118444},{2,2,0.3265405025883527},
{2,3,0.2853821005295374},{2,4,0.37084364350305604},{2,5,0.7712614412738392},
{2,6,0.48827568677332767},{2,7,0.6778218860461915},{2,8,0.5938059733191623},
{2,9,0.6820028950182039},{2,10,0.9044095994342638},{3,1,0.2955159160404255},
{3,2,0.5282493625208143},{3,3,0.42845419636264376},{3,4,0.7218240641305156},
{3,5,0.59547766033399},{3,6,0.7227340557713245},{3,7,0.8060361002558039},
{3,8,0.9680116790953004},{3,9,0.8313974333929903},{3,10,1.1225949195098706},
{4,1,0.7020888353902941},{4,2,0.8193550748524754},{4,3,1.206425919729158},
{4,4,0.9896977126579931},{4,5,1.461727863156306},{4,6,1.4379027623315819},
{4,7,1.5563304057362273},{4,8,1.997674273351458},{4,9,1.8489958482998055},
{4,10,2.134463575861322}} *)

nlm = NonlinearModelFit[
data2, ({a1, a2, a3, a4}.{ch1, ch2, ch3, ch4}) Exp[b x],
{a1, a2, a3, a4, b},
{ch1, ch2, ch3, ch4, x}];
nlm // Normal
(* (0.2021089728245912 ch1 + 0.30462023006749483 ch2 +
0.41223483223929674 ch3 + 0.8096828300770483 ch4)
E^(0.09673863260484315 x) *)