2
$\begingroup$

I am trying to fit the experimental data to the following equation: (2*π/377)*Sum[S[i]^2*Y[i]*w^2/(W0[i]^2 - w^2)^2 + Y[i]^2*w^2), {i, 1, n}].

The data I have is for n=4 with W01=0 and so far, I've done the following:

data = Import["https://pastebin.com/raw/SDzdNXrd", "Table"]
model = (2 π/377) (S1^2*Y1*w^2/((0 - w^2)^2 + Y1^2*w^2) + 
 S2^2*Y2*w^2/((W02^2 - w^2)^2 + Y2^2*w^2) + 
 S3^2*Y3*w^2/((W03^2 - w^2)^2 + Y3^2*w^2) + 
 S4^2*Y4*w^2/((W04^2 - w^2)^2 + Y4^2*w^2));
fit = FindFit[data, {model, {S1 > 0, S2 > 0, S3 > 0, S4 > 0, W02 > 0, W03 > 0,
                             W04 > 0, Y1 > 0, Y2 > 0, Y3 > 0, Y4 > 0}},
              {S1, S2, S3, S4, W02, W03, W04, Y1, Y2, Y3, Y4}, w,
              AccuracyGoal -> 5, MaxIterations -> 10000]
Show[Plot[Evaluate[model /. fit], {w, 0, 60000}, AxesOrigin -> {0, 0}, PlotRange -> All,
     AxesLabel -> {"Frequency", "Conductivity"}], ListPlot[data]]

And I get an extremely poor fit with this.

I am very new to Mathematica so any help would be appreciated.

$\endgroup$
5
  • $\begingroup$ Have you tried specifying better starting values for your parameters? You can use the alternative syntax FindFit[data, model, {{par1, startval1}, {par2, startval2}, ...}, vars]. See also the docs for FindFit. Note also that you immediately obtain an excellent fit if you remove the restrictions on the values of your parameters. Most parameters still remain positive anyway, with the exception of W01, which becomes small and negative. Only you know whether that is acceptable... $\endgroup$
    – MarcoB
    Commented Apr 11, 2018 at 22:13
  • $\begingroup$ Hi, negative values aren't allowed for any of the parameters. The problem is even if I start with other starting values for the parameters, one of the sets of values {S, W0, Y} comes out to be very large compared to others, which doesn't at all make sense physically (according to the experimenter who gave me the data). $\endgroup$ Commented Apr 11, 2018 at 22:33
  • $\begingroup$ Is it possible that the model is not appropriate for the data then? Can you find a set of parameters that generates a shape similar to the data you have from the model you were given? $\endgroup$
    – MarcoB
    Commented Apr 11, 2018 at 22:51
  • $\begingroup$ Yes, in fact, if I remove the Y in the numerator and the factor (2[Pi]/377), it fits very smoothly(pasteboard.co/HgbCdfo.png). But the person who gave me the data is very clear that the model in question is the right one. So I've been dealing with it for a while now and without much success. $\endgroup$ Commented Apr 11, 2018 at 23:16
  • $\begingroup$ @NabinBhatta Your fit is quite good with less parameters. Adding Y's is not a good idea, since it is then you may over-fitting. Always prefer models as simple as possible, i.e. minimum # of parameters with minimum error. $\endgroup$ Commented Apr 12, 2018 at 3:59

2 Answers 2

4
$\begingroup$

Part of the problem is that there is a massive amount of data at low frequencies and a sparse amount of data at high frequencies making the fitting with linear frequency skewed.

There are 4379 data point. Below is a plot for the first 4329 (black) data points and the data points from 4330 to 4379.

ListPlot[{data[[1 ;; 4329]], data[[4330 ;; -1]]}, 
 PlotStyle -> {Black, Red}, PlotRange -> {{0, 50000}, {0, 8000}}]

Mathematica graphics

The shape of a function that has the form below

Mathematica graphics

is a peak whose amplitude is controlled by S4 and center by W4. The parameter Y4 also affects the amplitude but also controls the sharpness of the peak. See the Manipulate below.

Manipulate[

 Show[
  ListPlot[data, PlotStyle -> Black, 
   PlotRange -> {{0, 50000}, {0, 8000}}],
  Plot[(2 \[Pi]/377) (S4^2*Y4*w^2/((W04^2 - w^2)^2 + Y4^2*w^2)),
   {w, 20000, 50000}, PlotStyle -> Red, PlotRange -> Full]
  ],

 {{S4, 64455}, 5000, 80000, Appearance -> "Open"},
 {{Y4, 43240}, 500, 50000, Appearance -> "Open"},
 {{W04, 43240}, 0, 50000, Appearance -> "Open"}
 ]

