3
$\begingroup$

I have a list of numbers that, taken in order, describe a spiral in the complex plane. Here is a small subset of it, sufficient to illustrate the problem:

list = {0.01263537001717932 + 0.07067349025848388*I, 
 0.047250267751479655 + 0.05209549861304419*I, 
 0.06596386316698546 + 0.02009477948782494*I, 
 0.06598609040888194 - 0.014961756210224754*I, 
 0.049803550815569794 - 0.04397149975599973*I, 
 0.023258224576407013 - 0.06099441316314905*I, 
 -0.006595117201610087 - 0.06383941766600312*I,
 -0.03328082763577616 - 0.053650337718323365*I,
 -0.05206972138609968 - 0.03391359218040549*I, 
 -0.060484300582113376 - 0.009278027786894835*I, 
 -0.05825268030262476 + 0.015525647453509116*I, 
 -0.04688894972912606 + 0.03649893985640911*I, 
 -0.029084475019386297 + 0.05086222796961101*I, 
 -0.008065197460458237 + 0.05723002594987413*I, 
 0.012975410250310286 + 0.05553656470450358*I, 
 0.031305061954391124 + 0.046792399011505864*I, 
 0.04492053689497534 + 0.03275238482001683*I, 
 0.05266179241274281 + 0.015563156240077383*I, 
 0.054197917274563315 - 0.002559520264756074*I, 
 0.04991942530288897 - 0.019591820266882193*I, 
 0.04077290783116659 - 0.03388637260674089*I, 
 0.028070627915692578 - 0.044278439184983934*I, 
 0.01330140973390842 - 0.05012128281294888*I, 
 -0.002038137442440774 - 0.051263817289083816*I, 
 -0.0165802068618294 - 0.04798609535183286*I, 
 -0.02918463847368017 - 0.04090829084508081*I, 
 -0.03899948965681603 - 0.030887325316889896*I, 
 -0.0454864418481062 - 0.01891284287902329*I, 
 -0.04841647916489085 - 0.006011387028409941*I, 
 -0.04784247104545147 + 0.0068352246318077806*I, 
 -0.044055588207469146 + 0.018753932199973615*I, 
 -0.03753212093089012 + 0.029025665281149138*I, 
 -0.028876474184390682 + 0.03711227916123358*I, 
 -0.018765076191573443 + 0.04266692897277542*I, 
 -0.007894807565850093 + 0.04553040163855704*I, 
 0.0030615501673451687 + 0.04571633159705858*I, 
 0.013491358125476596 + 0.043388326822319356*I, 
 0.022869033045436377 + 0.038831902057961665*I, 
 0.030773906534906574 + 0.032423827007817124*I, 
 0.0368996370049279 + 0.02460111005525975*I, 
 0.041056149623619684 + 0.015831403761543628*I, 
 0.04316530554941611 + 0.006586176629337748*I, 
 0.043251606732296555 - 0.0026824244302409914*I, 
 0.04142925049656888 - 0.011560475046661134*I, 
 0.03788678231239426 - 0.01968671709393703*I, 
 0.03287047720230047 - 0.026762590899545816*I, 
 0.026667429059332358 - 0.0325581087282407*I, 
 0.019589158683998575 - 0.03691397289036705*I, 
 0.011956378358853229 - 0.03974044128341975*I, 
 0.004085383055940278 - 0.04101349823952391*I, 
 -0.0037236170298382623 - 0.04076890659052196*I}; 
 
ListLinePlot[ReIm[list], AspectRatio -> 1]

enter image description here

The full data set is thousands of entries long, and produces the following ListLinePlot:

enter image description here

I'd like to use FindFit to generate a model of the form a*x^b*E^(c*x) (assuming that is the most appropriate method and model). But the problem is that

Table[Arg[list[[j]]], {j, Length[list]}]

{1.393880369, 0.8341309359, 0.2957019528, -0.222970991, -0.723286234,
-1.206498742, -1.67373906, -2.126029227, -2.564296235, -2.989383456,
2.881124903, 2.480154188, 2.090243925, 1.710800332, 1.341276188,
0.9811660967, 0.6300023216, 0.2873511408, -0.0471903747,
-0.3739972084, -0.6934189241, -1.005781912, -1.311391395,
-1.610533215, -1.903475446, -2.190469834, -2.471753094, -2.747548088,
-3.018064883, 2.999683597, 2.739139467, 2.483310929, 2.232030327,
1.985138787, 1.742485618, 1.503927758, 1.269329268, 1.038560868,
0.8114995072, 0.5880279702, 0.3680345156, 0.1514125389,
-0.06193973712, -0.2721195498, -0.4792198625, -0.6833296112,
-0.8845339342, -1.082914385, -1.278549132, -1.47151314, -1.661878347}

produces arguments that do not respect the 'winding number' of the spiral. (Why it produces some negative arguments is also baffling to me, but secondary). Since the values of Arg are restricted to an interval of 2*Pi, and are not in regular increments, I can't apply FindFit to more than a single circuit of the origin.

I could of course use FindFit[list, a*x^b*Cos[c*x] - I*a*x^b*Sin[c*x], {a, b, c}, x] as the model, but this produces an unevaluated result - and even if it did produce a result, it would be harder to work with.

