BUG: results of "scipy.dblquad" are confused!

[Leisures2]

Describe your issue.

I want to calculate integration of the pdf of bivariate_normal distribution (should be 1).
But I got an extremely small result using "dblquad" and don't know the reason.

Reproducing Code Example

import numpy as np
import scipy as sp
from scipy.integrate import dblquad  
Cov_BossGal1 = np.array(
    [[2.860520e-02, -4.939281e-02],
     [-4.939281e-02, 5.307187e-01],
     ])
x_BossGal1 = np.array([10.23, 25.00])
D_try = sp.stats.multivariate_normal(x_BossGal1,Cov_BossGal1)
print(D_try.cdf(np.array([np.infty,np.infty]))) # cdf = 1(There is no problem)
print(dblquad(lambda y, x: D_try.pdf(np.array([x, y])), -np.infty, np.infty, -np.infty, np.infty)) #this one is so small!

Error message

scipy.integrate.dblquad

SciPy/NumPy/Python version information

scipy.integrate

[tupui]

Hi @Leisures2, thank you for reporting. The issue is that the PDF mode is very localised, so if you integrate numerically over a very large domain it is very easy to miss a spike. Just changing the bounds solves it:

dblquad(lambda x, y: D_try.pdf(np.array([x, y])), 0, 100, 0, 100)
# (0.9969079885378376, 1.4899985581153049e-08)

If that's not an option, then you need to play with the tolerances, e.g.

dblquad(lambda x, y: D_try.pdf(np.array([x, y])), -np.infty, np.infty, -np.infty, np.infty, epsabs=1e-50)
# (0.999999999999939, 5.4931195671689236e-09)

As you discretise the problem further, it takes more time.

I am closing now as I don't think there is any issue here.