Uniroot：QF 的端点不是相反符号（虽然它们有）

Question

我正在尝试手动编写一个分位数函数，为了说明我的观点，我提出了下一个 CDF（实际上与 pnorm(x,mean=1,标准差=0)):

cdf = function(x){
  integrate(function(x) {dnorm(x, mean=1)},lower = -Inf, upper=x)$value 
}

然后我在分位数函数的代码中使用 root 函数，并测试它生成下一条错误消息：

qf = function(q){
  uniroot(function(x)
    {cdf(x)-q},
    interval=c(-Inf, Inf))$root
}

qf(0.5)
-----------------------------------------------
Error in uniroot(function(x) { : 
  f() values at end points not of opposite sign

一开始以为是无限大的间隔，所以改成相反符号的大值如：c(-1000, 1000).

我能做些什么来修复它，或者 使用不同的方式编程我的分位数函数？

Answer 1

让我引用integrate的文档：

Like all numerical integration routines, these evaluate the function on a finite set of points. If the function is approximately constant (in particular, zero) over nearly all its range it is possible that the result and error estimate may be seriously wrong.

这就是你问题的核心。看看这个：

integrate(function(x) {dnorm(x, mean=1)},lower = -Inf, upper = -Inf)
#1 with absolute error < 1.6e-05

这显然是错误的。来自 uniroot 的错误就是由此产生的。如果可以避免，则不应对密度函数进行数值积分。

Answer 2

除了 Roland 指出的 integrate 和 lower = -Inf, upper = -Inf 的不良行为外，它还对 space 进行了糟糕的采样，以便在端点非常接近时放弃之前进行整合到 0:

pdf <- function(x, mean = 0, sd = 1) {
  p <- dnorm(x, mean, sd)
  print(c(x = x, pdf = p))
  return(p)
}

cdf <- function(x){
  integrate(pdf, lower = -Inf, upper=x)
}

cdf(-1000)
#>        x1        x2        x3        x4        x5        x6        x7        x8 
#> -1001.000 -1233.065 -1000.004 -1038.299 -1000.026 -1013.800 -1000.072 -1006.738 
#>        x9       x10       x11       x12       x13       x14       x15      pdf1 
#> -1000.148 -1003.832 -1000.261 -1002.366 -1000.423 -1001.525 -1000.656     0.000 
#>      pdf2      pdf3      pdf4      pdf5      pdf6      pdf7      pdf8      pdf9 
#>     0.000     0.000     0.000     0.000     0.000     0.000     0.000     0.000 
#>     pdf10     pdf11     pdf12     pdf13     pdf14     pdf15 
#>     0.000     0.000     0.000     0.000     0.000     0.000
#> 0 with absolute error < 0
cdf(1000)
#>       x1       x2       x3       x4       x5       x6       x7       x8 
#> 999.0000 766.9348 999.9957 961.7012 999.9739 986.2000 999.9275 993.2621 
#>       x9      x10      x11      x12      x13      x14      x15     pdf1 
#> 999.8516 996.1681 999.7390 997.6339 999.5774 998.4754 999.3441   0.0000 
#>     pdf2     pdf3     pdf4     pdf5     pdf6     pdf7     pdf8     pdf9 
#>   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 
#>    pdf10    pdf11    pdf12    pdf13    pdf14    pdf15 
#>   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
#> 0 with absolute error < 0

一个解决方法是在密度远离 0 的值处拆分积分，然后在 uniroot 中的平均值周围设置一个较宽的有限搜索区间：

cdf <- function(x, mean = 0, sd = 1) {
  if (x > mean) {
    integrate(function(x) dnorm(x, mean, sd), lower = -Inf, upper = mean)$value +
      integrate(function(x) dnorm(x, mean, sd), lower = mean, upper = x)$value
  } else {
    integrate(function(x) dnorm(x, mean, sd), lower = -Inf, upper = x)$value
  }
}

qf <- function(q, mean = 0, sd = 1) {
  uniroot(function(x) cdf(x, mean, sd) - q, interval = c(-sd, sd)*100 + mean)$root
}

qf(0.5, 1)
#> [1] 1
qnorm(0.5, 1)
#> [1] 1
qf(0.6, 1)
#> [1] 1.253346
qnorm(0.6, 1)
#> [1] 1.253347

Uniroot：QF 的端点不是相反符号（虽然它们有）

Uniroot: End points not of opposite signs (while they have) for QF

r

statistics

quantile