Mathematica graphics

By playing with the Manipulate some reasonable starting parameters were determined.

High Frequency data

Parameters for the hump in the high frequency were determined in a separate fitting process.

model4 = (2 \[Pi]/377) (S4^2*Y4*w^2/((W04^2 - w^2)^2 + Y4^2*w^2))

nlm4 = NonlinearModelFit[data[[4330 ;; -1]],
         {model4, {S4 > 0, Y4 > 0, W04 > 0}}, {S4, Y4, W04}, w]

nlm4["BestFitParameters"]

{S4 -> 64456.4, Y4 -> 43241.5, W04 -> 42416.4}

Show[
 ListPlot[data[[4330 ;; -1]], PlotStyle -> Black, 
  PlotRange -> {{0, 50000}, {0, 8000}}],
 Plot[nlm4[w], {w, 20000, 50000}, PlotStyle -> Red, PlotRange -> Full]
 ]

Mathematica graphics

Full Model

I found it necessary to constrain the W0n parameters. In particular the late frequency data needed to be constrained to a high value or it would be swamped by the large amount of early frequency data (alternatively using a weighting function related to frequency might be the way to go).

I used the Manipulate above to get reasonable starting values for parameters in group 1 and 2 and allowed group 3 to be free.

 model = (2 π/377) (
 S1^2*Y1*w^2/((0 - w^2)^2 + Y1^2*w^2) + 
 S2^2*Y2*w^2/((W02^2 - w^2)^2 + Y2^2*w^2) + 
 S3^2*Y3*w^2/((W03^2 - w^2)^2 + Y3^2*w^2) + 
 S4^2*Y4*w^2/((W04^2 - w^2)^2 + Y4^2*w^2)
 )

nlm = NonlinearModelFit[
  data, 
  {
   model,
   {S1 > 0, Y1 > 0,
    S2 > 0, Y2 > 0, W02 > 0.5,
    S3 > 0, Y3 > 0, 10000 > W03 > 10,
    S4 > 0, Y4 > 0, W04 > 30000
    }
   },
  {
   {S1, 40000}, {Y1, 6000},
   {S2, 13000}, {Y2, 700}, {W02, 10},
   S3, Y3, W03,
   {S4, 65000}, {Y4, 50000}, {W04, 45000}
   },
  w
  ]

The parameters are approximately

nlm["BestFitParameters"]

{S1 -> 40942.1, Y1 -> 6071.72, S2 -> 13472.8, Y2 -> 759.179, 
 W02 -> 4.52962, S3 -> 354.288, Y3 -> 3.10132, W03 -> 518.807, 
 S4 -> 67601.3, Y4 -> 50274.8, W04 -> 46349.4}

Below is a comparison of the model and data

Show[
 ListPlot[data, PlotStyle -> Black, 
  PlotRange -> {{0, 50000}, {0, 8000}}],
 Plot[nlm[w], {w, 0, 50000}, PlotStyle -> Red, PlotRange -> Full]
 ]

Mathematica graphics

The residuals (difference between the fit and data) are shown below

Quiet@ListPlot[ Transpose[{data[[All, 1]], nlm["StandardizedResiduals"]}], PlotStyle -> Black, PlotRange -> Full]

Mathematica graphics

Weighted model

A plot of the data versus logarithm of the frequency

ListLogLinearPlot[data, PlotStyle -> Black, PlotRange -> Full]

Mathematica graphics

makes it appear that the low frequency data is poorer quality than the high frequency data.

Let's try the same fit weighting the data by the log of the frequency.

nlmW = NonlinearModelFit[
  data, 
  {
   model,
   {S1 > 0, Y1 > 0,
    S2 > 0, Y2 > 0, W02 > 0.5,
    S3 > 0, Y3 > 0, 10000 > W03 > 10,
    S4 > 0, Y4 > 0, W04 > 30000
    }
   },
  {
   {S1, 40000}, {Y1, 6000},
   {S2, 13000}, {Y2, 700}, {W02, 10},
   S3, Y3, W03,
   {S4, 65000}, {Y4, 50000}, {W04, 45000}
   },
  w,
  Weights -> Log[data[[All, 1]]]
  ]

The parameters change with the weighting

nlmW["BestFitParameters"]

{S1 -> 40014.8, Y1 -> 7761.67, S2 -> 15067.5, Y2 -> 929.066, 
 W02 -> 361.469, S3 -> 11880.2, Y3 -> 3343.04, W03 -> 2073.73, 
 S4 -> 52602.7, Y4 -> 29854.2, W04 -> 42160.2}

