1
$\begingroup$

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?

| improve this question | |
$\endgroup$
  • 1
    $\begingroup$ It looks like you are trying to perform simultaneous fitting. Check these threads: (1), (2). $\endgroup$ – Alexey Popkov Jun 28 '16 at 19:45
  • $\begingroup$ 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. $\endgroup$ – John Doty Jun 28 '16 at 19:59
3
$\begingroup$

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}]

enter image description here

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}]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Doing it that way loses all the handy machinery behind the FittedModel object that NonlinearModelFit yields. $\endgroup$ – John Doty Jun 28 '16 at 19:50
  • 1
    $\begingroup$ see edit, although I think I just reinvented one of the linked answers. Maybe this question should be marked as a dup. $\endgroup$ – george2079 Jun 28 '16 at 20:09
  • $\begingroup$ 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. $\endgroup$ – John Doty Jun 28 '16 at 20:25
  • $\begingroup$ see the last edit, its not so bad if you build up the expression programmatically. $\endgroup$ – george2079 Jun 28 '16 at 21:01
  • $\begingroup$ 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 ツ $\endgroup$ – John Doty Jun 28 '16 at 21:44
2
$\begingroup$

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) *)
| improve this answer | |
$\endgroup$

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.