多站点负指数函数全局参数优化（R语言实现）

Got it, let's walk through how to optimize those k and k2 parameters for your double negative exponential model in R. Here's a step-by-step approach tailored to your dataset:

Using cbind creates a matrix with mixed data types (factors + numbers), which can cause issues down the line. Let's convert it to a proper data frame instead:

# Build a clean data frame (way easier to work with than a mixed matrix)
data_n <- data.frame(
  Study = as.factor(c(1, 1, 1, 1, 2, 2, 2, 3, 3)),
  Time = as.numeric(c(0, 0.08, 0.16, 0.24, 0, 0.05, 0.88, 0, 0.99)),
  Remaining = as.numeric(c(100, 80, 69, 45, 100, 60, 35, 0, 25))
)

We want to find the k and k2 values that make our predicted Remaining values as close as possible to the observed data. The standard way to do this is to minimize the residual sum of squares (RSS)—the total squared difference between predicted and actual values. Here's the function:

# Calculate RSS for given k and k2 parameters
rss_function <- function(params, data) {
  k <- params[1]
  k2 <- params[2]
  
  # Your double negative exponential model (41 is fixed as you specified)
  predicted <- 41 * exp(-k * data$Time) + 59 * exp(-k2 * data$Time)
  
  # Sum of squared residuals
  rss <- sum((data$Remaining - predicted)^2)
  return(rss)
}

R's built-in optim() function is perfect for this kind of unconstrained (or constrained) minimization. We just need to start with an initial guess for k and k2—pick values that make sense for your data (e.g., small positive numbers since we're dealing with decay):

# Initial guesses for k and k2 (adjust these if optimization doesn't converge)
initial_guess <- c(k = 1, k2 = 0.5)

# Run the optimization
optim_results <- optim(
  par = initial_guess,
  fn = rss_function,
  data = data_n,
  method = "L-BFGS-B",  # Good for bounded parameters
  lower = c(0, 0)       # Ensure k and k2 are positive (decay can't have negative rates!)
)

Once the optimization finishes, you can pull out the optimal parameters and the minimum RSS:

# Print the optimized values
cat("Optimal k:", round(optim_results$par[1], 3), "\n")
cat("Optimal k2:", round(optim_results$par[2], 3), "\n")
cat("Minimum RSS:", round(optim_results$value, 2), "\n")

To make sure the model fits your data, plot the observed vs. predicted values:

# Generate predicted values using the optimal parameters
data_n$Predicted <- 41 * exp(-optim_results$par[1] * data_n$Time) + 
  59 * exp(-optim_results$par[2] * data_n$Time)

# Plot observed points (colored by study) and predicted line
plot(data_n$Time, data_n$Remaining, 
     col = data_n$Study, pch = 16,
     xlab = "Time", ylab = "Remaining",
     main = "Observed vs. Predicted Remaining Values")
lines(data_n$Time, data_n$Predicted, col = "red", lwd = 2)
legend("topright", legend = levels(data_n$Study), col = 1:3, pch = 16)

• Initial guesses matter: If the optimization doesn't converge (check optim_results$convergence—a value of 0 means success), try different starting values for k and k2.
• Parameter bounds: We added lower = c(0,0) because decay rates can't be negative—this prevents nonsensical results.
• Model flexibility: If you ever need to adjust the fixed 41 constant later, you can add it as a third parameter to optimize (just update the function and initial guess).

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