R/fit_distributions.R

In dendrometry: Forest Estimations and Dendrometric Computations

#' Unified Parameter Estimation for Probability Distributions
#'
#' @description
#' Estimates parameters of probability distributions using various methods:
#' Maximum Likelihood (MLE), Maximum Product Spacing (MPS),
#' or Method of Moments (MOM).
#'
#' @param data numeric vector of observed data.
#' @param dist character string specifying the distribution. Options include:
#'   "normal", "exponential", "gamma", "weibull", "weibull3", "lognormal",
#'   "lognormal3", "beta", or "custom".
#' @param method character string specifying estimation method. Options:
#'   "mle" (Maximum Likelihood - default), "mps" (Maximum Product Spacing),
#'   "mom" (Method of Moments).
#' @param start named list or numeric vector of initial parameter values.
#'   Required for "custom" distributions with "mps" or "mle" methods.
#' @param lower named list or numeric vector of lower bounds for parameters.
#' @param upper named list or numeric vector of upper bounds for parameters.
#' @param optim_method optimization method passed to \code{\link{optim}}. Default is "Nelder-Mead".
#' @param custom_functions list containing custom distribution functions (for dist="custom"):
#'   \itemize{
#'     \item \code{pdf}: probability density function f(x, params) (for MLE and ties correction)
#'     \item \code{cdf}: cumulative distribution function F(x, params) (for MPS)
#'     \item \code{param_names}: character vector of parameter names
#'     \item \code{start_fn}: function to compute starting values (optional)
#'   }
#' @param tol_spacing numeric tolerance for spacings/densities to avoid log(0). Default is 1e-16.
#' @param tol_param numeric tolerance for parameter lower bounds. Default is 1e-6.
#' @param ties_method character string for handling ties in MPS. Options:
#'   "cheng_amin" (default - Cheng & Amin 1983), "none", "cheng_stephens" (Cheng & Stephens 1989).
#'   Only applicable when method = "mps".
#' @param ... Additional arguments passed to \code{\link{optim}}.
#'
#' @return A list with class "fit_dist" containing:
#' \item{estimate}{Named vector of parameter estimates}
#' \item{vcov}{Variance-covariance matrix (for mps/mle)}
#' \item{se}{Standard errors (for mps/mle)}
#' \item{loglik}{Log-likelihood value}
#' \item{aic}{Akaike Information Criterion}
#' \item{bic}{Bayesian Information Criterion}
#' \item{objective}{Maximum value of objective function (logspacing for mps, loglik for mle)}
#' \item{ks_statistic}{Kolmogorov-Smirnov test statistic}
#' \item{ks_pvalue}{Kolmogorov-Smirnov test p-value}
#' \item{convergence}{Convergence code from optim (0 indicates success)}
#' \item{message}{Convergence message from optim}
#' \item{data}{Original data (sorted)}
#' \item{dist}{Distribution name}
#' \item{method}{Estimation method used}
#' \item{n}{Sample size}
#' \item{k}{Number of parameters}
#' \item{tol_spacing}{Tolerance used for spacings/densities}
#' \item{tol_param}{Tolerance used for parameter bounds}
#' \item{ties_method}{Ties correction method (for MPS only)}
#'
#' @references
#' Cheng, R. C. H., & Amin, N. A. K. (1983). Estimating parameters in continuous
#' univariate distributions with a shifted origin. Journal of the Royal Statistical
#' Society: Series B, 45(3), 394-403.
#'
#' Cheng, R. C. H., & Stephens, M. A. (1989). A goodness-of-fit test using Moran's
#' statistic with estimated parameters. Biometrika, 76(2), 385-392.
#'
#' @examples
#' # MLE estimation (default)
#' set.seed(123)
#' x <- rweibull(100, shape = 2.5, scale = 1.5)
#' fit1 <- fit_dist(x, dist = "weibull")
#' print(fit1)
#'
#' # MPS estimation with Cheng-Amin ties correction
#' fit2 <- fit_dist(x, dist = "weibull", method = "mps")
#'
#' # L-Moments estimation (under development)
#' # fit3 <- fit_dist(x, dist = "weibull", method = "lm")
#'
#' # Method of Moments
#' fit4 <- fit_dist(x, dist = "weibull", method = "mom")
#'
#' # Compare fits
#' cat("AIC - MLE:", fit1$aic, "MPS:", fit2$aic, "\n")
#'
#' @export
fit_dist <- function(data, dist = "normal", method = "mle",
                     start = NULL, lower = NULL, upper = NULL,
                     optim_method = "Nelder-Mead", custom_functions = NULL,
                     tol_spacing = 1e-16, tol_param = 1e-6,
                     ties_method = "cheng_amin", ...) {
  # Input validation
  if (!is.numeric(data) || length(data) < 2) {
    stop("'data' must be a numeric vector with at least 2 observations.")
  }

  if (!method %in% c("mps", "mle", "lm", "mom")) {
    stop("'method' must be one of: 'mle', 'mps', 'lm', 'mom'")
  }

  if (!ties_method %in% c("none", "cheng_amin", "cheng_stephens")) {
    stop("'ties_method' must be one of: 'none', 'cheng_amin', 'cheng_stephens'")
  }

  if (any(is.na(data))) {
    warning("NA values removed from data")
    data <- data[!is.na(data)]
  }

  # Sort data
  x <- sort(data)
  n <- length(x)

  # Route to appropriate estimation method
  result <- switch(method,
    mps = ,
    mle = fit_optimization(
      x, dist, method, start, lower, upper,
      optim_method, custom_functions,
      tol_spacing, tol_param, ties_method, ...
    ),
    lm = fit_lmoments(x, dist),
    mom = fit_moments(x, dist, tol_param),
    stop("Invalid method")
  )

  # Add common elements
  result$data <- x
  result$dist <- dist
  result$method <- method
  result$n <- n
  result$k <- length(result$estimate)
  result$tol_spacing <- tol_spacing
  result$tol_param <- tol_param
  result$ties_method <- if (method == "mps") ties_method else NA

  # Clean up for non-optimization methods
  if (method %in% c("lm", "mom")) {
    result$vcov <- NULL
    result$se <- NULL
    result$objective <- NA
    result$convergence <- NA
    result$message <- NULL
  }

  class(result) <- "fit_dist"
  return(result)
}

