当我使用“集成”功能时，R 中的集成问题答案

【问题标题】：Integration problem in R when I use the function "integrate"当我使用“集成”功能时，R 中的集成问题
【发布时间】：2020-04-17 11:13:35
【问题描述】：

我正在尝试使用生成的数据集计算一种基尼指数。但是，我在最后一个集成功能中遇到了问题。如果我尝试集成名为 f1 的函数， R 说

Error in integrate(Q, 0, p) : length(upper) == 1 is not TRUE

我的代码是

# set up parameters b>a>1 and the number of observations n
n <- 1000
a <- 2
b <- 4

# generate x and y
# where x follows beta distribution
# y = 10x+3
x <- rbeta(n,a,b)
y <- 10*x+3

# the starting point of the integration having problem
Q <- function(q) {
  quantile(y,q)
}

# integrate the function Q from 0 to p
G <- function(p) {
  integrate(Q,0,p)
}

# compute a function
L <- function(p) {
  numer <- G(p)$value
  dino <- G(1)$value
  numer/dino
}

# the part having problem
d <- 3
f1 <- function(p) {
  ((1-p)^(d-2))*L(p)
}
integrate(f1,0,1) # In this integration, the aforementioned error appears

我认为，重复集成可能会产生问题，但我不知道确切的问题是什么。请帮帮我！

【问题讨论】：

f1(0.1) （例如）会引发错误，因此尝试集成本身需要调试的函数似乎有点奇怪。 L(0.1) 也会引发错误。为什么不在尝试定义使用它的其他函数之前调试一个函数？
这是一个很奇怪的问题。 Q 似乎在整合时表现得病态，尽管 curve(Q(x),0,1) 显示了一个相当温和的图表。我认为默认的quantile 函数给出了顺序统计的分段线性插值，所以这种行为看起来很奇怪，甚至可能是其中一个 R 函数中的错误。
玩弄这个，如果我注释掉你的最后一行，源它，并评估integrate(Q,0,1) 大约三分之一的时间我收到错误消息Error in integrate(Q, 0, 1) : maximum number of subdivisions reached 但其他时候我源它该错误不会发生。

标签： r integrate

【解决方案1】：

正如@John Coleman 所提到的，integrate 需要具有矢量化函数和适当的subdivisions 选项才能完成积分任务。即使您已经为积分提供了矢量化函数，在integrate(...,subdivisions = ) 中正确设置subdivisions 有时也很棘手。

为了解决您的问题，我推荐 integral 来自包 pracma，其中您仍然是积分的矢量化函数（请参阅我对函数 G 和 L 所做的），但无需设置细分手动，即，

library(pracma)

# set up parameters b>a>1 and the number of observations n
n <- 1000
a <- 2
b <- 4

# generate x and y
# where x follows beta distribution
# y = 10x+3
x <- rbeta(n,a,b)
y <- 10*x+3

# the starting point of the integration having problem
Q <- function(q) {
  quantile(y,q)
}

# integrate the function Q from 0 to p
G <- function(p) {
  integral(Q,0,p)
}

# compute a function
L <- function(p) {
  numer <- Vectorize(G)(p)
  dino <- G(1)
  numer/dino
}

# the part having problem
d <- 3
f1 <- function(p) {
  ((1-p)^(d-2))*L(p)
}

res <- integral(f1,0,1)

你会得到

> res
[1] 0.1283569

【讨论】：

不错。知道为什么 integrate 与 OP 的 Q 函数有问题（至少有时）？
@JohnColeman integrate 在我认为的某些情况下高度依赖 subdivisions，即使 OP 已经对函数进行了矢量化，我也不确定这是否是原因
非常感谢！但是，我遇到了一个新问题。这个操作需要很长时间，实际上我必须重复这个计算大约 50 次。就你而言，没有耗时的问题吗？？
@minchul 是的，对我来说同样的问题。它是由Vectorize 引起的，但我不知道如何规避它。也许您可以在数学上简化积分过程并在可能的情况下重试
@ThomasIsCoding 好的，再次感谢您！你的解决方案太好了！

【解决方案2】：

您报告的错误是由于integrate中的函数必须矢量化而integrate本身没有矢量化。

来自帮助 (?integrate)：

f 必须接受输入向量并生成函数向量在这些点进行评估。 Vectorize 函数可能有助于将 f 转换为这种形式。

因此，一个“修复”是将您对 f1 的定义替换为：

f1 <- Vectorize(function(p) {
  ((1-p)^(d-2))*L(p)
})

但是当我运行生成的代码时，我总是得到：

Error in integrate(Q, 0, p) : maximum number of subdivisions reached

一种解决方案可能是组装大量分位数，然后将其平滑并使用它而不是您的Q，尽管这里的错误让我觉得很奇怪。

【讨论】：