自定义分布CDF逆函数绘制QQ图遇uniroot端点符号错误求助

错误原因 #

报错f() values at end points not of opposite sign是因为uniroot要求求解区间的端点函数值符号相反，但当前设置的区间无法满足：

• 当p接近1时，CDF在upper=1000处的取值可能已趋近于1，和p的差值符号一致；
• 默认uniroot不会自动扩展区间，初始区间不包含根时直接报错。

修复步骤 #

1. 优化分位数函数 #

修改qf函数，添加区间自动扩展参数，同时处理边界分位数避免无解：

qf<-function(p,b){
  # 限制p在(0.001, 0.999)，避免边界值导致的无解问题
  p <- pmax(pmin(p, 0.999), 0.001)
  # 开启自动扩展区间，返回根的数值而非列表
  uniroot(function(x) cdf(x,b) - p, lower=0, upper=1000, extendInt = "yes")$root
}

2. 简化代码复用 #

复用已定义的cdf函数到pakash，避免重复编写代码：

pakash <- function(q,b) cdf(q,b)

3. 完整可运行代码 #

# 定义Akash分布CDF
cdf<-function(x,b){
  1 - (1+(b*x*(b*x+2)/(b^2+2)))*exp(-b*x)
}

# 修复后的分位数函数
qf<-function(p,b){
  p <- pmax(pmin(p, 0.999), 0.001)
  uniroot(function(x) cdf(x,b) - p, lower=0, upper=1000, extendInt = "yes")$root
}

# 样本数据
DATA1<-c(14, 13, 19, 22, 11,17)

# 定义Akash分布的密度、CDF、分位数函数
dakash <- function(x,b) (b^3/(b^2+2))*(1+x^2)*exp(-b*x)
pakash <- function(q,b) cdf(q,b)
qakash <- function(p,b) qf(p,b)

# 拟合分布
fitgumbel <- fitdist(DATA1, "akash", start=list(b=5))

# 绘制QQ图
qqcomp(list(fitgumbel), legendtext=c("akash"))

额外提示 #

• 若拟合时出现收敛问题，可尝试调整start的初始值（比如从b=1开始）；
• 手动测试分位数函数：qf(0.5, fitgumbel$estimate$b)，确认能正常返回结果。

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