#' Unified Optimization-based Estimation (MPS and MLE)
#'
#' @description
#' Internal function that performs parameter estimation using optimization-based
#' methods (Maximum Likelihood or Maximum Product Spacing).
#'
#' @param x numeric vector of sorted data.
#' @param dist character string specifying the distribution.
#' @param method character string: "mle" or "mps".
#' @param start initial parameter values.
#' @param lower lower bounds for parameters.
#' @param upper upper bounds for parameters.
#' @param optim_method optimization method for \code{\link{optim}}.
#' @param custom_functions list of custom distribution functions.
#' @param tol_spacing numeric tolerance for spacings/densities.
#' @param tol_param numeric tolerance for parameter bounds.
#' @param ties_method character string for ties correction method.
#' @param ... additional arguments for \code{\link{optim}}.
#'
#' @return A list containing parameter estimates, standard errors, objective value,
#'   log-likelihood, AIC, BIC, KS test results, and convergence information.
#'
#' @keywords internal
fit_optimization <- function(x, dist, method, start, lower, upper, optim_method,
                             custom_functions, tol_spacing, tol_param, ties_method, ...) {
  n <- length(x)

  # Validate custom distributions
  if (dist == "custom") {
    if (is.null(custom_functions)) {
      stop("For custom distributions, provide 'custom_functions'")
    }

    if (method == "mps" && is.null(custom_functions$cdf)) {
      stop("For custom distributions with MPS, provide 'cdf' in custom_functions")
    }
    if (method == "mle" && is.null(custom_functions$pdf)) {
      stop("For custom distributions with MLE, provide 'pdf' in custom_functions")
    }

    # Check for PDF when using Cheng-Amin ties correction
    if (method == "mps" && ties_method == "cheng_amin" && is.null(custom_functions$pdf)) {
      stop("For custom distributions with ties_method='cheng_amin', provide 'pdf' in custom_functions")
    }

    param_names <- custom_functions$param_names
  } else {
    param_names <- get_param_names(dist)
  }

  # Create objective function
  obj_func <- create_objective(method, x, dist, tol_spacing, ties_method, custom_functions)

  # Get starting values
  if (is.null(start)) {
    if (dist == "custom" && !is.null(custom_functions$start_fn)) {
      start <- custom_functions$start_fn(x)
    } else {
      start <- get_start_values(x, dist, tol_param)
    }
  } else if (is.list(start)) {
    start <- unlist(start)
  }

  # Perform optimization
  opt_result <- optimize_params(
    obj_func, start, lower, upper, optim_method,
    dist, x, tol_param, ...
  )

  names(opt_result$par) <- param_names

  # Compute variance-covariance matrix and standard errors with robust error handling
  vcov_mat <- se <- NULL
  if (!is.null(opt_result$hessian)) {
    vcov_mat <- tryCatch(solve(opt_result$hessian), error = function(e) NULL)

    if (!is.null(vcov_mat)) {
      se <- tryCatch(
        {
          se_vals <- sqrt(diag(vcov_mat))
          # Check for NaN (negative variance estimates)
          if (any(is.nan(se_vals))) {
            warning("Negative variance estimates detected; SE set to NULL")
            NULL
          } else {
            names(se_vals) <- param_names
            se_vals
          }
        },
        error = function(e) {
          warning("Error computing standard errors: ", e$message)
          NULL
        }
      )

      if (!is.null(vcov_mat)) {
        colnames(vcov_mat) <- rownames(vcov_mat) <- param_names
      }
    }
  }

  # Compute log-likelihood
  loglik <- sum(log(pmax(get_pdf(x, opt_result$par, dist, custom_functions), tol_spacing)))

  # Compute information criteria
  k <- length(opt_result$par)
  aic <- -2 * loglik + 2 * k
  bic <- -2 * loglik + k * log(n)

  # Perform Kolmogorov-Smirnov test
  ks_result <- tryCatch(
    {
      ks_test <- ks.test(x, function(q) get_cdf(q, opt_result$par, dist, custom_functions))
      list(statistic = as.numeric(ks_test$statistic), p.value = ks_test$p.value)
    },
    error = function(e) {
      warning("KS test failed: ", e$message)
      list(statistic = NA, p.value = NA)
    }
  )

  list(
    estimate = opt_result$par,
    vcov = vcov_mat,
    se = se,
    loglik = loglik,
    aic = aic,
    bic = bic,
    objective = -opt_result$value,
    ks_statistic = ks_result$statistic,
    ks_pvalue = ks_result$p.value,
    convergence = opt_result$convergence,
    message = opt_result$message
  )
}

