r - Uniroot 函数查找自定义 CDF 的逆

Question

我正在尝试找到自定义 CDF 函数的逆函数

在这样做的过程中，我偶然发现了寻求解决反函数的 uniroot 函数。但是，我遇到了一些错误，非常感谢您对此的帮助。

我的目标是按照 R 的基本函数qnorm qexp等编写一个函数，但使用qcustom.

# Error 1: (Error in f(x) : argument "p" is missing, with no default)

inverse <- function(f, lower=0.01, upper=500000) {
  function(y) {
    uniroot(function(x) f(x) - y, lower=lower, upper=upper)[1]
  }
}

qcustom <- inverse(function(x, p, a, b) {
  x^2+a+b+p
  return(result)
})

qcustom(0.5) ##when i add in further parameters, e.g p = 0.5 it says it's an unused parameter##

标签： r

Answer 1

让我们从头开始：

custom_cdf <- function(x, p, a, b) {
   p*exp(-a*x) + (1-p)*exp(-b*x)
}

这有参数c(p, a, b)，让我们试着找到一个对你说你想要c(a,b)的有意义的配对。range(x) = c(0, 500e3)

plot.new()
curve(custom_cdf(x, p = 0.5, a = 1e-4, b = 1e-5), xlim = c(0, 500e3))

一组合理的值可以是：

a <- 1e-4
b <- 1e-5

让我们尝试多个值p：

first <- TRUE
plot.new()
for (p in seq.default(0.1, 1, length.out = 25)) {
  curve(custom_cdf(x, p = p, a = 1e-4, b = 1e-5), xlim = c(0, 500e3), add = !first)
  first <- FALSE
}
rm('p')

对于 CDF 的倒数，我们可以：

xx <- seq.default(0, 500e3, length.out = 2500)
# xx <- seq.default(0, 500e3, by = 10)
# xx <- seq.default(0, 500e3, by = 1)
yy <- custom_cdf(xx, p = 0.5, a = a, b = b)

plot.new()
plot(yy, xx, main = "inverse of CDF = Quantile Function",type = 'l')

即我们可以反转绘图。所以现在，这是反函数的图片。

让我们用它uniroot来获得一个反值。

uniroot(
  function(x) custom_cdf(x, p = 0.5, a = a, b = b) - 0.2,
  interval = c(0, 100e50),
  extendInt = "yes"
) -> found_20pct_point
#' 
#' Let's plot that point:
#' 
points(custom_cdf(found_20pct_point[[1]], p = 0.5, a = a, b = b), found_20pct_point[[1]])

好的，所以现在我们了解了一切，我们可以做你想做的事情：

custom_quantile <- function(p, a, b) {
  particular_cdf <- function(x) custom_cdf(x, p, a, b)
  function(prob) {
    uniroot(
      function(x) particular_cdf(x) - prob,
      interval = c(0, 100e50),
      extendInt = "yes"
    )$root
  }
}
custom_quantile_specific <- 
  custom_quantile(p = 0.5, a = a, b = b)
custom_quantile_specific(0.2)

这很有效，因为它产生了接近期望值的东西。

但是 R 中的常规分位数函数确实取决于分布的参数，所以这就足够了：

custom_quantile <- function(prob, p, a, b) {
  particular_cdf <- function(x) custom_cdf(x, p, a, b)
  uniroot(
    function(x) particular_cdf(x) - prob,
    interval = c(0, 100e50),
    extendInt = "yes"
  )$root
}
custom_quantile(0.2, p = 0.5, a = a, b = b)