首次在R中使用optimize()函数：我的实现是否正确？

问题背景 #

本人非专业领域出身，通过谷歌和ChatGPT编写了R代码，试图替代Excel求解器寻找毕达哥拉斯期望模型中最小MSE对应的最优指数系数c。不确定代码逻辑是否正确，尤其疑惑：

• optimize()函数的参数传递逻辑
• 目标函数中predicted_xPts_perc、se、mean(se)三个值的合理性

提供的示例数据与代码如下：

data <- structure(list(xg_for_per_match = c(2.431, 2.256, 2.399, 1.86777777777778, 
                                    2.56, 1.776, 1.59222222222222, 1.665, 0.907, 1.691, 1.52555555555556, 
                                    0.56, 1.46555555555556, 0.754, 1.955, 1.285, 1.38428571428571, 
                                    1.282, 1.86625, 1.5075), xg_against_per_match = c(0.757, 0.682, 
                                                                                      0.739, 1.40555555555556, 0.849, 1.409, 1.17777777777778, 1.822, 
                                                                                      1.938, 2.489, 1.70666666666667, 3.136, 2.37555555555556, 2.985, 
                                                                                      0.78375, 1.56, 1.43, 0.812, 0.945, 0.94875), actual_points_perc = c(0.866666666666667, 
                                                                                                                                                          0.766666666666667, 0.733333333333333, 0.666666666666667, 0.6, 
                                                                                                                                                          0.6, 0.592592592592593, 0.566666666666667, 0.4, 0.4, 0.37037037037037, 
                                                                                                                                                          0.2, 0.111111111111111, 0, 0.791666666666667, 0.625, 0.619047619047619, 
                                                                                                                                                          0.6, 0.541666666666667, 0.416666666666667), se = c(0.0020194397996264, 
                                                                                                                                                                                                             0.0223794165299103, 0.0323996328399904, 0.000796304748754376, 
                                                                                                                                                                                                             0.0905484326863811, 0.00018817673840017, 0.00288909832711458, 
                                                                                                                                                                                                             0.0124545498241967, 0.0485423096732081, 0.0070889390780061, 0.0054423086667482, 
                                                                                                                                                                                                             0.0285940157162698, 0.027082829359042, 0.00359736139900119, 0.00488179512531288, 
                                                                                                                                                                                                             0.048737642866935, 0.0183025693952632, 0.0129244422758347, 0.0646460987510766, 
                                                                                                                                                                                                             0.0897732500472823)), row.names = c(NA, -20L), class = c("tbl_df", 
                                                                                                                                                                                                                                                                      "tbl", "data.frame"))

# 目标函数
objective_function <- function(c) {
  predicted_xPts_perc <- data$xg_for_per_match^c / (data$xg_for_per_match^c + data$xg_against_per_match^c)
  se <- (predicted_xPts_perc - data$actual_points_perc)^2
  mean(se)
}

# 优化目标
result <- optimize(f = objective_function, interval = c(0, 5), lower = 1, upper = 3)
optimal_c <- result$minimum
print(optimal_c)

代码验证与解答 #

1. 参数传递逻辑 #

• optimize()会自动在指定的区间内（代码中lower=1和upper=3优先级高于interval=c(0,5)，实际优化范围是1到3），不断传入不同的c值到objective_function，计算对应MSE值，最终找到使MSE最小的c。
• 目标函数定义为接受单个参数c，完全符合optimize()对一维优化目标函数的要求，逻辑正确。

2. 核心值合理性验证 #

• predicted_xPts_perc: 这是毕达哥拉斯期望的标准变形，用预期进球（xg）替代实际进球计算预期积分占比，完全契合模型核心逻辑——球队积分占比由攻防预期表现的比值决定。
• se: 计算每个样本的预测值与实际值的平方误差，这是求解MSE的必要步骤，平方误差能放大较大偏差的影响，符合最小化MSE的优化目标。
• mean(se): 对所有样本的平方误差取均值得到整体均方误差，这正是我们要最小化的目标函数值，逻辑正确。

3. 代码小优化建议 #

• 同时设置interval和lower/upper会造成冗余，optimize()优先使用lower和upper，建议删除interval参数，修改为：

result <- optimize(f = objective_function, lower = 1, upper = 3)

• 数据框中自带的se列是冗余的，目标函数会重新计算平方误差，建议删除该列，避免歧义。

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