#' Create Objective Function for MPS or MLE
#'
#' @description
#' Creates the objective function for optimization-based estimation methods.
#' For MLE, returns negative log-likelihood. For MPS, returns negative log-product-spacing
#' with optional ties correction.
#'
#' @param method character string: "mle" or "mps".
#' @param x numeric vector of sorted data.
#' @param dist character string specifying the distribution.
#' @param tol_spacing numeric tolerance for spacings/densities.
#' @param ties_method character string for ties correction.
#' @param custom_functions list of custom distribution functions.
#'
#' @return A function that computes the objective value given parameters.
#'
#' @keywords internal
create_objective <- function(method, x, dist, tol_spacing, ties_method, custom_functions) {
  n <- length(x)

  if (method == "mle") {
    # MLE objective: negative log-likelihood
    return(function(params) {
      densities <- get_pdf(x, params, dist, custom_functions)
      densities <- pmax(densities, tol_spacing)
      return(-sum(log(densities)))
    })
  } else {
    # MPS objective: negative log-product-spacing
    if (ties_method == "none") {
      return(function(params) {
        F_vals <- get_cdf(x, params, dist, custom_functions)
        F_vals <- c(0, F_vals, 1)
        spacings <- diff(F_vals)
        spacings <- pmax(spacings, tol_spacing)
        return(-sum(log(spacings)))
      })
    } else if (ties_method == "cheng_amin") {
      # Cheng & Amin (1983): Replace zero spacings with PDF values
      return(function(params) {
        F_vals <- get_cdf(x, params, dist, custom_functions)
        F_vals <- c(0, F_vals, 1)
        spacings <- diff(F_vals)

        # Identify tied observations (zero or near-zero spacings)
        tied_idx <- which(spacings < tol_spacing)

        if (length(tied_idx) > 0) {
          # Get PDF values for tied spacing correction
          pdf_vals <- get_pdf(x, params, dist, custom_functions)

          # For tied observations, use PDF values instead of spacings
          for (idx in tied_idx) {
            if (idx <= n) { # Not the final spacing
              spacings[idx] <- max(pdf_vals[idx], tol_spacing)
            }
          }
        }

        spacings <- pmax(spacings, tol_spacing)
        return(-sum(log(spacings)))
      })
    } else if (ties_method == "cheng_stephens") {
      # Cheng & Stephens (1989): Modified spacing calculation
      return(function(params) {
        F_vals <- get_cdf(x, params, dist, custom_functions)

        # Calculate spacings with Moran's statistic adjustment
        spacings <- numeric(n + 1)
        spacings[1] <- F_vals[1]
        spacings[n + 1] <- 1 - F_vals[n]

        for (i in 2:n) {
          spacings[i] <- F_vals[i] - F_vals[i - 1]

          # Apply Cheng & Stephens correction for small spacings
          if (spacings[i] < tol_spacing) {
            # Use average of adjacent spacings
            spacings[i] <- mean(c(
              spacings[i - 1],
              if (i < n) F_vals[i + 1] - F_vals[i] else 1 - F_vals[i]
            ))
          }
        }

        spacings <- pmax(spacings, tol_spacing)
        return(-sum(log(spacings)))
      })
    }
  }
}

