在嵌套函数中传递参数和默认参数
Passing arguments and default arguments in nested functions
我有一个函数 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
}
请注意,鉴于我对积分下限和上限的默认设置,我希望当函数为 运行 且仅指定参数时,它将 return 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 或参数有关。
这里似乎有两个问题,都与传递参数有关:第一个传递的参数太多,第二个传递的参数太少。
首先,在您的 x_pdf
定义中,您使用了一个接受单个参数的匿名函数 (function(X){dfun(x=X)}
),但您还尝试传递其他参数(params
list) 到带有 do.call
的匿名函数,这将引发错误。该部分应该看起来像这样:
do.call(dfun, c(list(x = X), params))
现在,您已经将 x_pdf
定义为 require 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-02-17。
如果您已经读到这里,您可能还会喜欢 Hadley Wickham 的 Advanced R。
更多额外位
@andrewH 我从你的评论中了解到你可能正在寻找截断分布的方法,例如求高于整个分布均值的部分分布的均值。
要做到这一点,仅将第一时刻的被积函数从平均值向上积分是不够的:您还必须重新缩放被积函数中的 PDF,使其在截断后再次成为正确的 PDF (弥补丢失的概率质量,如果你愿意的话,用 "hand wave-y" 修辞格)。您可以通过将原始 PDF 的积分除以截断后的 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
我有一个函数 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
}
请注意,鉴于我对积分下限和上限的默认设置,我希望当函数为 运行 且仅指定参数时,它将 return 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 或参数有关。
这里似乎有两个问题,都与传递参数有关:第一个传递的参数太多,第二个传递的参数太少。
首先,在您的 x_pdf
定义中,您使用了一个接受单个参数的匿名函数 (function(X){dfun(x=X)}
),但您还尝试传递其他参数(params
list) 到带有 do.call
的匿名函数,这将引发错误。该部分应该看起来像这样:
do.call(dfun, c(list(x = X), params))
现在,您已经将 x_pdf
定义为 require 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-02-17。
如果您已经读到这里,您可能还会喜欢 Hadley Wickham 的 Advanced R。
更多额外位
@andrewH 我从你的评论中了解到你可能正在寻找截断分布的方法,例如求高于整个分布均值的部分分布的均值。
要做到这一点,仅将第一时刻的被积函数从平均值向上积分是不够的:您还必须重新缩放被积函数中的 PDF,使其在截断后再次成为正确的 PDF (弥补丢失的概率质量,如果你愿意的话,用 "hand wave-y" 修辞格)。您可以通过将原始 PDF 的积分除以截断后的 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