Simultaneous plotting, Nonlinear model fit is not converging for certain data sets

Question

I had asked a question regarding simultaneous plotting here.

My Program will fit real part and imaginary part of a data simultaneously and gives out the value something like this

I column has frequency, II column has Amplitud, III column has Phase.

Here is my modified program from the responses I got from previous question.

Qfit[data1_, fmin_, fmax_] := Module[{da1 = {}, da2 = {}, L, A, B, kint, kport, f, f0, index, mini,
res, i, r, myF, nlm, data = {}, allData}, Clear[A, B, kint, kport, f, f0, index, mini, res, L]; data = Select[data1, fmin <= #[[1]] <= fmax &];da1 = Transpose[{data[[All, 1]], 
 data[[All, 2]]*Cos[data[[All, 3]]]}];da2 = Transpose[{data[[All, 
  1]], (data[[All, 2]]*Sin[data[[All, 3]]])}]; mini = Flatten[Position[da1[[All, 2]], Min[da1[[All, 2]]]]][[1]];res = da1[[mini, 1]];allData =Join[{1, Sequence @@ #} & /@ da1, {2, Sequence @@ #} & /@ da2]; i[f_] := A*B*(4 (f - f0)^2 + (kint - kport) (kint + kport))/(4 (f - f0)^2 + (kint + kport)^2);r[f_] := A*(1 - B (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2));myF[index_, f_] := KroneckerDelta[index - 1] i[f] + KroneckerDelta[index - 2] r[f]; nlm = NonlinearModelFit[allData, myF[index, 
 f], {{kint, 0.001}, {kport, 0.001}, {f0, res}, {A, 0.1}, {B, 
  0.05}}, {index, f}, MaxIterations -> Infinity];fitparams = nlm["BestFitParameters"]; Print[fitparams];Set @@@ fitparams;Print["Internal Quality factor (Loss inside cavity)  =  ", f0/kint, "    External Quality factor (Input port 1 loss)  =  ", f0/kport, "    Total loss or Quality factor   =  ", f0/(kint + kport)]; Print["The graphs are ", Show[ListPlot[{da1, da2}], 
Plot[{nlm[1, f], nlm[2, f]}, {f, fmin, fmax}], 
ImageSize -> Large],Show[ListPlot[Transpose@{Last /@ da1, Last /@ da2}, Joined -> True,
  AxesLabel -> {"Re", "Im"}], 
ParametricPlot[{nlm[1, f], nlm[2, f]}, {f, Min@(First /@ da1), 
  Max@(First /@ da1)}, PlotStyle -> Red], PlotRange -> All, 
ImageSize -> Large]]]

This usually work for all most all the values, but for some data like this,here is a data from experiment.

For this we have to compute as follows Qfit[data,10.45,10.475].

It is not converging at all. I do not why it is happening I appreciate any quick responses, and sorry for my ignorance.

Answer 1

This is an extended comment rather than an answer.

I've downloaded the problem data and I believe the trouble comes from the common parameters not having anywhere near the same values between the two subsets of the data. (However, there is a confounding issue described below.)

I took your code and modified it so that it estimates the two regressions separately:

Qfit[data1_, fmin_, fmax_] := 
 Module[{da1 = {}, da2 = {}, L, A, B, kint, kport, f, f0, index, mini,
    res, i, r, myF, nlm, data = {}, allData}, 
  Clear[A, B, kint, kport, f, f0, index, mini, res, L]; 
  data = Select[data1, fmin <= #[[1]] <= fmax &]; 
  da1 = Transpose[{data[[All, 1]], 
     data[[All, 2]]*Cos[data[[All, 3]]]}]; 
  da2 = Transpose[{data[[All, 1]], (data[[All, 2]]*
       Sin[data[[All, 3]]])}]; 
  mini = Flatten[Position[da1[[All, 2]], Min[da1[[All, 2]]]]][[1]]; 
  res = da1[[mini, 1]];
  allData = 
   Join[{1, Sequence @@ #} & /@ da1, {2, Sequence @@ #} & /@ da2]; 
  ListPointPlot3D[allData]; 
  i[f_] := A*
    B*(4 (f - f0)^2 + (kint - kport) (kint + 
          kport))/(4 (f - f0)^2 + (kint + kport)^2); 
  r[f_] := A*(1 - 
      B (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2)); 
  myF[index_, f_] := 
   KroneckerDelta[index - 1] i[f] + KroneckerDelta[index - 2] r[f]; 
  allData1 = Select[allData, #[[1]] == 1 &]; 
  allData2 = Select[allData, #[[1]] == 2 &];
  allData1 = allData1[[All, {2, 3}]];
  allData2 = allData2[[All, {2, 3}]];
  nlm1 = NonlinearModelFit[allData1, 
    myF[1, f], {{kint, 0.001}, {kport, 0.001}, {f0, res}, {A, 
      0.1}, {B, 0.05}}, f, MaxIterations -> Infinity];
  nlm2 = NonlinearModelFit[allData2, 
    myF[2, f], {{kint, 0.001}, {kport, 0.001}, {f0, res}, {A, 
      0.1}, {B, 0.05}}, f, MaxIterations -> Infinity];

  fitparams1 = nlm1["BestFitParameters"];
  fitparams2 = nlm2["BestFitParameters"]; 
  Print[{fitparams1, fitparams2}];
  Print[nlm1["CorrelationMatrix"]];
  Print[nlm2["CorrelationMatrix"]];
  Print["The graphs are ", 
   Show[ListPlot[{da1, da2}], 
    Plot[{nlm1[f], nlm2[f]}, {f, fmin, fmax}], ImageSize -> Large], 
   Show[ListPlot[Transpose@{Last /@ da1, Last /@ da2}, Joined -> True,
      AxesLabel -> {"Re", "Im"}], 
    ParametricPlot[{nlm1[f], nlm2[f]}, {f, Min@(First /@ da1), 
      Max@(First /@ da1)}, PlotStyle -> Red], PlotRange -> All, 
    ImageSize -> Large]]]
Qfit[data, 10.45, 10.6]

While there was a warning about "The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum.", the fits look fine:

But the (cleaned-up) solutions look very different:

sol1 = {kint -> -0.0027306463823735722`, 
  kport -> -0.0013994544277047404`,
  f0 -> 10.508996991971562`, 
  A -> 0.9737627424533896`, 
  B -> 0.6095062568309071`};
sol2 = {kint -> 0.032106337286351976`, 
  kport -> -0.02754995227500653`,
  f0 -> 10.50900015588828`, 
  A -> 0.0036249688207138334`,
  B -> 9.185735216137683`};

Now fitting i[f] we can't estimate both A and B but only their product. When the product A B is examined for both fits, the products are very different.

The confounding factor is that the estimated correlations among the parameter estimators are many times close to +1 or -1 which can cause trouble with the estimation procedure. It's possible that a reformulation of the model with maybe ksum = kint + kport and kdiff = kint - kport might make the estimation process more stable.

If the data sets that do work fine have high correlations among the parameter estimators, then I would highly recommend some sort of reparameterization of the model.

Answer 2

Maybe this is the answer as to why there is instability and why the extreme correlation between some of the parameter estimators.

Look at the limits of the two functions being fit (i[f] and r[f]) as $f \rightarrow \infty$ and $f \rightarrow -\infty$:

Limit[A*B*(4 (f - f0)^2 + (kint - kport) (kint + kport))/(4 (f - f0)^2 + (kint + kport)^2),
  f -> ∞]
(* A B *)

Limit[A*B*(4 (f - f0)^2 + (kint - kport) (kint + kport))/(4 (f - f0)^2 + (kint + kport)^2), 
  f -> -∞]
(* A B *)

Limit[A*(1 - B (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2)), f -> ∞]
(* A *)

Limit[A*(1 - B (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2)), f -> -∞]
(* A *)

If we look at the data, the value for A in the r[f] curve would appear to be zero or very close to zero. That forces B to be very large.

Is there a fix? A re-parameterization might do it. Rather than having A and B we use a and ab:

i[f_] := ab*(4 (f - f0)^2 + (kint - kport) (kint + kport))/(4 (f - f0)^2 + (kint + kport)^2); 
r[f_] := a - ab (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2);

Here's the full updated code:

Qfit[data1_, fmin_, fmax_] := 
 Module[{da1 = {}, da2 = {}, L, A, B, kint, kport, f, f0, index, mini,
    res, i, r, myF, nlm, data = {}, allData}, 
  Clear[A, B, kint, kport, f, f0, index, mini, res, L]; 
  data = Select[data1, fmin <= #[[1]] <= fmax &]; 
  da1 = Transpose[{data[[All, 1]], 
     data[[All, 2]]*Cos[data[[All, 3]]]}]; 
  da2 = Transpose[{data[[All, 1]], (data[[All, 2]]*
       Sin[data[[All, 3]]])}]; 
  mini = Flatten[Position[da1[[All, 2]], Min[da1[[All, 2]]]]][[1]]; 
  res = da1[[mini, 1]]; 
  allData = 
   Join[{1, Sequence @@ #} & /@ da1, {2, Sequence @@ #} & /@ da2]; 
  i[f_] := ab*(4 (f - f0)^2 + (kint - kport) (kint + 
          kport))/(4 (f - f0)^2 + (kint + kport)^2); 
  r[f_] := a - 
    ab (4 kport (f - f0))/(4 (f - f0)^2 + (kint + kport)^2); 
  myF[index_, f_] := 
   KroneckerDelta[index - 1] i[f] + KroneckerDelta[index - 2] r[f]; 
  nlm = NonlinearModelFit[allData, 
    myF[index, 
     f], {{kint, 0.001}, {kport, 0.001}, {f0, res}, {a, 0.1}, {ab, 
      0.05}}, {index, f}, MaxIterations -> Infinity]; 
  fitparams = nlm["BestFitParameters"]; Print[fitparams]; 
  Set @@@ fitparams; 
  Print["Internal Quality factor (Loss inside cavity)  =  ", f0/kint, 
   "    External Quality factor (Input port 1 loss)  =  ", f0/kport, 
   "    Total loss or Quality factor   =  ", f0/(kint + kport)]; 
  Print["The graphs are ", 
   Show[ListPlot[{da1, da2}], 
    Plot[{nlm[1, f], nlm[2, f]}, {f, fmin, fmax}], 
    ImageSize -> Large], 
   Show[ListPlot[Transpose@{Last /@ da1, Last /@ da2}, Joined -> True,
      AxesLabel -> {"Re", "Im"}], 
    ParametricPlot[{nlm[1, f], nlm[2, f]}, {f, Min@(First /@ da1), 
      Max@(First /@ da1)}, PlotStyle -> Red], PlotRange -> All, 
    ImageSize -> Large]]]

Now when we run the analysis we get the following:

Qfit[data, 10.45, 10.475]