#' L-Moments Estimation
#'
#' @description
#' Estimates distribution parameters using the method of L-moments.
#'
#' @param x numeric vector of sorted data.
#' @param dist character string specifying the distribution.
#'
#' @return A list containing parameter estimates and placeholder values for
#'   other components.
#'
#' @keywords internal
fit_lmoments <- function(x, dist) {
  stop("L-moments estimation not yet implemented (under development)...")
}

#' Method of Moments Estimation
#'
#' @description
#' Estimates distribution parameters using the classical method of moments,
#' matching sample moments to theoretical moments.
#'
#' @param x numeric vector of sorted data.
#' @param dist character string specifying the distribution.
#' @param tol_param numeric tolerance for parameter bounds.
#'
#' @return A list containing parameter estimates and placeholder values for
#'   other components.
#'
#' @keywords internal
fit_moments <- function(x, dist, tol_param) {
  m1 <- mean(x)
  m2 <- var(x)

  estimates <- switch(dist,
    normal = c(mean = m1, sd = sqrt(m2)),
    exponential = c(rate = 1 / m1),
    gamma = c(shape = m1^2 / m2, rate = m1 / m2),
    weibull = {
      shape <- m1 / sqrt(m2)
      scale <- m1 / gamma(1 + 1 / shape)
      c(shape = shape, scale = scale)
    },
    lognormal = {
      mu <- log(m1^2 / sqrt(m2 + m1^2))
      sigma <- sqrt(log(1 + m2 / m1^2))
      c(meanlog = mu, sdlog = sigma)
    },
    lognormal3 = {
      loc <- min(x) - 0.01 * diff(range(x))
      x_shifted <- x - loc
      m1_s <- mean(x_shifted)
      m2_s <- var(x_shifted)
      mu <- log(m1_s^2 / sqrt(m2_s + m1_s^2))
      sigma <- sqrt(log(1 + m2_s / m1_s^2))
      c(meanlog = mu, sdlog = sigma, loc = loc)
    },
    beta = {
      a <- m1 * (m1 * (1 - m1) / m2 - 1)
      b <- (1 - m1) * (m1 * (1 - m1) / m2 - 1)
      c(shape1 = max(a, tol_param), shape2 = max(b, tol_param))
    },
    stop("Method of moments not available for this distribution")
  )

  # Compute log-likelihood, AIC, BIC for comparison
  n <- length(x)
  k <- length(estimates)
  loglik <- sum(log(pmax(get_pdf(x, estimates, dist, NULL), 1e-16)))
  aic <- -2 * loglik + 2 * k
  bic <- -2 * loglik + k * log(n)

  # Perform KS test
  ks_result <- tryCatch(
    {
      ks_test <- ks.test(x, function(q) get_cdf(q, estimates, dist, NULL))
      list(statistic = as.numeric(ks_test$statistic), p.value = ks_test$p.value)
    },
    error = function(e) {
      list(statistic = NA, p.value = NA)
    }
  )

  list(
    estimate = estimates, vcov = NULL, se = NULL,
    loglik = loglik, aic = aic, bic = bic, objective = NA,
    ks_statistic = ks_result$statistic, ks_pvalue = ks_result$p.value,
    convergence = 0, message = "Method of moments estimation"
  )
}

