5
$\begingroup$

So, where to begin. I have a certain dataset that I want to fit to a function. However, the function is rather complicated (matrix inversions, linearsolve) and does in general not have a closed form expression. I understand that this is quite a daunting task, but I'm going to attempt to go through with this, and see what is possible.

To begin with, I suppose I should set up a working example. The function is as follows

GammaEqCTFit[\[Omega]_?NumericQ, LL_, CC_, RR_, CCJ_, CC\[Kappa]_] := 
 Module[{ain = 1., K, 
   C\[Kappa]t, \[Kappa]1, \[Gamma]1, \[Gamma]2, \[Gamma]3, MC, 
   iC, \[Omega]1, \[Omega]2, \[Omega]3, J12, J23, J13, Dm, Av, r, Zc, 
   L1, L2, L3, CJ1, CJ2, CJ3, R1, R2, R3, C\[Kappa], C1, C2, C3},
  Zc = 50;
  L1 = LL; L2 = LL; L3 = LL;
  CJ1 = CCJ; CJ2 = CCJ; CJ3 = CCJ;
  R1 = RR; R2 = RR; R3 = RR;
  C\[Kappa] = CC\[Kappa];
  C1 = CC; C2 = CC; C3 = CC;
  K = 1 + (1/(Zc C\[Kappa] \[Omega]))^2; 
  C\[Kappa]t = C\[Kappa] (1 - 1/K);

  MC = ( {
     {C1 + C\[Kappa]t + CJ1 + CJ3, -CJ1, -CJ3},
     {-CJ1, C2 + CJ1 + CJ2, -CJ2},
     {-CJ3, -CJ2, C3 + CJ2 + CJ3}
    } );
  iC = Inverse[MC];

  \[Kappa]1 = 1/(K Zc ) iC[[1, 1]]; \[Gamma]1 = 
   1/R1  iC[[1, 1]]; \[Gamma]2 = 1/R2  iC[[2, 2]]; \[Gamma]3 = 
   1/R3  iC[[3, 3]];
  \[Omega]1 = Sqrt[1/L1 iC[[1, 1]]]; \[Omega]2 = Sqrt[
   1/L2 iC[[2, 2]]]; \[Omega]3 = Sqrt[1/L3 iC[[3, 3]]];
  J12 = 1/2 iC[[1, 2]]/Sqrt[iC[[1, 1]] iC[[2, 2]]]
     Sqrt[\[Omega]1 \[Omega]2];
  J23 = 1/2 iC[[2, 3]]/Sqrt[iC[[2, 2]] iC[[3, 3]]]
     Sqrt[\[Omega]2 \[Omega]3];
  J13 = 1/2 iC[[1, 3]]/Sqrt[iC[[1, 1]] iC[[3, 3]]]
     Sqrt[\[Omega]1 \[Omega]3];
  Dm = \!\(\*
TagBox[
RowBox[{"(", "", GridBox[{
{
RowBox[{"\[Omega]1", "-", "\[Omega]", "+", 
FractionBox[
RowBox[{"I", " ", "\[Gamma]1"}], "2"], "+", 
FractionBox[
RowBox[{"I", " ", "\[Kappa]1"}], "2"]}], "J12", "J13"},
{"J12", 
RowBox[{"\[Omega]2", "-", "\[Omega]", "+", 
FractionBox[
RowBox[{"I", " ", "\[Gamma]2"}], "2"]}], "J23"},
{"J13", "J23", 
RowBox[{"\[Omega]3", "-", "\[Omega]", "+", 
FractionBox[
RowBox[{"I", " ", "\[Gamma]3"}], "2"]}]}
},
GridBoxAlignment->{
         "Columns" -> {{Center}}, "ColumnsIndexed" -> {}, 
          "Rows" -> {{Baseline}}, "RowsIndexed" -> {}, "Items" -> {}, 
          "ItemsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {
Offset[0.27999999999999997`], {
Offset[0.7]}, 
Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
Offset[0.2], {
Offset[0.4]}, 
Offset[0.2]}, "RowsIndexed" -> {}, "Items" -> {}, 
          "ItemsIndexed" -> {}}], "", ")"}],
Function[BoxForm`e$, 
MatrixForm[BoxForm`e$]]]\);
  Av = I Sqrt[{\[Kappa]1, 0, 0}] ain ;
  r = LinearSolve[Dm, Av];
  \[CapitalGamma] = 1 - Sqrt[{\[Kappa]1, 0, 0}].r/ain;
  Return[\[CapitalGamma]]
  ]

it looks very ugly when copy pasting, here is a picture instead (right click view for a bigger version). here

Alright, so, this beast is fully functional, but clearly not all that fast. I don't see an obvious way of speeding it up however, so, yeah. Next, we should perhaps generate a sample data set:

Freq = Range[6, 6.2, 0.005];
S11 = Table[
   GammaEqCTFit[2*Pi*Freq[[n]], 1.74, 300*10^-6, 20000, 30*10^-6, 
    1*10^-6], {n, 1, Length[Freq]}];
AbsS11 = Abs[S11];
DataTable = Table[{Freq[[n]], AbsS11[[n]]}, {n, 1, Length[Freq]}];
Export["CTTestData.dat", DataTable, "Table"]
ListPlot[Table[{Freq[[n]], AbsS11[[n]]}, {n, 1, Length[Freq]}], 
 PlotRange -> All, 
 AxesLabel -> {"Freq (Ghz)", "Abs (\[CapitalGamma])"}]

which looks like this

enter image description here

I then load the data, and use the most obvious fitting function I could find

AbsS11Data = SetPrecision[Import["CTTestData.dat", "Table"], 25];
nlm = NonlinearModelFit[AbsS11Data, 
  Abs[GammaEqCTFit[2*Pi*\[Omega], 1.74, 300*10^-6, 20000, 30*10^-6, 
    CC\[Kappa]]], {{CC\[Kappa], 5*10^-6}}, \[Omega], 
  WorkingPrecision -> 25]

I've changed this part, incorporating the comments to this post. This example works, and gives the correct output. A problem is however that once I leave more than 1 parameter open, it again stops functioning. Say we take

nlm = NonlinearModelFit[AbsS11Data, 
  Abs[GammaEqCTFit[2*Pi*\[Omega], 1.74, 300*10^-6, 20000, CCJ, 
    CC\[Kappa]]], {{CCJ, 10*10^-6}, {CC\[Kappa], 5*10^-6}}, \[Omega], 
  WorkingPrecision -> 25]

This no longer gives correct answers. It produces the message "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."

I've tried setting these parameters to higher values, like 25, but to no avail. So I'm not sure what to do from this point on.

I should add that I've looked at a few other posts on this topic, such as Problem with NonlinearModelFit and some others. However, not only do the suggestions there not improve the result (using differential evolution and such), the problem there is often that no initial value is given. Now, I'm cheating here and putting it very close to the actual initial value, and even then it seems to fail, which is why I opened a new question.

$\endgroup$
  • $\begingroup$ Sorry to any readers before the edit, NonlinearModelFit can be used without a close form expression, that was my bad. $\endgroup$ – user129412 Jun 12 '15 at 15:49
  • $\begingroup$ You should post here the InputForm of your expressions (no Boxes) $\endgroup$ – Dr. belisarius Jun 12 '15 at 16:10
  • $\begingroup$ Hm, I'm not sure how that works. If I put InputForm around the expression, it returns null? $\endgroup$ – user129412 Jun 12 '15 at 16:32
  • $\begingroup$ Thinking about the issue for a little while, I should probably take a different route. Maybe I can fit the J, omega, kappa and gamma values (no matrix inversion required) and work back from there. $\endgroup$ – user129412 Jun 12 '15 at 17:06
  • 1
    $\begingroup$ You will have to set the input data to an arbitrarily higher precision as well; you could use AbsS11Data = SetPrecision[Import["CTTestData.dat", "Table"], 25];. You should also set the precision of the numerical parameters you fed to the model function, e.g. 1.74`25; and of the initial guess: 0.5`25*10^-6, and then you add WorkingPrecision -> 25 as an option to NonlinearModelFit. $\endgroup$ – MarcoB Jun 12 '15 at 17:49
3
$\begingroup$

I think the issue is indeed about precision. And just one main change seems to make it all work. The line

Module[{ain = 1., K, 

needs to be changed to

Module[{ain = 1, K, 

Otherwise the function GammaEqCTFit does not return a number with enough precision. Also, the parameters for GammaEqCTFit should be entered as rational numbers to keep the high precision. Below I've added in a bit of random noise to be a more challenging test:

Freq = Range[6, 62/10, 5/1000];
AbsS11 = Abs[
   Table[GammaEqCTFit[2*Pi*Freq[[n]], 174/100, 300*10^-6, 20000, 30*10^-6, 1*10^-6],
    {n, 1, Length[Freq]}]];
(* Add in some random noise *)
SeedRandom[12345];
AbsS11 = AbsS11 + RandomVariate[NormalDistribution[0, 1/30000], Length[Freq], 
  WorkingPrecision -> 25];
data = Transpose[{Freq, AbsS11}];

nlm = NonlinearModelFit[data, Abs[GammaEqCTFit[2*Pi*ω, 174/100, 300*10^-6, 20000, CCJ, CCκ]],
  {{CCJ, 32*10^-6}, {CCκ, 8^-7}}, ω, WorkingPrecision -> 25];
nlm["BestFitParameters"]
(* {CCJ -> 0.00002999980726002366869708664, CCκ -> 9.983725353089902821210579*10^-7} *)

Show[ListPlot[data, PlotRange -> All, 
  AxesLabel -> {"Freq (Ghz)", "Abs (Γ)"}],
 Plot[nlm[x], {x, 6, 6.2}, PlotRange -> All]]

Curve fit and data

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.