在R中寻找数据生成型非多项式函数的最大根

针对你遇到的函数可能存在1-2个根、uniroot()因区间端点符号相同失败的问题，结合已知的h→∞时函数趋向负无穷的特性，推荐以下两种稳健方案：

方案1：使用rootSolve包批量找根后取最大值 #

rootSolve包的uniroot.all()函数可以直接找出指定区间内的所有根，无需手动处理区间符号问题，之后只需取最大的根即可：

步骤： #

1. 安装并加载包
2. 调用uniroot.all()获取所有根
3. 提取最大值作为目标根

# 安装包（仅首次需要）
install.packages("rootSolve")
library(rootSolve)

# 你的目标函数
h.root <- function(x, h) {
  n <- length(x)
  f2 <- function(y) {
    h2 <- h*h
    (sum(((x-y)^2/h2 - 1) * dnorm((x-y)/h)) / (n*h*h2))^2
  }
  f2.int <- integrate(Vectorize(f2), -Inf, Inf)$value
  (0.5 / sqrt(pi) / f2.int / n)^0.2 - h
}

# 示例使用
set.seed(36)
x <- rnorm(50)

# 查找区间内所有根（区间设置足够覆盖可能的根范围）
all_roots <- uniroot.all(h.root, interval = c(0.01, 10), x = x)
# 提取最大根
max_root <- max(all_roots)
print(max_root)

方案2：手动扫描符号变化区间再用uniroot()求解 #

如果不想依赖第三方包，可以先扫描h的取值，找到函数由正变负的区间（这是最大根所在的区间，因为h→∞时函数为负），再用uniroot()精准求解：

# 你的目标函数
h.root <- function(x, h) {
  n <- length(x)
  f2 <- function(y) {
    h2 <- h*h
    (sum(((x-y)^2/h2 - 1) * dnorm((x-y)/h)) / (n*h*h2))^2
  }
  f2.int <- integrate(Vectorize(f2), -Inf, Inf)$value
  (0.5 / sqrt(pi) / f2.int / n)^0.2 - h
}

# 定义扫描函数，寻找函数由正转负的区间
find_max_root_interval <- function(x, start_h = 0.01, step = 0.01, max_h = 10) {
  prev_val <- h.root(x, start_h)
  current_h <- start_h
  
  while (current_h < max_h) {
    current_h <- current_h + step
    curr_val <- h.root(x, current_h)
    # 找到正转负的区间（对应最大根）
    if (prev_val > 0 && curr_val < 0) {
      return(c(current_h - step, current_h))
    }
    prev_val <- curr_val
  }
  stop("未找到函数由正转负的区间，请检查参数或函数特性")
}

# 示例使用
set.seed(36)
x <- rnorm(50)

# 获取目标区间
target_interval <- find_max_root_interval(x)
# 求解最大根
max_root <- uniroot(h.root, target_interval, x = x)$root
print(max_root)

补充说明： #

• 你的临时方案使用extendInt="downX"存在隐患：当函数存在两个根时，该参数可能会扩展区间找到左侧的小根，而非目标最大根。
• 可以优化h.root函数的计算效率，比如避免重复计算h*h，或简化积分逻辑，减少多次调用时的耗时。

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