#' Optimization Wrapper
#'
#' @description
#' Wrapper function for \code{\link{optim}} that handles bounded and unbounded
#' optimization methods with appropriate default bounds.
#'
#' @param obj_func objective function to minimize.
#' @param start initial parameter values.
#' @param lower lower bounds for parameters.
#' @param upper upper bounds for parameters.
#' @param optim_method optimization method.
#' @param dist character string specifying the distribution.
#' @param x numeric vector of data.
#' @param tol_param numeric tolerance for parameter bounds.
#' @param ... additional arguments for \code{\link{optim}}.
#'
#' @return Result from \code{\link{optim}}.
#'
#' @keywords internal
optimize_params <- function(obj_func, start, lower, upper, optim_method, dist, x, tol_param, ...) {
  # Set default bounds for bounded methods
  if (optim_method %in% c("L-BFGS-B", "Brent")) {
    if (is.null(lower)) {
      lower <- get_default_bounds(dist, "lower", x, tol_param)
    } else if (is.list(lower)) {
      lower <- unlist(lower)
    }

    if (is.null(upper)) {
      upper <- get_default_bounds(dist, "upper", x, tol_param)
    } else if (is.list(upper)) {
      upper <- unlist(upper)
    }

    opt_result <- optim(
      par = start, fn = obj_func, method = optim_method,
      hessian = TRUE, lower = lower, upper = upper, ...
    )
  } else {
    opt_result <- optim(
      par = start, fn = obj_func, method = optim_method,
      hessian = TRUE, ...
    )
  }

  opt_result
}

#' Get PDF for Standard and Custom Distributions
#'
#' @description
#' Evaluates the probability density function for standard distributions
#' or uses custom PDF function.
#'
#' @param x numeric vector of quantiles.
#' @param params numeric vector of parameter values.
#' @param dist character string specifying the distribution.
#' @param custom_functions list of custom distribution functions (optional).
#'
#' @return Numeric vector of density values.
#'
#' @keywords internal
get_pdf <- function(x, params, dist, custom_functions = NULL) {
  if (dist == "custom") {
    if (is.null(custom_functions) || is.null(custom_functions$pdf)) {
      stop("PDF function not available for custom distribution")
    }
    return(custom_functions$pdf(x, params))
  }

  switch(dist,
    normal = dnorm(x, mean = params[1], sd = params[2]),
    exponential = dexp(x, rate = params[1]),
    gamma = dgamma(x, shape = params[1], rate = params[2]),
    weibull = dweibull(x, shape = params[1], scale = params[2]),
    weibull3 = dweibull(x - params[3], shape = params[1], scale = params[2]),
    lognormal = dlnorm(x, meanlog = params[1], sdlog = params[2]),
    lognormal3 = dlnorm(x - params[3], meanlog = params[1], sdlog = params[2]),
    beta = dbeta(x, shape1 = params[1], shape2 = params[2]),
    stop("Distribution not supported")
  )
}

#' Get CDF for Standard and Custom Distributions
#'
#' @description
#' Evaluates the cumulative distribution function for standard distributions
#' or uses custom CDF function.
#'
#' @param x numeric vector of quantiles.
#' @param params numeric vector of parameter values.
#' @param dist character string specifying the distribution.
#' @param custom_functions list of custom distribution functions (optional).
#'
#' @return Numeric vector of probabilities.
#'
#' @keywords internal
get_cdf <- function(x, params, dist, custom_functions = NULL) {
  if (dist == "custom") {
    if (is.null(custom_functions) || is.null(custom_functions$cdf)) {
      stop("CDF function not available for custom distribution")
    }
    return(custom_functions$cdf(x, params))
  }

  switch(dist,
    normal = pnorm(x, mean = params[1], sd = params[2]),
    exponential = pexp(x, rate = params[1]),
    gamma = pgamma(x, shape = params[1], rate = params[2]),
    weibull = pweibull(x, shape = params[1], scale = params[2]),
    weibull3 = pweibull(x - params[3], shape = params[1], scale = params[2]),
    lognormal = plnorm(x, meanlog = params[1], sdlog = params[2]),
    lognormal3 = plnorm(x - params[3], meanlog = params[1], sdlog = params[2]),
    beta = pbeta(x, shape1 = params[1], shape2 = params[2]),
    stop("Distribution not supported")
  )
}