How do I go about fitting list to a meaningful model?

$\endgroup$
3
  • $\begingroup$ Should we assume that the sample points are equally spaced in terms of x? If so, then modelling the real and imaginary components separately might be the first step to find the form of the model. Maybe first taking a look at dataRe = Transpose[{Range[Length[list]], Re[list]}]; dataIm = Transpose[{Range[Length[list]], Im[list]}]; ListLinePlot[{dataRe, dataIm}]. (And a link to the complete dataset might be helpful.) $\endgroup$ – JimB Apr 23 at 16:30
  • $\begingroup$ Can we assume that each turn of the spiral has at least one data point? $\endgroup$ – A.G. Apr 23 at 17:51
  • $\begingroup$ Do you have a fit criteria in mind? (least squares, max distance, etc.) $\endgroup$ – A.G. Apr 23 at 17:52
9
$\begingroup$

It will be easier to fit the data if it does not wind around the origin. So, let's unwind the spiral. Applying the Arg function to the data list gives

ListLinePlot[Arg/@list]

enter image description here

We see that arguments are decreasing, except where they jump up by $2\pi$. We want to unwind the arguments by removing the jumps. Here is a function that will unwind the arguments

ClearAll[unwind]
unwind[arg_List] := 
 Block[{n, xlist = Thread[{Most[arg], Rest[arg]}]},
  n = First@
    FirstPosition[xlist, SelectFirst[xlist, Last[#] > First[#] &]];
  If[NumericQ[n],
   Join[arg[[;; n]], unwind[-2 \[Pi] + arg[[n + 1 ;;]]]],
   arg]
  ]
ListLinePlot[unwind[Arg /@ list], ImageSize -> Small]

enter image description here

Now we do the fit and plot the results

arg = unwind[Arg /@ list];
abs = Abs /@ list;
fit = FindFit[Thread[{arg, abs}], b x + c, {b, c}, x]

(* {b -> 0.00192386, c -> 0.0671975}  *)

ParametricPlot[ReIm[(b x + c /. fit) Exp[I x]], {x, -\[Pi], 5 \[Pi]},
 Epilog -> {PointSize[.02], Red, Point[ReIm[list]]},
 ImageSize -> Small]

enter image description here

This model uses a linear $b x + c$ relation to relate the radial distance to the angular coordinate. Other relationships may be possible.

Note that the unwind function only works with decreasing data that has $2\pi$ jumps. If we Reverse the data list we would have increasing data with $-2\pi$ jumps. unwind requires 2 modifications to handle increasing data (instead of decreasing data). First, change the > comparison to < and second, change $-2\pi$ to $2\pi$.

$\endgroup$
2
  • $\begingroup$ Very nice. If one uses a quadratic rather than just a linear model, the fit looks perfect to the eye. But using $AIC_c$ as a guide, a 14-degree polynomial is even better and supported statistically. $\endgroup$ – JimB Apr 24 at 4:47
  • $\begingroup$ Thanks to all of you. @JimB, any chance you could share that model? $\endgroup$ – Richard Burke-Ward Apr 24 at 11:15
2
$\begingroup$

@LouisB supplied the answer and this is just an extended comment illustrating how other models using that answer might be selected.

First, rather than a linear model with a polynomial of order 1, consider a quadratic (2nd degree polynomial).

arg = unwind[Arg /@ list];
abs = Abs /@ list;
lm2 = LinearModelFit[Transpose[{arg, abs}], {x, x^2}, x];

ParametricPlot[ReIm[lm2[x] Exp[I x]], {x, -5 π, π}, 
  Epilog -> {PointSize[.02], Red, Point[ReIm[list]]}, ImageSize -> Medium]

Data and 2nd degree polynomial fit

An absolute measure of the fit is the root mean square error for predicting abs from arg:

lm2["EstimatedVariance"]^0.5
(* 0.0000451409 *)

If by this measure the fit is good enough, one can stop. But it can pay to look at the residuals to see if there are any additional patterns that might be better fit with additional terms in the model.

ListPlot[Transpose[{lm2["PredictedResponse"], lm2["FitResiduals"]}]]

Predicted response vs fit residuals

Clearly there is more that could be explained possibly by adding in a 3rd order term.

lm3 = LinearModelFit[Transpose[{arg, abs}], {x, x^2, x^3}, x];
lm3["EstimatedVariance"]^0.5
(* 3.69339*10^-6 *)

This seems a much better fit.

This process of adding polynomial terms can continue until the $AIC_c$ is minimized. For this particular dataset, a 14-degree polynomial has the miniminum $AIC_c$ value. This is unlikely to hold for the larger dataset, however. (I mostly deal with biological data and I have never ever seen anything more that a 6-th degree polynomial that might be appropriate.)

aicc[p_] := LinearModelFit[Transpose[{arg, abs}], Table[x^i, {i, p}], x]["AICc"]
ListPlot[Table[{p, aicc[p]}, {p, 20}], Frame -> True, FrameLabel -> {"Degree of polynomial", "AICc"}]

AICc vs degree of polynomial

$\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.