BUG: dblquad and args kwarg

[RichardAFrazin]

Describe your issue.

There is a bug with scipy.dblquad that seems to be involved with using the args kwarg. The example below is the integral of a 2D Normal distribution, which should be unity. I am using the args kwarg to pass the coords of the center of the Normal. Note that the boundary of the integral is also using the same variable that are passed in args (although I don't know why that might matter).

Reproducing Code Example

< import numpy as np; from scipy.integrate import dblquad
< sig = 1.3; x0=5.; y0=-3.;
< g = lambda x,y,xc,yc : (1/(2*np.pi*sig*sig))*np.exp(- ((x-xc)**2 + (y-yc)**2)/(2*sig*sig) )
< print([g(x0,y0,x0,y0), g(x0+1,y0+1,x0,y0)])
[0.09417452253958304, 0.05211400420209555]
< print([g(x0,y0,x0-5*sig,y0), g(x0+5*sig,y0,x0,y0)])
[3.509557831511321e-07, 3.509557831511321e-07]

< dblquad(g,x0-5*sig,x0+5*sig,y0-5*sig,y0+5*sig,args=(x0,y0))
> (0.015445922091068943, 6.2668582074865746e-09)
this answer is not close enough to unity for my taste.

Error message

no error message, just the wrong answer!

SciPy/NumPy/Python version information

np 1.16.2 scipy 1.2.1

[tupui]

Hi @RichardAFrazin, thank you for reporting. Boundaries and tolerances are very important as it will impact how the domain is discretised and walked through. See #17554.

Unrelated (I think), consider updating SciPy and NumPy. These are very old versions and we pushed in both literally thousands of bug fixes and improvements.

[tupui]

Seeing this is exactly the same case/distribution as the other issue. I did not check your formula, there might as well be an issue there.

Closing for now as I am pretty certain there is no bug. If you still find something I am happy to reopen. Thanks again 😃

[RichardAFrazin]

It looks like the problem in 17554 was due trying to integrate over a narrow spike in a large domain. In my case, the boundaries are tight, extending only to +- 5 standard deviations from the center in each direction. My formula for the normal is correct. I will try updating scipy and numpy

…

On Thu, Jan 12, 2023 at 9:22 AM Pamphile Roy ***@***.***> wrote: Hi @RichardAFrazin <https://github.com/RichardAFrazin>, thank you for reporting. Boundaries and tolerances are very important as it will impact how the domain is discretised and walked through. See #17554 <#17554>. Unrelated (I think), consider updating SciPy and NumPy. These are very old versions and we pushed in both literally thousands of bug fixes and improvements.Message ID: ***@***.***>

[tupui]

Did you try adjusting the tolerance and other parameters of the function?

[RichardAFrazin]

Yes, I did.

…

On Thu, Jan 12, 2023 at 3:08 PM Pamphile Roy ***@***.***> wrote: Did you try adjusting the tolerance and other parameters of the function? — Reply to this email directly, view it on GitHub <#17776 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH47V5QFFPDD3MOIQDCTH33WSBQDRANCNFSM6AAAAAATZJWZSA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[mdhaber]

I modified your code so that:

def g(x,y,xc,yc):
    print(xc, yc)
    return (1/(2*np.pi*sig*sig))*np.exp(- ((x-xc)**2 + (y-yc)**2)/(2*sig*sig) )

It prints lots of (5.0, -3.0). So I don't think args is the problem.
However, I haven't gotten the correct result from the integral yet.

[RichardAFrazin]

I know you can pass the 2nd parameters via the args kwarg, but something screwy is happening.

…

On Thu, Jan 12, 2023 at 3:09 PM Matt Haberland ***@***.***> wrote: I modified your code so that: def g(x,y,xc,yc): print(xc, yc) return (1/(2*np.pi*sig*sig))*np.exp(- ((x-xc)**2 + (y-yc)**2)/(2*sig*sig) ) It prints lots of (5.0, -3.0). So I don't think args is the problem. However, I haven't gotten the correct result from the integral yet. —Message ID: ***@***.***>

[mdhaber]

< dblquad(g,x0-5*sig,x0+5*sig,y0-5*sig,y0+5*sig,args=(x0,y0))

Check the order of your arguments:

Here is code that shows that dblquad is not the problem here:

import numpy as np
from scipy import stats, integrate
sig, x0, y0 = 1.3, 5., -3
cov = stats.Covariance.from_diagonal([sig**2, sig**2])

def g(y, x, xc, yc):
    z = [x, y]
    mean = [xc, yc]
    return stats.multivariate_normal(mean, cov).pdf(z)    

res = integrate.dblquad(g, x0-5*sig, x0+5*sig, y0-5*sig, y0+5*sig, 
                        args=(x0, y0))
print(res)  # (0.9999988533940415, 8.668315648607816e-10)

I find the order (y, x) surprising, too, but it's not something we can change.

[RichardAFrazin]

I did check them carefully. Please point out the error in my original bug report if you think I did not.

…

On Thu, Jan 12, 2023 at 3:14 PM Matt Haberland ***@***.***> wrote: < dblquad(g,x0-5*sig,x0+5*sig,y0-5*sig,y0+5*sig,args=(x0,y0)) Check the order of your arguments: [image: image] <https://user-images.githubusercontent.com/6570539/212171149-4fb3300e-3d07-4acf-ab4a-9c8bc5a8dc67.png>Message ID: ***@***.***>

[mdhaber]

Make it g = lambda y, x,xc,yc : (1/(2*np.pi*sig*sig))*np.exp(- ((x-xc)**2 + (y-yc)**2)/(2*sig*sig) ). Reverse x and y.

[tupui]

I am pretty sure the function is not correct. I tried to plotted it and it's flat.

[mdhaber]

I think the function is fine:

g(y0, x0, x0, y0) - stats.multivariate_normal([x0, y0], cov=np.eye(2)*sig**2).pdf([x0, y0])  # -1.3877787807814457e-17

The main problem is the order that x and y are passed into the function. dblquad expects y, x, but the original code was x, y.

[RichardAFrazin]

If you plotted it with a line going through the center, it's not flat.

…

On Thu, Jan 12, 2023 at 3:18 PM Pamphile Roy ***@***.***> wrote: I am pretty sure the function is not correct. I tried to plotted it and it's flat. — Reply to this email directly, view it on GitHub <#17776 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH47V5XGS643UIO4O3LHMHDWSBRI7ANCNFSM6AAAAAATZJWZSA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[tupui]

https://www.wolframalpha.com/input?i=z+%3D+%281%2F%282*pi*1.3*1.3%29%29*exp%28-+%28%28x-5%29**2+%2B+%28y%2B3%29**2%29%2F%282*1.3*1.3%29+%29

[tupui]

[RichardAFrazin]

Yes, the problem does indeed seem to be y,x instead of x,y. Thank you. I appreciate the help!

…

On Thu, Jan 12, 2023 at 3:21 PM Matt Haberland ***@***.***> wrote: I think the function is fine: g(y0, x0, x0, y0) - stats.multivariate_normal([x0, y0], cov=sig).pdf([x0, y0]) # 0.0 It's just the order that x and y are passed into the function. dblquad <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html> expects y, x, but the original code was x, y. — Reply to this email directly, view it on GitHub <#17776 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH47V5WHXZQHGVAUZDC4QDTWSBRTBANCNFSM6AAAAAATZJWZSA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[tupui]

Actually, if you have no correlation, shouldn't the right term just be exp(.../2)? As 1-p **2 would be 1.

[mdhaber]

I was missing that cov should have been np.eye(2) * sig**2 instead of np.eye(2) * sig before. Not sure what is going wrong with the plot @tupui but the function is OK aside from the argument order.

[RichardAFrazin]

Yes, I agree.

…

On Thu, Jan 12, 2023 at 3:41 PM Matt Haberland ***@***.***> wrote: I was missing that cov should have been np.eye(2) * sig**2 instead of np.eye(2) * sig before. Not sure what is going wrong with the plot @tupui <https://github.com/tupui> but the function is OK aside from the argument order. Message ID: ***@***.***>