#' Get Starting Values for Parameter Estimation
#'
#' @description
#' Computes reasonable starting values for optimization based on sample moments
#' and distribution-specific heuristics.
#'
#' @param x numeric vector of sorted data.
#' @param dist character string specifying the distribution.
#' @param tol_param numeric tolerance for parameter bounds.
#'
#' @return Numeric vector of starting parameter values.
#'
#' @keywords internal
get_start_values <- function(x, dist, tol_param) {
  m <- mean(x)
  v <- var(x)

  switch(dist,
    normal = c(m, sd(x)),
    exponential = c(1 / m),
    gamma = c(m^2 / v, m / v),
    weibull = c(1.5, m),
    weibull3 = c(
      1.5, m - min(x) + 0.01 * diff(range(x)),
      min(x) - 0.01 * diff(range(x))
    ),
    lognormal = c(mean(log(x)), sd(log(x))),
    lognormal3 = {
      loc <- min(x) - 0.01 * diff(range(x))
      x_shifted <- x - loc
      c(mean(log(x_shifted)), sd(log(x_shifted)), loc)
    },
    beta = {
      a <- m * (m * (1 - m) / v - 1)
      b <- (1 - m) * (m * (1 - m) / v - 1)
      c(max(a, tol_param), max(b, tol_param))
    },
    stop("Distribution not supported")
  )
}

#' Get Default Parameter Bounds
#'
#' @description
#' Returns default lower or upper bounds for distribution parameters
#' used in bounded optimization methods.
#'
#' @param dist character string specifying the distribution.
#' @param bound_type character string: "lower" or "upper".
#' @param x numeric vector of data.
#' @param tol_param numeric tolerance for parameter bounds.
#'
#' @return Numeric vector of parameter bounds.
#'
#' @keywords internal
get_default_bounds <- function(dist, bound_type, x, tol_param) {
  if (bound_type == "lower") {
    switch(dist,
      normal = c(-Inf, tol_param),
      exponential = c(tol_param),
      gamma = c(tol_param, tol_param),
      weibull = c(tol_param, tol_param),
      weibull3 = c(tol_param, tol_param, -Inf),
      lognormal = c(-Inf, tol_param),
      lognormal3 = c(-Inf, tol_param, -Inf),
      beta = c(tol_param, tol_param),
      c(-Inf)
    )
  } else {
    switch(dist,
      normal = c(Inf, Inf),
      exponential = c(Inf),
      gamma = c(Inf, Inf),
      weibull = c(Inf, Inf),
      weibull3 = c(Inf, Inf, min(x) - tol_param),
      lognormal = c(Inf, Inf),
      lognormal3 = c(Inf, Inf, min(x) - tol_param),
      beta = c(Inf, Inf),
      c(Inf)
    )
  }
}

#' Get Parameter Names for Standard Distributions
#'
#' @description
#' Returns the parameter names for standard distributions.
#'
#' @param dist character string specifying the distribution.
#'
#' @return Character vector of parameter names.
#'
#' @keywords internal
get_param_names <- function(dist) {
  switch(dist,
    normal = c("mean", "sd"),
    exponential = c("rate"),
    gamma = c("shape", "rate"),
    weibull = c("shape", "scale"),
    weibull3 = c("shape", "scale", "loc"),
    lognormal = c("meanlog", "sdlog"),
    lognormal3 = c("meanlog", "sdlog", "loc"),
    beta = c("shape1", "shape2"),
    stop("Distribution not supported")
  )
}

#' Extract Variance-Covariance Matrix
#'
#' @description
#' S3 method to extract the variance-covariance matrix from a fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments (not used).
#'
#' @return Variance-covariance matrix of parameter estimates.
#'
#' @export
vcov.fit_dist <- function(object, ...) {
  if (is.null(object$vcov)) {
    stop("Variance-covariance matrix not available for this method")
  }
  object$vcov
}

#' Extract Model Coefficients
#'
#' @description
#' S3 method to extract parameter estimates from a fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments (not used).
#'
#' @return Named numeric vector of parameter estimates.
#'
#' @export
coef.fit_dist <- function(object, ...) {
  object$estimate
}

