如何编写通用R函数拟合分位数函数至给定三分位数？

错误根源分析 #

1. 分位数函数调用方式错误：qfun(1/3, par)无法适配通用分位数函数（如qgamma需要shape/rate这类命名参数），直接传入向量par会导致参数解析失败。
2. optim参数传递混乱：内部函数gapfn重复定义外部已有的参数，加上optim的参数匹配规则，导致t2参数传递丢失。
3. 闭包特性未利用：内部函数可直接访问外部函数的t1/t2/qfun变量，无需重复传参，避免冗余和匹配错误。

修正后的通用函数代码 #

tertile_fit <- function(t1, t2, qfun, par_start) {
  # 目标函数：通过闭包直接使用外部变量，无需重复传参
  gapfn <- function(par) {
    # 用do.call将par转为参数列表，适配任意分位数函数的参数格式
    q1 <- do.call(qfun, args = c(list(p = 1/3), as.list(par)))
    q2 <- do.call(qfun, args = c(list(p = 2/3), as.list(par)))
    # 计算平方误差和作为优化目标
    (t1 - q1)^2 + (t2 - q2)^2
  }
  
  # 调用optim执行优化，仅需传入初始参数和目标函数
  result <- optim(par = par_start, fn = gapfn)
  return(result)
}

测试验证 #

以gamma分布为例，和原特定函数结果对比：

# 测试通用函数
fit_result <- tertile_fit(t1 = 0.5, t2 = 2, qfun = qgamma, par_start = c(1, 10))
print(fit_result)

# 和原特定函数对比验证一致性
t1 <- 0.5; t2 <- 2
gapfn_old <- function(par, t1, t2) {
  a <- t1 - qgamma(1/3, shape = par[1], rate = par[2])
  b <- t2 - qgamma(2/3, shape = par[1], rate = par[2])
  return(a^2 + b^2)
}
old_result <- optim(par = c(1,10), gapfn_old, t1, t2)
print(all.equal(fit_result$par, old_result$par)) # 输出TRUE，说明结果一致

扩展说明 #

• 适配其他分布：比如正态分布qnorm，只需传入对应初始参数：

  tertile_fit(t1 = -0.43, t2 = 0.43, qfun = qnorm, par_start = c(0,1))
  

• 参数约束：若分位数函数有参数范围限制（如gamma的参数需大于0），可在optim中添加约束：

  optim(par = par_start, fn = gapfn, method = "L-BFGS-B", lower = rep(1e-6, length(par_start)))
  

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