6
$\begingroup$

How would I fit the following parametric function describing a spiral to a set of data points?

spiral = t*{Cos[t + x], Sin[t + x]}

x is the unknown parameter which I want to get from the fit.

Since this is a geometric problem, the correct solution would minimize the sum of squares of distances between the data points and the curve.

The following code creates a possible set of data points:

parameterToBeFound = x -> 0;
noise = 0.1;
numOfPoints = 100;
tMin = 2 Pi;
tMax = 6 Pi;
dataPoints = Table[(spiral /. parameterToBeFound) + {RandomReal[t*noise] - 
  t*noise/2, RandomReal[t*noise] - t*noise/2}, {t, tMin, tMax, (tMax - tMin)/(numOfPoints - 1)}];
ListPlot[dataPoints, AspectRatio -> Automatic]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ I assume these are not your real data points. For the real ones, do you have the t values that created them, or just a set of x and y coordinates? $\endgroup$ – Michael Seifert May 8 '17 at 13:58
  • $\begingroup$ @Michael Seifert: You are correct; these are no real data points. They have been generated by the provided code. $\endgroup$ – Gerhard Wonner May 10 '17 at 4:37
  • $\begingroup$ @Michael Seifert: I just wanted to get a hint for my more complex real data Problem. $\endgroup$ – Gerhard Wonner May 10 '17 at 5:05
6
$\begingroup$

Essentially one needs to fitt two separate curves but with a common parameter. Also, one needs to recognize what sort of error structure might be involved. (It's not just about the regression parameters but also about the error structure.)

Suppose that the model generating the data is of the following form:

$$x = t \cos{(t+x)} + e_x$$ $$y = t \sin{(t+x)} + e_y$$

with $e_x$ and $e_y$ being independent with a common normal distribution having mean zero and standard deviation $\sigma$. The model that we fit will also be of this form.

First generate some data (not with a uniform distribution of errors):

spiral = t*{Cos[t + x], Sin[t + x]};
parameterToBeFound = x -> 0;
σ = 1;
numOfPoints = 100;
tMin = 2 Pi;
tMax = 6 Pi;
dataPoints = Table[{t, (spiral[[1]] /. parameterToBeFound) + 
  RandomVariate[NormalDistribution[0, σ], 1][[1]],
  (spiral[[2]] /. parameterToBeFound) + RandomVariate[NormalDistribution[0, σ], 1][[1]]},
  {t, tMin, tMax, (tMax - tMin)/(numOfPoints - 1)}];

Now we restructure the data so that we can use NonlinearModelFit. We create two dummy variables to indicate which curve is being fit.

data = Flatten[Table[{{dataPoints[[i, 1]], 1, 0, dataPoints[[i, 2]]},
  {dataPoints[[i, 1]], 0, 1, dataPoints[[i, 3]]}}, {i, numOfPoints}], 1];

Execute NonlinearModelFit:

nlm = NonlinearModelFit[data, c t Cos[t + x] + s t Sin[t + x], x, {t, c, s}];
Show[ListPlot[dataPoints[[All, {2, 3}]]],
  ParametricPlot[{nlm[t, 1, 0], nlm[t, 0, 1]}, {t, tMin, tMax}], AspectRatio -> 1]

Spiral fit and data

If the errors for the two curves have different variances and/or are correlated, then something else needs to be used and the above might only be able to get one starting values for that other procedure.

$\endgroup$
  • $\begingroup$ Pretty slick method; I like the use of the extra dummy parameters to simultaneously fit both the $x$ and $y$ data. That said, I'm not sure that it can be extended to work if you don't have the $t$ values for the data points. (See my question in the comments below the original question; hopefully the OP will return some time to clarify.) $\endgroup$ – Michael Seifert May 8 '17 at 17:45
  • $\begingroup$ Agreed: gotta have the $t$ values or at least the order of the data points and spacing of the $t$ values (which I guess is pretty much the same as knowing the $t$ values). $\endgroup$ – JimB May 8 '17 at 17:52
  • $\begingroup$ Unfortunately, in my real data problem I do not have the t-values. However, maybe I can guess them somehow. $\endgroup$ – Gerhard Wonner May 10 '17 at 5:10
  • $\begingroup$ One can probably perform an iterative method when the t-values are not available (although as a statistician I would have concerns about obtaining appropriate measures of precision for the estimator of x). One would guess at a value of x. Then for each data point choose the value of t on the curve that is closest to the data point. Then the sum of squares of distances between the data points and associate points on the curve is calculated. Find the value of x that minimizes the sum of squares. $\endgroup$ – JimB May 11 '17 at 16:25
  • $\begingroup$ That is exactly what I need. Since this is a geometric problem, the correct solution would minimize the sum of squares of distances between the data points and the curve. I added this information to my question. Also added geometry tags. $\endgroup$ – Gerhard Wonner May 12 '17 at 5:07
5
$\begingroup$

In what follows, I used $x = 1 + 1/e \approx 1.36787944...$, since $x = 0$ and $x = 1$ are often "special" values in curve-fitting routines and I didn't want to come up with a result that doesn't work for values other than these "special" ones.

We can transfer the data into polar coordinates. In terms of this, we have the parametric curve $r = t$ and $\theta = t + x$, or $\theta = r + x$. However, $\theta$ is only defined between $-\pi$ and $\pi$, so really we have $$ \theta(r) = (r + x) \mod 2\pi, $$ where we define the values of $\text{mod}$ to lie between $-\pi$ and $\pi$. The code to do this transfer is:

polardata = CoordinateTransform["Cartesian" -> "Polar", #] & /@ dataPoints;
ListPlot[polardata]

enter image description here

I attempted to use the pre-fab function NonLinearModelFit to find the value of $x$ from this graph; but it tended to have trouble, presumably because of the discontinuity in the $\text{mod}$ function. I ended up defining my own least-squares fitting function, modulo $2 \pi$. The squared error in a given data point $(r_i, \theta_i)$ will be $$ \{ [\theta_i - (r_i + x)] \, \text{mod} \, 2 \pi \}^2, $$ and the total squared error will be the sum of this over all data points. We want to find the value of $x$ that minimizes this total squared error.

pointsquareerror[point_, x_] := Mod[(point[[2]] - (point[[1]] + x)), 2 π, -π]^2
totalsquareerror[x_] := Total[pointsquareerror[#, x] & /@ polardata];
Plot[totalsquareerror[x], {x, 0, 2 π}, PlotPoints -> 200]

enter image description here

We can see that there is a minimum at approximately $x = 1.4$, and using FindMinimum we can get it a little more accurately:

fitx = FindMinimum[totalsquareerror[x], {x, 0, 2 π}]
(* {11.1933, {x -> 1.41496}} *)

Comparing:

dataplot = ListPlot[dataPoints, AspectRatio -> Automatic];
fitplot = ParametricPlot[t*{Cos[t + x], Sin[t + x]} /. minx[[2]], {t, 2 π, 6 π}, PlotStyle -> Red];
Show[dataplot, fitplot]

enter image description here

$\endgroup$
  • $\begingroup$ If anyone can get NonLinearModelFit to work reliably, I'd be interested to see it. In my experience, sometimes it worked fine, and sometimes it failed utterly (even when given an initial "guess" for $x$.) $\endgroup$ – Michael Seifert May 8 '17 at 15:18
  • $\begingroup$ I hope my answer shows how to use NonlinearModelFit for this data and model. $\endgroup$ – JimB May 8 '17 at 17:40

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.