Fitting parametric Spiral Function to Data Points

Question

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]

Answer 1

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]

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.

Answer 2

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]

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]

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]