如何用R语言optim()函数求解参数？结果异常求助

Let’s break down why your current optim() call is giving wrong results and walk through the correct approach to estimate your alpha parameter.

What’s Wrong With Your Original Code? #

Your implementation has three critical issues that throw off the optimization:

1. Bad loss function choice: You’re returning -sum(rate - pred) as the objective. This is equivalent to maximizing sum(rate - pred), which doesn’t have a meaningful minimum (e.g., as alpha grows to infinity, pred approaches 0, and sum(rate - pred) flattens to a constant). The BFGS algorithm can’t converge to a valid solution here.
2. Extra parameter dimension: You passed an initial value c(0,0) but your model only uses one parameter (alpha). Optim tries to optimize two parameters, with the second one completely unused—this confuses the algorithm and leads to garbage results.
3. Non-convex objective: Summing raw errors instead of squared errors creates an objective function that’s not convex, making it hard for optim to find the global minimum.

The Correct Approach #

We’ll refactor the objective function to use standard regression loss, fix the parameter count, and use a better initial value.

1. Define a Valid Loss Function #

For nonlinear regression, least squares (sum of squared errors) is the standard choice. If you want to account for different group sizes (n), use weighted least squares (weight each error by the group’s sample size):

# Standard least squares loss: minimize sum of squared errors
ls_loss <- function(alpha, data) {
  pred <- with(data, 1 / (1 + alpha * (1 - real)/real))
  sum((data$rate - pred)^2)
}

# Weighted least squares: accounts for group sample size n
weighted_ls_loss <- function(alpha, data) {
  pred <- with(data, 1 / (1 + alpha * (1 - real)/real))
  sum(data$n * (data$rate - pred)^2)
}

2. Pick a Better Initial Value #

Use the average of your uniroot.all results as the starting point—this is already close to the true solution, helping optim converge faster:

init_alpha <- mean(scaler)

3. Run optim() Correctly #

Now pass a single initial value and the valid loss function:

# Standard least squares optimization
opt_ls <- optim(par = init_alpha, fn = ls_loss, data = opt, method = "BFGS", hessian = TRUE)

# Weighted least squares optimization
opt_weighted <- optim(par = init_alpha, fn = weighted_ls_loss, data = opt, method = "BFGS", hessian = TRUE)

4. Check the Results #

View the estimated alpha and compare predicted vs. actual rates:

# Standard LS alpha estimate
cat("Standard LS alpha:", opt_ls$par, "\n")

# Predicted rates from standard LS
pred_ls <- with(opt, 1 / (1 + opt_ls$par * (1 - real)/real))
cbind(Actual_Rate = opt$rate, Predicted_Rate = pred_ls)

# Weighted LS alpha estimate
cat("\nWeighted LS alpha:", opt_weighted$par, "\n")

# Predicted rates from weighted LS
pred_weighted <- with(opt, 1 / (1 + opt_weighted$par * (1 - real)/real))
cbind(Actual_Rate = opt$rate, Predicted_Weighted = pred_weighted)

Alternative: Use nls() for Nonlinear Least Squares #

For this kind of problem, R’s nls() function is more intuitive—it handles the loss function and optimization internally:

# Fit the nonlinear model directly
nls_model <- nls(rate ~ 1/(1 + alpha*(1 - real)/real), 
                 data = opt, 
                 start = list(alpha = init_alpha))

# View model summary
summary(nls_model)

# Get predicted values
predict(nls_model)

Verify Against Your Original Validation #

You can check if the weighted average of predicted rates matches the weighted average of actual rates:

# Weighted average of actual rates
actual_weighted_mean <- sum(opt$rate * opt$n) / sum(opt$n)

# Weighted average of predicted rates (standard LS)
pred_weighted_mean_ls <- sum(pred_ls * opt$n) / sum(opt$n)

# Compare
round(actual_weighted_mean, 3) == round(pred_weighted_mean_ls, 3)

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