Findfit doesn't give the good fit; Changing the starting values will not change the results

Question

I have some complex data:

data={{0.45, 0. + 0.0119333 I}, {0.46, 0. + 0.0100877 I}, {0.47, 0. + 0.00584566 I}, {0.48, 0.00829149 + 0.0000687535 I}, {0.49, 0.0146198 + 0.000213785 I}, {0.5, 0.0202219 + 0.000409093 I}, {0.51,0.0257641 + 0.000664232 I}, {0.52, 0.0314458 + 0.000989817 I}, {0.53, 0.0373584 + 0.0013976 I}, {0.54, 0.0435542 + 0.00190058 I}, {0.55, 0.0500677 + 0.00251309 I}, {0.56, 0.0569241 + 0.00325092 I}, {0.57, 0.0641438 + 0.0041315 I}, {0.58,0.0717442 + 0.005174 I}, {0.59, 0.0797407 + 0.00639954 I}, {0.6, 0.0881481 + 0.00783143 I}, {0.61, 0.0969803 + 0.00949535 I}, {0.62, 0.106251 + 0.0114197 I}, {0.63, 0.115974 + 0.0136358 I}, {0.64, 0.126162 + 0.0161785 I}, {0.65, 0.136829 + 0.0190864 I}, {0.66, 0.147988 + 0.0224023 I}, {0.67, 0.159653 + 0.0261742 I}, {0.68, 0.171837 + 0.0304556 I}, {0.69, 0.184553 + 0.0353065 I}, {0.7, 0.197814 + 0.0407944 I}, {0.71, 0.21163 + 0.0469957 I}, {0.72, 0.226011 + 0.0539967 I}, {0.73, 0.240966 + 0.0618958 I}, {0.74, 0.2565 + 0.0708055 I}, {0.75, 0.272612 + 0.0808548 I}, {0.76, 0.289298 + 0.0921927 I}, {0.77, 0.306543 + 0.104992 I}, {0.78, 0.324321 + 0.119453 I}, {0.79, 0.342588 + 0.135812 I}, {0.8, 0.361278 + 0.154344 I}, {0.81, 0.380289 + 0.175377 I}, {0.82, 0.39947 + 0.199295 I}, {0.83, 0.418604 + 0.226557 I}, {0.84, 0.43737 + 0.257703 I}, {0.85, 0.455307 + 0.293371 I}, {0.86, 0.471745 + 0.3343 I}, {0.87, 0.485713 + 0.381328 I}, {0.88, 0.495804 + 0.435363 I}, {0.89, 0.499993 + 0.497286 I}, {0.9, 0.495388 + 0.567757 I}, {0.91, 0.477963 + 0.646802 I}, {0.92, 0.442359 + 0.733064 I}, {0.93, 0.382045 + 0.822555 I}, {0.94, 0.290457 + 0.906983 I}, {0.95, 0.163934 + 0.972362 I}, {0.96,0.00673389 + 0.999955 I}, {0.97, -0.16425 + 0.972252 I}, {0.98, -0.320906 + 0.883431 I}, {0.99, -0.435064 + 0.746414 I}, {1., -0.492227 + 0.587822 I}}

which want to fit the following function to them:

function=b + a/((-0.977727 + 0.0601085 I) + s);

Without Using any starting values:

ans = FindFit[data, function, {a, b}, s]

{a -> -0.0650086 - 0.0165728 I, b -> -0.0131795 - 0.0513961 I}

The results for the real part of the fitted function is:

 Show[Plot[Re[Evaluate[function /. ans]], {s, 0.45, 1.008}, 
  PlotStyle -> Directive[Red, Thickness[0.003]], 
  PlotRange -> {{0.45, 1.1}, {-0.5, 0.53}}, Frame -> True, 
  BaseStyle -> FontSize -> 16], 
 ListPlot[Re[data], PlotStyle -> Black]]

which is not good enough. Using Manipulate for choosing good starting values doesn't change the result:

func[s1_, a1_, b1_] := function /. s -> s1 /. a -> a1 /. b -> b1

