Integrate() in R Gives Terribly Wrong Answer

R integrate: returns wrong solution (is using wrong quadrature points?)

you could shift the integration variable to centre the peak,

wrapper <- function(x, func_mean_v, ...)
   integrandFunc_F(x+func_mean_v, func_mean_v=func_mean_v, ...)

integrate(wrapper, rel.tol = 1e-8, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
          func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
          func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)
# 1 with absolute error < 1.3e-09

Precision issue of integrate function in R

You are advised to use Inf,-Inf. Also it seems, interestingly, that integrate() is converting -Inf to Inf when used as upper limit. :

> integrate(p, -Inf, -1*.Machine$double.xmax)
0 with absolute error < 0
> integrate(p, -Inf, -2*.Machine$double.xmax)
1.772454 with absolute error < 4.3e-06

This didn't result in what we would expect. Lets try to split the integral in two, first:

> integrate(p, 0, Inf)
0.8862269 with absolute error < 2.2e-06
> integrate(p, 0, 1e100)
0 with absolute error < 0
> integrate(p, 0, 1e2)
0.8862269 with absolute error < 9.9e-11

This seems perfectly consistent with the advice of using Inf,-Inf, but notice what happens, when switching signs:

> integrate(p, 0, -Inf)
0.8862269 with absolute error < 2.2e-06
> integrate(p, 0, -1e100)
0 with absolute error < 0
> integrate(p, 0, -1e2)
-0.8862269 with absolute error < 9.9e-11

At last you can always change tolerance, subdivisions, but this will not help you in this case.

EDIT: As @Ben pointed out, it's a minor bug because integrate checks whether the upper limit is finite and not if it's Inf (or -Inf).

Integration problem in R when I use the function integrate

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

Integrate function in R returns error

Your default tolerance (.Machine$double.eps^0.25) is too large, so we change it:

> tail <- function(x) {585.1961*x^-2.592484}
> integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^0.5 )
0.0004727982 with absolute error < 1.5e-09

and this is approximately the same as:

> integrate(tail, 1000, Inf)$val-integrate(tail, 1000, 5000)$val
[1] 0.0004726847

Of course you can just set your tolerance to .Machine$double.eps, but this comes at a cost of time:

> library(microbenchmark)
> a<-function(x) for(i in 1:50) integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^x )
> microbenchmark(a(0.5),a(0.7),a(1))
Unit: milliseconds
   expr      min       lq   median       uq      max neval
 a(0.5) 10.44027 10.97920 11.12981 11.40529 19.70019   100
 a(0.7) 13.02904 13.69813 13.95942 14.89460 23.02422   100
   a(1) 15.14433 15.96499 16.12595 16.38194 26.27847   100

eg. a increase in time around 50 pct.

---

Related Topics