R/ransac_nls.R

In RANSAC: Robust Model Fitting Using the RANSAC Algorithm

#' Robust Nonlinear Model Fitting via RANSAC
#'
#' Fits a robust nonlinear model (`nls`) using the RANSAC algorithm.
#'
#' @param formula Model formula.
#' @param data Data frame containing the model variables.
#' @param start Named list of initial parameter estimates (as required by `nls`).
#' @param n_min Minimum number of points to fit the model.
#' @param n_iter Number of iterations (higher values make the model more robust).
#' @param tol Absolute tolerance to consider a point as an inlier.
#' @param verbose If `TRUE`, shows progress messages.
#'
#' @return An `nls` model fitted only to the inliers, with an additional class `"ransac_nls"` and an `"inliers"` attribute.
#'
#' @examples
#' set.seed(123)
#' D <- seq(10, 50, by = 5)
#' H <- seq(15, 30, length.out = length(D))
#' V <- 0.01 * D^2 * H + rnorm(length(D), sd = 5)
#' V[c(3, 7)] <- V[c(3, 7)] + 50  # add outliers
#' data <- data.frame(D = D, H = H, V = V)
#'
#' model <- ransac_nls(V ~ a * D^b * H^c, data = data,
#'                     start = list(a = 0.01, b = 2, c = 1),
#'                     n_min = 4, n_iter = 100, tol = 10)
#' summary(model)
#'
#' @export
#' @importFrom stats nls predict nls.control
ransac_nls <- function(formula, data, start, n_min, n_iter = 100, tol = 0.2, verbose = FALSE) {
  best_inliers <- c()
  best_model <- NULL

  for (i in 1:n_iter) {
    sample_idx <- sample(1:nrow(data), n_min)
    sample_data <- data[sample_idx, ]

    model_try <- try(nls(formula, data = sample_data, start = start, control = nls.control(warnOnly = TRUE)), silent = TRUE)
    if (inherits(model_try, "try-error")) next

    pred <- try(predict(model_try, newdata = data), silent = TRUE)
    if (inherits(pred, "try-error")) next

    y_var <- all.vars(formula)[1]
    error <- abs(data[[y_var]] - pred)
    inliers <- which(error < tol)

    if (length(inliers) > length(best_inliers)) {
      best_inliers <- inliers
      best_model <- model_try
      if (verbose) message("Iteration ", i, ": ", length(inliers), " inliers")
    }
  }

  final_model <- nls(formula, data = data[best_inliers, ], start = start, control = nls.control(warnOnly = TRUE))
  class(final_model) <- c("ransac_nls", class(final_model))
  attr(final_model, "inliers") <- best_inliers
  return(final_model)
}