Uniroot 函数查找自定义 CDF 的逆答案

【问题标题】：Uniroot function to find the inverse of a custom CDFUniroot 函数查找自定义 CDF 的逆
【发布时间】：2021-10-26 08:48:54
【问题描述】：

我正在尝试找到自定义 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##

【问题讨论】：

看起来像force 问题。目前inverse 函数的uniroot 部分被分配给qcustom 而不是评估它。见adv-r.hadley.nz/function-factories.html。您能否提供一些您期望函数返回的示例？

标签： r

【解决方案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)

【讨论】：