#' Confidence Intervals for Parameters
#'
#' @description
#' S3 method to compute confidence intervals for distribution parameters
#' based on asymptotic normality of maximum likelihood estimates.
#'
#' @param object an object of class "fit_dist".
#' @param parm character vector of parameter names or numeric vector of indices.
#'   If missing, all parameters are considered.
#' @param level confidence level (default: 0.95).
#' @param ... additional arguments (not used).
#'
#' @return A matrix with columns giving lower and upper confidence limits
#'   for each parameter.
#'
#' @export
confint.fit_dist <- function(object, parm, level = 0.95, ...) {
  if (is.null(object$se)) {
    warning("Standard errors not available; confidence intervals cannot be computed")
    return(NULL)
  }

  cf <- coef(object)
  pnames <- names(cf)
  if (missing(parm)) {
    parm <- pnames
  } else if (is.numeric(parm)) parm <- pnames[parm]

  a <- (1 - level) / 2
  a <- c(a, 1 - a)
  pct <- paste(format(100 * a, trim = TRUE, scientific = FALSE, digits = 3), "%")

  ci <- array(NA,
    dim = c(length(parm), 2L),
    dimnames = list(parm, pct)
  )

  ses <- object$se[parm]
  ci[] <- cf[parm] + ses %o% qnorm(a)
  ci
}

#' Extract Log-Likelihood
#'
#' @description
#' S3 method to extract the log-likelihood value from a fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments (not used).
#'
#' @return Log-likelihood value.
#'
#' @export
logLik.fit_dist <- function(object, ...) {
  if (is.na(object$loglik)) {
    stop("Log-likelihood not available")
  }
  val <- object$loglik
  attr(val, "df") <- object$k
  attr(val, "nobs") <- object$n
  class(val) <- "logLik"
  val
}

#' Extract AIC
#'
#' @description
#' S3 method to extract the Akaike Information Criterion from a fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments (not used).
#' @param k penalty per parameter (default: 2 for AIC).
#'
#' @return AIC value.
#'
#' @export
AIC.fit_dist <- function(object, ..., k = 2) {
  if (is.na(object$loglik)) {
    stop("AIC not available")
  }
  -2 * object$loglik + k * object$k
}

#' Extract BIC
#'
#' @description
#' S3 method to extract the Bayesian Information Criterion from a fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments (not used).
#'
#' @return BIC value.
#'
#' @export
BIC.fit_dist <- function(object, ...) {
  if (is.na(object$loglik)) {
    stop("BIC not available")
  }
  object$bic
}

#' Print Method for Fitted Distributions
#'
#' @description
#' S3 method to print a summary of the fitted distribution.
#'
#' @param x an object of class "fit_dist".
#' @param ... additional arguments (not used).
#'
#' @return The object invisibly.
#'
#' @export
print.fit_dist <- function(x, ...) {
  cat("Distribution Parameter Estimation\n")
  cat("==================================\n")
  cat("Method:", toupper(x$method), "\n")
  cat("Distribution:", x$dist, "\n")
  cat("Sample size:", x$n, "\n")
  if (!is.na(x$ties_method)) {
    cat("Ties method:", x$ties_method, "\n")
  }
  cat("\n")

  cat("Parameter estimates:\n")
  if (!is.null(x$se)) {
    ci <- confint(x)
    out <- cbind(Estimate = x$estimate, SE = x$se, ci)
    print(out)
  } else {
    print(x$estimate)
  }
  cat("\n")

  if (!is.na(x$loglik)) {
    cat("Log-likelihood:", round(x$loglik, digits = 3), "\n")
  }

  if (!is.na(x$aic)) {
    cat("AIC:", round(x$aic, digits = 3), "  ")
    cat("BIC:", round(x$bic, digits = 3), "\n")
  }

  if (!is.na(x$objective)) {
    obj_name <- if (x$method == "mps") "Log-spacing" else "Log-likelihood (objective)"
    cat(obj_name, ":", round(x$objective, digits = 3), "\n")
  }

  if (!is.na(x$ks_statistic)) {
    cat("\nKolmogorov-Smirnov goodness-of-fit test:\n")
    cat("  D-statistic =", round(x$ks_statistic, digits = 4), "\n")
    cat("  p-value     =", format.pval(x$ks_pvalue, digits = 4), "\n")
  }

  if (!is.null(x$convergence) && !is.na(x$convergence)) {
    cat("\nConvergence:", ifelse(x$convergence == 0, "successful",
      paste("failed (code", x$convergence, ")")
    ), "\n")
  }

  invisible(x)
}

#' Summary Method for Fitted Distributions
#'
#' @description
#' S3 method to provide a detailed summary of the fitted distribution.
#'
#' @param object an object of class "fit_dist".
#' @param ... additional arguments passed to \code{\link{print.fit_dist}}.
#'
#' @return The object invisibly.
#'
#' @export
summary.fit_dist <- function(object, ...) {
  print(object, ...)
}