Integrate() Gives Totally Wrong Number

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).

Error in integrate : evaluation of function gave a result of wrong length

You can try:

integrate(Vectorize(f),0,1)$value

see the manual of integrate: f should an R function taking a numeric first argument and returning a numeric vector of the same length. Vectorize will make f a such function that returns the same length output as input.

Scipy.integrate gives odd results; are there best practices?

As with many numerical methods, it can be very sensitive to asymptotes, zeros, etc. The only choice is to keep giving it 'hints' if it will accept them.

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

Unexpected behaviour of scipy.integrate

It is not a bug, it has to do with the numerical precision of the integration, and the fact that you are integrating a function that is (almost) 0 in most of the interval.
From the docs:

Be aware that pulse shapes and other sharp features as compared to the
size of the integration interval may not be integrated correctly using
this method.

Based on your output, the function is using only two (last=2) intervals, evaluating values for rlist=(2.06145321e-045, 0.00000000e+000, ..) on each (see the docs for more detail on the output)

You can add points to the interval to force the routine to use points closer to the left limit.

a = quad(lambda x: np.exp(-x), 0, 1e9, points=np.logspace(-10,3,10))
print(a)
(0.9999999999999997, 2.247900608926337e-09)

Adding to the explanation (thanks to @norok2): Note that points is a sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). In this case, I'm not using it to point out discontinuities, but rather to force quad to perform more integration steps near the left boundary, using a log-spaced interval since a I have an exponential function (this is of course arbitrary and for this function, since I know its shape).