rnumerical-integration

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


I'm trying to compute a kind of Gini index using a generated dataset. But, I got a problem in the last integrate function. If I try to integrate the function named f1, R says

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

My code is

# 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

I think, the repeated integrate could make a problem but I have no idea what is the exact problem. Please help me!


Solution

  • As mentioned by @John Coleman, integrate needs to have a vectorized function and a proper subdivisions option to fulfill the integral task. Even if you have already provided a vectorized function for integral, it is sometimes tricky to properly set the subdivisions in integrate(...,subdivisions = ).

    To address your problem, I recommend integral from package pracma, where you still a vectorized function for integral (see what I have done to functions G and L), but no need to set subdivisions manually, i.e.,

    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)
    

    then you will get

    > res
    [1] 0.1283569