Manipulate[
 Show[ListPlot[Re[data]], 
  Plot[Re[func[s1, a1, b1]], {s1, 0.45, 1.008}, PlotStyle -> Red], 
  Frame -> True, Axes -> None], {a1, -0.2 - 0.1 I, 
  0.3 + 0.2 I}, {b1, -0.2 + 0.2 I, 0.2 - 0.2 I}]


 FindFit[data, 
  function, {{a, -0.10951986312866213` - 
     0.04571191787719727` I}, {b, -0.10639984130859376` + 
     0.10639984130859376` I}}, s]

{a -> -0.0650086 - 0.0165728 I, b -> -0.0131795 - 0.0513961 I}

How can I find a better fit?

Answer 1

My answer below is mostly inspired by the answers of @Matariki and @george2079, so thank you both for your efforts and insight. I also got some ideas (about minimizing least square) from here (warning: it's a Matlab SO).

My approach is to minimize the least square errors/residuals/differences between the data and the fit function in terms of both real and imaginary. For example, for the first data pairs {0.45, 0. + 0.0119333 I}, s = 0.45, real y-value = 0, and imaginary y-value = 0.0119333. With those 3 numbers, I could construct 2 equations, by plugging in s on the RHS of your function and either the real part or the imaginary part on the LHS. I then take the LHS i.e. data minus RHS i.e. model, like so:

(* 1: Real data - real model *)
Re[data[[1, 2]]] - Re[function /. s -> data[[1, 1]]]
(* 0. - Re[(-1.87065 - 0.213068 I) a + b] *)

(* 2: Imaginary data - imaginary model *)
Im[data[[1, 2]]] - Im[function /. s -> data[[1, 1]]]
(* 0.0119333 - Im[(-1.87065 - 0.213068 I) a + b] *)

I would do this for all 56 data points provided to get 112 equations, 56 real and 56 imaginary. Then, I would use NMinimize to minimize the least squares of these LHS - RHS differences. The goal is to make the least square residuals of these differences to be as close to 0 as possible.

equationsForFunction = 
  Join[Re@data[[All, 2]] - Re[function /. s -> data[[All, 1]]] (* all 56 real equations *), 
   Im@data[[All, 2]] - Im[function /. s -> data[[All, 1]]]] (* all 56 imaginary equations *);
ansForFunction = NMinimize[#.# &@equationsForFunction, {a, b}] (* minimize the least square residual *)
(* {0.918251, {a -> -0.0658589, b -> -0.0721499}} *)

As you could see the final residual (0.918251) is pretty small, and a and b were found. This was indeed @Oleksandr R. answer in his comment.

However, a and b need not be real--they could be complex--and the complex number in the denominator of your function---0.977727 + 0.0601085 I--needs not have those values for real and imaginary parts; as Matariki pointed out, that number is the reason why you got the poor fit no matter how much you optimize your a and b.

As a result, I constructed 2 more fitting functions, one with a and b as complex numbers (Rea + Ima I and Reb + Imb I), and the other with complex a & b and the complex number in the denominator as some yet-to-be-determined c (or more accurately, Rec + Imc I) to see if they would fit the data better. Below are these 3 functions:

function = b + a/((-0.977727 + 0.0601085 I) + s); (* real a, b *) 
functionComplex = 
  Reb + I Imb + (Rea + I Ima)/((-0.977727 + 0.0601085 I) + s); (* complex a, b *)
functionBetter = Reb + I Imb + (Rea + I Ima)/((Rec + I Imc) + s) (* complex a, b, c *);

And here is my use of NMinimize to find the parameters of these equations:

{equations, equationsComplex, equationsBetter} =
  Join[Re@data[[All, 2]] - Re[# /. s -> data[[All, 1]]], 
     Im@data[[All, 2]] - Im[# /. s -> data[[All, 1]]]] & /@ {function,
     functionComplex, functionBetter};
ans = NMinimize[#.# &@#1, #2] & @@@ {{equations, {a, 
      b}}, {equationsComplex, {Rea, Ima, Reb, 
      Imb}}, {equationsBetter, {Rea, Ima, Reb, Imb, Rec, Imc}}};
Grid[Join[{{"Least square residuals", "Parameters"}}, ans], 
 Frame -> All]

As you can see the equation with complex a, b, and c gives the best fit because it has the least square residual out of the 3 equations. Notice that my solution found from functionComplex (i.e. complex a & b) matches up with @george2079's answer using Norm to minimize both real and imaginary differences between fit and data.

Now let's plot the real part of these equations and compare them with the original data (I also plot Matariki's equation for comparison):

Show[ListPlot[Re[data], PlotStyle -> Gray],
 Plot[{Re[function /. ans[[1, 2]]], 
   Re[functionComplex /. ans[[2, 2]]], 
   Re[functionBetter /. ans[[3, 2]]], Re[functionMatariki]}, {s, 0, 
   1}, PlotRange -> Full, 
  PlotLegends -> {"function: real a, b", 
    "functionComplex: complex a, b", 
    "functionBetter: complex a, b, c", "functionMatariki"}]]

And plotting the imaginary parts gives:

Show[ListPlot[Transpose[{data[[All, 1]], Im[data][[All, 2]]}], 
  PlotStyle -> Gray], 
 Plot[{Im[function /. ans[[1, 2]]], 
   Im[functionComplex /. ans[[2, 2]]], 
   Im[functionBetter /. ans[[3, 2]]], Im[functionMatariki]}, {s, 0, 
   1}, PlotRange -> Full, 
  PlotLegends -> {"function: real a, b", 
    "functionComplex: complex a, b", 
    "functionBetter: complex a, b, c", "functionMatariki"}]]

As you can see from the plots, the function with complex a, b, and c does look like it fits the data the best.

Answer 2

Your model function isn't allowing for a better fit. Are you sure that the model is correct? Below is a plot using a different model which allows to find better start values for a fit. Mind you I haven't done a fit.

Show[ListPlot[Re[data]], 
 Plot[Re[b + a/((-1.028 + 0.0492 I) + s)] /. {a -> -0.248 - 0.306 I, 
    b -> -0.403 + 0.106 I}, {s, 0.45, 1.008}, PlotStyle -> Red], 
 Frame -> True, Axes -> None] 

Answer 3

If you fit the real part of the data directly:

  FindFit[ Re[data] , 
      Re[b + a/((-0.977727 + 0.0601085 I) + s)] , {{a, -.1}, {b, -.1}}, s]

{a -> -0.0791299, b -> -0.121803}

Gives a pretty good fit,

and surprisingly not too bad in the complex part either.

Better ... <-- not..

This explicitly does the least squares minimization over the complex values:

   NMinimize[ 
       Total[ Norm[(b + ib I ) +
             (a + ia I)/((-0.977727 + 0.0601085 I) + #[[1]]) - #[[2]]]^2 & /@
                  data] , {a, ia, b, ib}]

{a -> -0.0650086, ia -> -0.0165728, b -> -0.0131795, ib -> -0.0513963}

the resulting fits are pretty similar to above (ok, not great) so i wont bother posting..

edit.. of course @Oleksandr is right this is precisely the same result as obtained by FindFit.

In any case you may find this interesting:

 Manipulate[
    f[s_] = br + bi I  + (ar + I ai) /(-0.977727 + 0.0601085 I + s) ;
    Show[{
     Graphics[Line[{Re[#[[2]]], Im[#[[2]]]} & /@ data]],
     Graphics[
        Line[{{Re[#[[2]]], Im[#[[2]]]}, {Re[#], Im[#]} &@f[#[[1]]]}] & /@ 
         data],
    ParametricPlot[{Re[#], Im[#]} &@f[s], {s, .5, 1}]}],
      {{br, .1}, -.1, .1}, {{bi, .0015}, -.1, .1},
      {{ar,-.059}, -.1, .1}, {{ai, -.0125}, -.1, .1}]

Your functional form it turns out can fit your path in the complex plane very well, if the parameter 's' doesn't need to align with the arc length parameter of the data.