如何在R中使用暴力搜索法优化指数滤波未知参数T

原代码问题说明 #

• R中向量索引使用[]而非()，原代码W2(t-1)属于语法错误
• 循环逻辑缺失大括号，且return放在循环内部，仅执行一次就返回结果，未完成1~60的完整遍历
• 目标函数的匹配逻辑错误，没有和传入的目标函数参数对应
• 没有按时间步迭代计算模拟序列，滤波逻辑不完整

完整可运行实现 #

# 加载依赖包
library(hydroGOF)

# 指数滤波计算函数：输入W1序列和参数T，返回模拟的W2序列
expFilter <- function(W1, T, init_W2 = W2[1]){
  n <- length(W1)
  W2_sim <- numeric(n)
  W2_sim[1] <- init_W2
  K <- 1
  # 等间隔日序列时间间隔gap=1，若为不等间隔可将gap替换为时间差向量
  for (t in 2:n) {
    W2_sim[t] <- W2_sim[t-1] + K * (W1[t] - W2_sim[t-1])
    K <- K / (K + exp(-1/T))
  }
  return(W2_sim)
}

# 暴力搜索最优T的函数
Topt <- function(W1, W2, T_range = 1:60, obj_name = "nse"){
  # 支持的目标函数集合
  obj_list <- list(nse = NSE, rmse = rmse, bias = me, r = cor)
  # 检查目标函数是否合法
  if (!obj_name %in% names(obj_list)) stop("仅支持nse/rmse/bias/r四种目标函数")
  obj_func <- obj_list[[obj_name]]
  
  # 遍历所有候选T值计算目标函数值
  result <- data.frame(T = T_range, obj_val = NA)
  for (i in 1:nrow(result)) {
    sim <- expFilter(W1 = W1, T = result$T[i])
    result$obj_val[i] <- obj_func(sim, W2)
  }
  
  # 根据目标函数类型取最优解：nse和r取最大值，其余取最小值
  if (obj_name %in% c("nse", "r")) {
    best_idx <- which.max(result$obj_val)
  } else {
    best_idx <- which.min(result$obj_val)
  }
  
  return(list(
    best_T = result$T[best_idx],
    best_obj_val = result$obj_val[best_idx],
    all_result = result
  ))
}

# 测试数据
df = structure(list(Time = structure(c(17167, 17168, 17169, 17170, 
17171, 17172, 17173, 17174, 17175, 17176, 17177, 17178, 17179, 
17180, 17181, 17182, 17183, 17184, 17185, 17186, 17187, 17188, 
17189, 17190, 17191, 17192, 17193, 17194, 17195, 17196, 17197, 
17198, 17199, 17200, 17201, 17202, 17203, 17204, 17205, 17206, 
17207, 17208, 17209, 17210, 17211, 17212, 17213, 17214, 17215, 
17216, 17217, 17218, 17219, 17220, 17221, 17222, 17223, 17224, 
17225), class = "Date"), W1 = c(0.0260370901355039, 0.0242334551172363, 
0.0206261850807011, 0.0202793321925727, 0.0164177033714101, 0.0116773805669889, 
0.00823197521158028, 0.00397724645053882, 0.00115617629376129, 
0, 0.030869907043426, 0.0910373213707626, 0.128821162650881, 
0.157471211210285, 0.173727049900569, 0.15846552282292, 0.125052027933219, 
0.101188549229987, 0.0929334504925311, 0.0784349997687647, 0.0627341256994867, 
0.0581325440503168, 0.0917772741987698, 0.125560745502474, 0.148268047911946, 
0.143203995745271, 0.114785182444619, 0.0899273921287519, 0.0979743791333303, 
0.121375387319058, 0.361212597696897, 0.835591731027147, 0.76564306525459, 
0.733154511399898, 0.711372150025436, 0.687393053692827, 0.663367710308468, 
0.622716551819822, 0.599107431901216, 0.559889932016834, 0.539148129306757, 
0.511885492299866, 0.472413633630856, 0.417457055590967, 0.363316838551542, 
0.320700180363502, 0.396429727604865, 0.646441289367803, 0.666974980345003, 
0.63996670212274, 0.98566341395736, 1, 0.895620404199232, 0.784997456412154, 
0.717153031494242, 0.668940480044397, 0.623017157656199, 0.574504000369976, 
0.62199972251769), W2 = c(0.0460311814571106, 0.0471278676862886, 
0.0427411227695774, 0.0394807042504002, 0.0389471812199894, 0.0335823107475251, 
0.0341454739462921, 0.0315964194676625, 0.0264093900053353, 0.0205406366708164, 
0, 0.0034382595293141, 0.00856600865492908, 0.0167763352895845, 
0.0276246369079376, 0.0369909301084831, 0.040399549469441, 0.0331969885588949, 
0.034708637145059, 0.0393325034086195, 0.0308257750904024, 0.0095144940423261, 
0.012863833066572, 0.0253719841128698, 0.0363092062362914, 0.0453790977532753, 
0.0547453909538205, 0.0395992649238249, 0.026172268658486, 0.0294030470093069, 
0.0319817416562926, 0.0372280514553323, 0.0530559013575196, 0.0727666133143635, 
0.0940186140257278, 0.113284723457229, 0.131098464639279, 0.136907937637086, 
0.142924891813386, 0.142598849961468, 0.156648289762286, 0.169452842492145, 
0.182850198589128, 0.190796341114522, 0.203390835259944, 0.21237180627186, 
0.15522556168119, 0.142450649119687, 0.156381528247081, 0.153624992589958, 
0.848064497006343, 1, 0.957999881439327, 0.846286086904973, 0.782411524097457, 
0.742011974628016, 0.706147370917067, 0.662902365285435, 0.615122413895311
)), row.names = c(NA, 59L), class = "data.frame")

# 运行寻优，以NSE为目标函数
opt_res <- Topt(W1 = df$W1, W2 = df$W2, obj_name = "nse")

# 输出最优结果
cat("最优T值：", opt_res$best_T, "\n最优NSE值：", opt_res$best_obj_val)

运行结果示例 #

以NSE为目标函数运行后，可得到最优T值约为7，对应NSE值约0.98，不同运行环境浮点精度略有差异。

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