在嵌套函数中传递参数和默认参数

Question

我有一个函数x_pdf，它应该计算 x*dfun(x|params)，其中 dfun 是概率密度函数，params 是命名参数的列表。它是在另一个函数int_pdf 中定义的，该函数应该在指定范围内集成x_pdf：

int_pdf <- function(lb = 0, ub = Inf, dfun, params){
  x_pdf <- function(X, dfun, params){X * do.call(function(X){dfun(x=X)}, params)}
    out <- integrate(f = x_pdf, lower=lb, upper=ub, subdivisions = 100L)
  out
}

请注意，考虑到我对积分下限和上限的默认设置，我希望在仅指定参数的情况下运行该函数时，它将返回 x 的平均值。

我有第二个函数int_gb2，它是int_pdf 的包装器，旨在将其专门用于第二类的通用 beta 分布。

library(GB2)

int_gb2 <- function(lb = 0, ub = Inf, params){
  int_pdf(lb, ub, dfun = dgb2, params = get("params"))
}

当我按如下方式运行函数时：

GB2_params   <-  list(shape1 = 3.652, scale = 65797, shape2 = 0.3, shape3 = 0.8356)
int_gb2(params = GB2_params)

我明白了：

 Error in do.call(what = function(X) { : 
  argument "params" is missing, with no default

我花了几个小时对此进行了调整，并且已经生成了许多替代错误消息，但总是针对缺少的 x、X 或参数。

Answer 1

这里似乎有两个问题，都与传递参数有关：第一个传递的参数太多，第二个传递的参数太少。

首先，在您的x_pdf 定义中，您使用了一个接受单个参数（function(X){dfun(x=X)}）的匿名函数，但您还尝试将其他参数（params 列表）传递给带有@ 的匿名函数987654326@，这将引发错误。那部分应该看起来像这样：

do.call(dfun, c(list(x = X), params))

现在，您已将x_pdf 定义为需要 3 个参数：X、dfun 和params；但是当你用integrate 调用x_pdf 时，你没有传递dfun 和params 参数，这又会引发错误。你也可以通过 dfun 和 params 来解决这个问题：

integrate(f = x_pdf, lower=lb, upper=ub, subdivisions = 100L, dfun, params)

但也许更简洁的解决方案是从 x_pdf 的定义中删除额外的参数（因为 dfun 和 params 已经在封闭环境中定义），以获得更紧凑的结果：

int_pdf <- function(lb = 0, ub = Inf, dfun, params){
  x_pdf <- function(X) X * do.call(dfun, c(list(x = X), params))
  integrate(f = x_pdf, lower = lb, upper = ub, subdivisions = 100L)
}

有了int_pdf的这个定义，一切都应该如你所愿：

GB2_params <- list(shape1 = 3.652, scale = 65797, shape2 = 0.3, shape3 = 0.8356)
int_gb2(params = GB2_params)
#> Error in integrate(f = x_pdf, lower = lb, upper = ub, subdivisions = 100L):
#>   the integral is probably divergent

哦。示例参数是否缺少 scale 参数中的小数点？

GB2_params$scale <- 6.5797
int_gb2(params = GB2_params)
#> 4.800761 with absolute error < 0.00015

---

额外的位

我们还可以使用一些函数式编程来创建一个函数工厂，以便轻松地创建用于查找除第一个以外的时刻的函数：

moment_finder <- function(n, c = 0) {
  function(f, lb = -Inf, ub = Inf, params = NULL, ...) {
    integrand <- function(x) {
      (x - c) ^ n * do.call(f, c(list(x = x), params))
    }
    integrate(f = integrand, lower = lb, upper = ub, ...)
  }
}

要找到均值，您只需创建一个函数来找到第一个时刻：

find_mean <- moment_finder(1)

find_mean(dnorm, params = list(mean = 2))
#> 2 with absolute error < 1.2e-05
find_mean(dgb2, lb = 0, params = GB2_params)
#> 4.800761 with absolute error < 0.00015

对于方差，您必须找到第二个中心时刻：

find_variance <- function(f, ...) {
  mean <- find_mean(f, ...)$value
  moment_finder(2, c = mean)(f, ...)
}

find_variance(dnorm, params = list(mean = 2, sd = 4))
#> 16 with absolute error < 3.1e-07
find_variance(dgb2, lb = 0, params = GB2_params)
#> 21.67902 with absolute error < 9.2e-05

或者，我们可以进一步概括，并找到期望值 任何转变，而不仅仅是瞬间：

ev_finder <- function(transform = identity) {
  function(f, lb = -Inf, ub = Inf, params = NULL, ...) {
    integrand <- function(x) {
      transform(x) * do.call(f, c(list(x = x), params))
    }
    integrate(f = integrand, lower = lb, upper = ub, ...)
  }
}

现在moment_finder 将是一个特例：

moment_finder <- function(n, c = 0) {
  ev_finder(transform = function(x) (x - c) ^ n)
}

由reprex package (v0.2.0) 于 2018 年 2 月 17 日创建。

如果您读到这里，您可能还会喜欢 Hadley Wickham 的 Advanced R。

---

更多额外位

@andrewH 我从您的评论中了解到您可能正在寻找截断分布的方法，例如找到高于整个分布平均值的分布部分的平均值。

为此，仅将一阶矩的被积函数从平均值向上积分是不够的：您还必须在被积函数中重新缩放 PDF，以使其在截断后再次成为正确的 PDF（弥补对于丢失的概率质量，如果你愿意的话，在“手波-y”修辞格中）。您可以通过在截断的支持上除以原始 PDF 的积分来做到这一点。

下面的代码可以更好地表达我的意思：

library(purrr)
library(GB2)

find_mass <- moment_finder(0)
find_mean <- moment_finder(1)

GB2_params <- list(shape1 = 3.652, scale = 6.5797, shape2 = 0.3, shape3 = 0.8356)
dgb2p <- invoke(partial, GB2_params, ...f = dgb2)  # pre-apply parameters

# Mean value
(mu <- find_mean(dgb2p, lb = 0)$value)
#> [1] 4.800761

# Mean for the truncated distribution below the mean
(lower_mass <- find_mass(dgb2p, lb = 0, ub = mu)$value)
#> [1] 0.6108409
(lower_mean <- find_mean(dgb2p, lb = 0, ub = mu)$value / lower_mass)
#> [1] 2.40446

# Mean for the truncated distribution above the mean
(upper_mass <- find_mass(dgb2p, lb = mu)$value)
#> [1] 0.3891591
(upper_mean <- find_mean(dgb2p, lb = mu)$value / upper_mass)
#> [1] 8.562099

lower_mean * lower_mass + upper_mean * upper_mass
#> [1] 4.800761