CODEWITHSUNDEEP
rnumerical-integration
M.C. Park
M.C. Park

Reputation: 305

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!

Upvotes: 3

Views: 1199

Answers (2)

ThomasIsCoding
ThomasIsCoding

Reputation: 102241

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

Upvotes: 2

John Coleman
John Coleman

Reputation: 52008

The error that you reported is due to the fact that the function in integrate must be vectorized and integrate itself isn't vectorized.

From the help (?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.

Thus one "fix" is to replace your definition of f1 by:

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

But when I run the resulting code I always get:

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

A solution might be to assemble a large number of quantiles and then smooth it out and use that rather than your Q, although the error here strikes me as odd.

Upvotes: 1

Related Questions