These data are interesting in that they hide a subtle interdependence of errors. It can be seen only after fitting the model correctly. Choosing the right model makes a large difference in the estimates.
Let each observation be a point $Z_i$ in the plane, written with a capital letter to remind us it will be modeled as a random variable. Specifically, let $\zeta$ be the circle's center, $\rho$ its radius, and let $\theta$ be a reference angle. Then
$$Z_i = \zeta + \rho(\cos(\delta i + \theta), \sin(\delta i+\theta)) + E_i$$
describes points $Z_i$ on the circle that have been displaced randomly by vectors $E_i$ ($i=1,2,\ldots, 21$ for these data). The angular increment $\delta$ is one degree for these data. The data themselves are considered realizations of the $Z_i$, written $z_i$.
A least-squares solution is appropriate when the $E_i$ are independent and each has a circularly symmetric distribution around $(0,0)$ (generalizing the assumption of approximate normality). It finds values of $\zeta = (\xi,\eta)$, $\rho$, and $\theta$ minimizing the sum of squared residuals
$$f(\zeta,\rho,\theta) = \sum_i |z_i - (\zeta + \rho(\cos(\delta i + \theta), \sin(\delta i+\theta)))|^2.$$
Although the gradient of $f$ is a nonlinear function of $\theta$, $f$ is readily minimized provided we can supply a reasonable starting estimate for the parameters. One such estimate can be obtained by selecting three widespread points that typify the data locations. I have taken the centroids of three groups of five points: the first five, the last five, and the middle five. The starting parameters describe the circumcircle of these three points. Let's call this the "geometric fit."
The figure displays the geometric fit and the subsequent least squares fit obtained by starting with the geometric fit and applying a nonlinear solver to improve the value of $f$.
Visually, the geometric fit is excellent: the circular arc seems to split the data points (shown as hollow dots) pretty evenly. The three centroids used to find this fit are shown as solid dots.
The least squares fit does not look quite as good. To understand why not, I have drawn it a little differently. Recall that the model prescribes a location for each point along known angles. Thus, to each of those (21) known angles corresponds a definite predicted location. These are shown as small solid red dots. The correspondence between each observation $z_i$ and its fit $\hat z_i$ is indicated by drawing a light arrow from one towards the other.
For the first time we begin to see the role played by the sequence of points. They aren't just some geometric jumble, as shown in the left hand plot, but really are a two-dimensional time series with an associated structure. And an interesting structure it is! Look at how the gray arrows seem to spin around as they track along the fitted curve. The residuals are not independent. In other words, the locational error is changing systematically over time, and not entirely randomly as assumed.
This is evident in a plot of the residuals:
The x-coordinates are shown with circles (and a Loess smooth in black) while the y-coordinates are shown with triangles (and a Loess smooth in red). There is a systematic variation in the x-coordinates: they start high, drop low, then go back high. (Geometrically, the points start to the right of the fitted circle, move to its left, and then back to its right. They wobble up and down, too, but those wobbles look random.)
Because of this, the least-squares fit is not entirely trustworthy. However, because the x-residuals have at least cycled through one phase of high-low-high, the least-squares fit is probably pretty good compared to the geometric fit (which pays no attention to the underlying angle associated with each point and is confused by the systematic variation in the x-direction). Here are the two sets of parameter estimates:
zeta(1) zeta(2) rho theta
Geometric 0.3683045 0.1630800 0.3936789 -3.145837
Least squares 1.1618824 0.5750505 1.2705868 -2.860703
The least-squares radius of 1.27 is three times the geometric radius of 0.40. That large difference attests to the importance of using the right model to make the fit.
For an improved fit, look into multidimensional time-series methods. A simple one would be to adopt a distributional assumption and correlation structure for the residuals and perform the fitting using maximum likelihood.
This R
code did the work.
#
# The point sequence.
# All arrays of points will be maintained as 2 x n matrices, one point
# per column.
#
z <- matrix(scan(text="-0.04 0.1964
-0.0309 0.1764
-0.0099 0.1484
-0.0226 0.135
-0.0233 0.1332
-0.0261 0.1125
0.002 0.0633
0.0048 0.0404
-0.0011 0.0474
-0.0149 0.0353
-0.0073 -0.0017
0.0253 -0.019
0.0317 -0.0369
0.065 -0.069
0.0652 -0.1065
0.098 -0.1035
0.0837 -0.1064
0.1072 -0.1292
0.1179 -0.1717
0.1839 -0.175
0.1691 -0.178"), nrow=2, dimnames=list(c("x", "y")))
#
# The nominal angles, in increments of one degree.
#
theta <- 1/180 * pi * (1:dim(z)[2])
#
# The objective function.
#
f <- function(parms) {
zeta <- parms[1:2] # Center
rho <- parms[3] # Radius
theta.0 <- parms[4] # Start angle, radians
theta.d <- theta + theta.0
w <- z - (zeta + rho*rbind(cos(theta.d), sin(theta.d)))
sum(w^2)
}
#
# Rough estimate of starting values: use a triangle to approximate the data.
#
n <- dim(z)[2] # Number of data points
k <- 2 # Will use clusters of 2*k+1 points
pts <- sapply(c(k+1, floor(n/2), n-k),
function(i) c(mean(z[1, (i-k):(i+k)]), mean(z[2, (i-k):(i+k)])))
#
# Find the center and radius of the circumcircle.
#
v <- apply(pts, 1, diff) # Direction vectors for the triangle's sides
x <- solve(cbind(c(-v[1,2], v[1,1]), c(v[2,2], -v[2,1])), apply(v, 2, mean))
zeta <- v[1, ]/2 - x[1] * c(v[1,2], -v[1,1]) + pts[, 1]
rho <- sqrt(sum((p <- pts[, 1]-zeta)^2))
theta.0 <- atan2(p[2], p[1]) - theta[1] # Guess at the angle of origin
#
# Plot this fit.
#
par(mfrow=c(1,2))
a <- seq(0, 360, 1) / 180 * pi # Angular coordinates on a circle
plot(t(z), asp=1, main="Geometric fit") # The points
lines(t(rho*rbind(cos(a), sin(a)) + zeta))# The fitted circle (clipped)
points(t(pts), pch=16) # The triangle's vertices
points(matrix(zeta, 1, 2), pch=16) # The circle's center (clipped)
#
# Refine this estimate.
#
fit <- nlm(f, c(zeta, rho, theta.0)) # Newton-Raphson improvement
rbind(c(zeta, rho, theta.0), fit$estimate)# Compare the two fits
prediction <- fit$estimate[3]*rbind(cos(theta+fit$estimate[4]),
sin(theta+fit$estimate[4])) + fit$estimate[1:2]
#
# Plot it.
#
plot(rbind(t(z), t(prediction)), type="n", asp=1,
main="Least squares fit")
invisible(sapply(1:n, function(i)
arrows(z[1,i], z[2,i], prediction[1,i], prediction[2,i],
length=0.1, angle=20, col="#a0a0a0")))
lines(t(prediction), col="Red", lty=3)
points(t(z))
points(t(prediction), col="Red", pch=16, cex=0.7)
#
# Plot the residuals.
#
par(mfrow=c(1,1))
residuals <- z - prediction
k <- max(abs(residuals[2,]))
plot(c(0,n+1)*k, range(residuals[2,]), type="n",
asp = diff(range(residuals[2,])) / (k*1.1),
xaxt="n", bty="n",
xlab="", ylab="", main="Residual Vectors")
abline(h=0, col="Gray")
invisible(sapply(1:n, function(i)
arrows(i*k, 0, i*k+residuals[1,i], residuals[2,i], length=0.1, angle=20)))
#
# Plot the residual components.
#
plot(residuals[1,], type="n", ylab="", main="Residuals")
lines(lowess(residuals[2,], f=0.4), col="Red", lwd=2)
lines(lowess(residuals[1,], f=0.4), lwd=2)
points(residuals[1,])
points(residuals[2,], col="Red", pch=2)