基于一阶差分方程的参数拟合方法咨询及优化方案问询

Hey there! As someone new to fitting recursive models, it’s totally normal to wonder about the right tools and approaches. Let’s tackle your questions one by one:

1. Is using nls.lm appropriate for this problem? #

Absolutely! nls.lm is built for nonlinear least squares optimization, which is exactly what you’re doing here—minimizing the gap between your recursively computed y2_mod and experimental y2_exp. The Levenberg-Marquardt algorithm it uses is robust for most nonlinear fitting scenarios, especially when you can define a clear residual function that calculates y2_mod via your recursive equation and returns the residuals (y2_mod - y2_exp).

The key here is ensuring your residual function correctly implements the full recursive sequence using y1_exp—as long as that’s solid, nls.lm is a perfectly valid and reliable choice.

2. Are there better curve fitting/parameter optimization methods? #

It depends on your data characteristics and goals:

• If residuals are normally distributed: Stick with nls.lm—least squares is optimal under the Gauss-Markov theorem, and it’s computationally efficient.
• If you have outliers or non-normal residuals: Try robust nonlinear fitting with the robustbase package’s nlrob function, which downweights outliers to avoid skewed parameter estimates.
• If you want parameter uncertainty or prior knowledge: Go with Bayesian methods using packages like brms or rstan. These let you specify prior beliefs about parameter values (e.g., "a should be between 0 and 1") and generate posterior distributions that show how certain you can be about your fitted parameters.
• If your model has many local optima: Use global optimization methods like genetic algorithms (GA package) or particle swarm optimization (pso package). These explore more of the parameter space to avoid getting stuck in suboptimal solutions, though they’re slower than Levenberg-Marquardt.

3. Better code approaches or packages for recursive models needing the full y1_exp sequence? #

You’re right that deSolve is focused on initial-value problems for differential equations, but there are flexible ways to handle your recursive model without specialized packages:

Custom Recursive Function (Most Flexible) #

Write a simple function to compute y2_mod using a loop or vectorized tools like Reduce (faster for long sequences). Here’s an example:

# Define the recursive computation
compute_y2_mod <- function(a, fixed_b, y1_exp, init_y2) {
  n <- length(y1_exp)
  y2_mod <- numeric(n)
  y2_mod[1] <- init_y2
  
  # Loop through the sequence
  for (i in 2:n) {
    y2_mod[i] <- a * y2_mod[i-1] + fixed_b * y1_exp[i]
  }
  return(y2_mod)
}

# Residual function for nls.lm
resid_fun <- function(params, y1_exp, y2_exp, init_y2, fixed_b) {
  a <- params[1]
  y2_mod <- compute_y2_mod(a, fixed_b, y1_exp, init_y2)
  return(y2_mod - y2_exp)
}

# Example usage with nls.lm
library(minpack.lm)
initial_guess <- c(a = 0.5) # Starting value for a
fit <- nls.lm(par = initial_guess, 
              fn = resid_fun, 
              y1_exp = y1_exp, 
              y2_exp = y2_exp, 
              init_y2 = y2_exp[1], # Use first experimental value as initial guess
              fixed_b = 1.2) # Your fixed b value

For faster computation with long sequences, replace the loop with Reduce:

compute_y2_mod <- function(a, fixed_b, y1_exp, init_y2) {
  Reduce(function(prev, curr) a * prev + fixed_b * curr, 
         y1_exp[-1], 
         init = init_y2, 
         accumulate = TRUE)
}

Using deSolve for Discrete Models #

You can also adapt deSolve to handle your difference equation by treating it as a discrete-time system:

library(deSolve)
diff_model <- function(t, y, params) {
  a <- params["a"]
  fixed_b <- params["fixed_b"]
  # t is 0-indexed, so y1_exp[t+1] is the next input value
  dy <- a * y + fixed_b * y1_exp[t+1] - y # y(t+1) - y(t) = dy
  return(list(dy))
}

times <- 0:(length(y1_exp)-1) # Match length of y1_exp
init_conditions <- c(y = y2_exp[1])
params <- c(a = 0.5, fixed_b = 1.2)

# Solve discrete model
y2_mod <- ode(y = init_conditions, 
              times = times, 
              func = diff_model, 
              parms = params, 
              method = "euler", 
              discrete = TRUE)[, "y"]

While there’s no one-size-fits-all package for every recursive nonlinear model, custom functions give you full control over how y1_exp is integrated into the recursion—this is usually the most practical approach for domain-specific models.

内容的提问来源于stack exchange，提问作者Mike