The higher frequency data now fits better with a corresponding degradation in the low frequency data

Show[
 ListPlot[data, PlotStyle -> Black, 
  PlotRange -> {{0, 50000}, {0, 8000}}],
 Plot[nlmW[w], {w, 0, 50000}, PlotStyle -> Red, PlotRange -> Full]
 ]

Mathematica graphics

Quiet@ListPlot[
  Transpose[{data[[All, 1]], nlmW["StandardizedResiduals"]}], 
  PlotStyle -> Black, PlotRange -> Full]

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ Thank you! Your answer helped a lot! Btw, is there any way that I can plot the 4 components of the model separately in the same frame so that it is easier to see what constituted the fitted model? $\endgroup$ Commented Apr 12, 2018 at 23:07
2
$\begingroup$

(This answer is not very detailed, I might elaborate it later on... I hope it helps ...)

The data is not evenly spaced across the x-axis so it is a good idea to sample it in a way that produces more uniform x-values distribution.

After this is done, we still do not get good fitting results with all data samples. That is why I made the search loop below. After running it with 100 iterations, a few times I got good fits, like this one:

enter image description here

Here is the function from one of the those good fits:

2/377 \[Pi] ((1.971*10^12 #1^2)/(1.24692*10^7 #1^2 + #1^4) + (
    1.19479*10^13 #1^2)/(7.89216*10^7 #1^2 + (61.0327 - #1^2)^2) + (
    6.63511*10^10 #1^2)/(415362. #1^2 + (135157. - #1^2)^2) + (
    7.35593*10^13 #1^2)/(
    7.93947*10^8 #1^2 + (1.75308*10^9 - #1^2)^2)) &

Let us check it against the data:

ordInds = Ordering[sRes[[All, 2]]]
nlmFit = sRes[[First[ordInds], 1]]
fitVals = {#, nlmFit["Function"][#]} & /@ data[[All, 1]];

ListPlot[Transpose[{data[[All, 
    1]], (fitVals[[All, 2]] - data[[All, 2]])/data[[All, 2]]}], 
 PlotRange -> All]

enter image description here

And here is the fit and the data together:

Show[{ListLinePlot[{#, nlmFit["Function"][#]} & /@ data[[All, 1]], 
   PlotRange -> All], 
  ListPlot[data, PlotStyle -> {PointSize[0.004], Red}]}]

enter image description here

Code

Clear[GetDataSample]
GetDataSample[data_?MatrixQ, nBreaks_Integer, nCellSamples_Integer] :=    
    Block[{qs, breaks, data2, t},
   qs = Quantile[data[[All, 1]], Range[0, 1, 0.05]];
   breaks = 
    Range[Min[data[[All, 1]]], 
     qs[[-2]], (qs[[-2]] - Min[data[[All, 1]]])/nBreaks];
   breaks = Append[breaks, qs[[-1]]];
   data2 = 
    Map[Function[{q}, t = Select[data, q[[1]] <= #[[1]] < q[[2]] &]; 
      If[q[[2]] == Max[data[[All, 1]]], t, 
       RandomSample[t, Min[nCellSamples, Length[t]]]]], 
     Partition[breaks, 2, 1]];
   data2 = SortBy[Join @@ data2, First]
   ];

sRes =
  Table[(
    Print[i];
    data2 = GetDataSample[data, 10, 10];
    nlmFit2 = 
     NonlinearModelFit[
      data2, {model, {S1 > 0, S2 > 0, S3 > 0, S4 > 0, W02 > 0, 
        W03 > 0, W04 > 0, Y1 >= 0, Y2 >= 0, Y3 >= 0, Y4 >= 0}}, {S1, 
       S2, S3, S4, W02, W03, W04, Y1, Y2, Y3, Y4}, w, 
      AccuracyGoal -> 5, MaxIterations -> 10000];
    fitVals = {#, nlmFit2["Function"][#]} & /@ data[[All, 1]];
    diffs = (fitVals[[All, 2]] - data[[All, 2]])/data[[All, 2]];
    t = {nlmFit2, Max[Abs[diffs]], 
      Quantile[diffs, Range[0, 1, 0.1]]};
    Print[t];
    Print[
     ListLogPlot[{data2, {#, nlmFit2["Function"][#]} & /@ 
        data2[[All, 1]]}, PlotRange -> All]];
    t), {i, 100}];
$\endgroup$
1
  • $\begingroup$ Thanks! It was really helpful. $\endgroup$ Commented Apr 12, 2018 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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