Проблема интеграции в R, когда я использую функцию «интегрировать»

Question

Я пытаюсь вычислить тип индекса Джини, используя сгенерированный набор данных. Но у меня возникла проблема с последней функцией интегрирования. Если я пытаюсь интегрировать функцию с именем 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

Как упомянуто @ John Coleman , integrate должна иметь векторизованную функцию и соответствующую опцию subdivisions для выполнения интегральной задачи. Даже если вы уже предоставили векторизованную функцию для интеграла, иногда сложно правильно установить subdivisions в integrate(...,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

Answer 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, хотя эта ошибка мне кажется нечетный.