Search code examples
rnumerical-integration

Inconsistent results between R2Cuba and integrate()

I have experienced an inconsistent results from integrate() and suave() from package "R2Cuba" when computing a two-fold integration.

where

A similar question related to this topic can be found here, let's set f(x)=6*sin(x) and g(x)=1 and the upper time T=3.

The following code use integrate() and suave() to compute the objective 2-fold integral:

library(R2Cuba)

integrand = function(x){6*sin(x)}

phi = function(x){integrate(integrand,lower=x,upper = 3)[["value"]]^2}

NDIM=1
NCOMP=1

phicuba= function(x){suave(NDIM,NCOMP,integrand,lower=x,upper=3)$value^2}

foldintegral = integrate(Vectorize(phi),lower = 0,upper = 3)

foldintegralcuba = suave(NDIM,NCOMP,phicuba,lower = 0,upper = 3)

The results:

> foldintegral
167.3934 with absolute error < 1.9e-12
> foldintegralcuba
integral: 6.365749 (+-0.0057)
nregions: 8; number of evaluations:  10000; probability:  1

which is not consistent. However, if we only compare phi and phicuba

> phi(2)
[1] 11.85476
> phicuba(2)
Iteration 1:  1000 integrand evaluations so far
[1] 3.44367 +- 0.0429822    chisq 0 (0 df)
Iteration 2:  2000 integrand evaluations so far
[1] 3.4436 +- 0.00866171    chisq 0.00513973 (2 df)
Iteration 3:  3000 integrand evaluations so far
[1] 3.44181 +- 0.00386046   chisq 5.38913 (5 df)
Iteration 4:  4000 integrand evaluations so far
[1] 3.44283 +- 0.00206345   chisq 23.0816 (8 df)
[1] 11.85309

We got a result that can be seen as consistent. Moreover, if we use replace the integrand inside suave()

> suave(NDIM,NCOMP,phi,lower = 0,upper = 3)
Iteration 1:  1000 integrand evaluations so far
[1] 167.426 +- 4.91983      chisq 0 (0 df)
Iteration 2:  2000 integrand evaluations so far
[1] 167.315 +- 0.771248     chisq 0.00208764 (2 df)
Iteration 3:  3000 integrand evaluations so far
[1] 167.357 +- 0.432941     chisq 5.50239 (5 df)
Iteration 4:  4000 integrand evaluations so far
[1] 167.362 +- 0.129012     chisq 5.51661 (8 df)
integral: 167.3621 (+-0.13)
nregions: 4; number of evaluations:  4000; probability:  0.2988006

We still have consistent result.

I wish I can use the suave() since it is much faster when dealing with complicated integral, but why this inconsistent exists?

=========================UPDATE===============================

It seems it is because the algorithm used in suave(), From the Cuba Documentation, suave use globally adaptive subdivision + importance sampling. If change suave() to cuhre(), which only use globally adaptive subdivision, everything should be consistent.

Solution

  • It seems it is because the algorithm used in suave(), From the Cuba Documentation, suave use globally adaptive subdivision + importance sampling. If change suave() to cuhre(), which only use globally adaptive subdivision, everything should be consistent.