Least squares non-linear estimate in R

lsnl.R

library("R.matlab")
library("tidyverse")
library("pracma")
library("ggforce")

coord_data <- readMat("cirkel300.mat")

tbl <- as_tibble(list(X=coord_data$X[,1], Y=coord_data$Y[,1]))
# center the data
# tbl$X <- tbl$X - mean(tbl$X)
# tbl$Y <- tbl$Y - mean(tbl$Y)
# find the least squares best fit
m1 <- lsqnonlin(function(p) {
  with(tbl, sqrt((X - p[1])^2 + (Y - p[2])^2) - p[3])
}, c(0,0,1))


max_x <- max(tbl$X)
min_x <- min(tbl$X)
max_y <- max(tbl$Y)
min_y <- min(tbl$Y)

# plot the data and add the circle
plt <- tbl %>% 
  ggplot(aes(x=X, y=Y)) +
  geom_point(aes(color="points")) +
  coord_fixed(ratio = 1) +
  geom_circle(aes(x0=m1$x[1], y0=m1$x[2], r=m1$x[3], color="lsq"), size=1) +
  geom_point(aes(x=m1$x[1], y=m1$x[2], color="centre"), shape=3, size=3) +
  scale_color_manual(
    breaks = c("points", "lsq", "centre"),
    values= c("points"="red", "lsq"="blue", "centre"="blue"),
    labels = c("Points", "Least squares", "Centre")
    ) +
  guides(
    color = guide_legend(
        order = 1,
        override.aes = list(
            shape = c(
                16,
                NA,
                3),
            linetype = c(
                "blank",
                "solid",
                "blank"
            )
        )
    )
  ) +
  xlim(min_x, max_x) +
  ylim(min_y, max_y) + 
  labs(x="X", y="Y", title="Circle fit") +
  theme_minimal()

# add center coordinate label to plot
plt <- plt + annotate("text", x=m1$x[1], y=m1$x[2] - 0.25, label=paste0("(", round(m1$x[1], 2), ", ", round(m1$x[2], 2), ")"), size=3)

# add radius label to plot
plt <- plt + annotate("text", x=max_x - 1, y=max_y - 1, label=paste0("Radius: ", round(m1$x[3], 2)), size=3)

print(plt)


# calculate the euclidean distance between the center and each point
tbl$dist <- sqrt((tbl$X - m1$x[1])^2 + (tbl$Y - m1$x[2])^2)
tbl$residual <- tbl$dist - m1$x[3]

# plot the distance from the center
tbl %>% 
  ggplot(aes(x=seq(1, length(dist)), y=dist)) +
  geom_point() +
  stat_smooth(method="lm", formula=y~x, se=FALSE) +
  geom_segment(aes(xend=seq(1, length(dist)), yend=m1$x[3], color="residuals"), size=0.5, alpha=0.5) +
  geom_hline(yintercept=mean(tbl$dist), color="red") +
  labs(x="Point", y="Distance", title="Euclidean Distance from Center") +
  annotate("text", x=1, y=mean(tbl$dist) + 0.01, label=paste0("Mean: ", round(mean(tbl$dist), 2)), size=3) +
  theme_minimal()

r$> m1
$x
[1] -5.953157  3.032570  5.025793

$ssq
[1] 23.99621

$ng
[1] 2.727779e-06

$nh
[1] 3.770407e-08

$mu
[1] 0.03864297

$neval
[1] 11

$errno
[1] 1

$errmess
[1] "Stopped by small x-step."