使用函数 "integrate" 时 R 中的集成问题

Question

我正在尝试使用生成的数据集计算一种 Gini 指数。但是，我在最后一个集成功能中遇到了问题。如果我尝试集成名为 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

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

Answer 1

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

来自帮助(?integrate):

f must accept a vector of inputs and produce a vector of function evaluations at those points. The Vectorize function may be helpful to convert f to this form.

因此 "fix" 是将您对 f1 的定义替换为：

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

但是当我运行结果代码时，我总是得到：

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

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

Answer 2

如所述，integrate需要有一个向量化函数和一个适当的subdivisions选项来完成积分任务。即使您已经为积分提供了矢量化函数，有时也很难在 integrate(...,subdivisions = ) 中正确设置 subdivisions。

为了解决您的问题，我推荐包 pracma 中的 integral，其中您仍然是积分的矢量化函数（请参阅我对函数 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" 时 R 中的集成问题

Integration problem in R when I use the function "integrate"

r

integrate