CODEWITHSUNDEEP
pythonmathscipyintegralmismatch
develarist
develarist

Reputation: 1365

Double integral solution using scipy.integrate.nquad doesn't match integrate.dblquad

The first function in the code below uses double integration with scipy.integrate.dblquad to calculate the differential entropy, c*np.log(c), of a copula density function c, which has one dependence parameter, theta, usually positive.

enter image description here

The second function in the code below tries to solve the same problem as above, but using the multiple integral solver scipy.integrate.nquad

from scipy import integrate
import numpy as np

def dblquad_(theta):
    "Double integration"
    c = lambda v, u: ((1+theta)*(u*v)**(-1-theta)) * (u**(-theta)+v**(-theta)-1)**(-1/theta-2)
    return -integrate.dblquad(
        lambda u,v: c(v,u)*np.log(c(v,u)), 
        0, 1, lambda u: 0, lambda u: 1
        )[0]

def nquad_(n,theta):
    "Multiple integration"
    c = lambda *us: ((1+theta)*np.prod(us)**(-1-theta)) * (np.sum(np.power(us,-theta))-1)**(-1/theta-2)
    return -integrate.nquad(
        func   = lambda *us : c(*us)*np.log(c(*us)), 
        ranges = [(0,1) for i in range(n)],
        args   = (theta,) 
        )[0] 

n=2
theta = 1
print(dblquad_(theta))
print(nquad_(n,theta))

The dblquad-based function gives an answer of -0.7127, while the nquad gives -0.5823 and noticeably takes longer. Why are the solutions different even though I have set both to solve an n=2-dimensional problem?

Upvotes: 1

Views: 688

Answers (1)

David M.
David M.

Reputation: 4588

With the values of n and theta you provided, the output of your code is:

-0.1931471805597395
0.17055845832017144,

not -0.7127 and -0.5823.

The first value (-0.1931471805597395) is the correct one (you can check it here for yourself).

The problem in nquad_ lies with the handling of the theta parameter. Thanks @mikuszefski for providing the explanation; I copy it here for clarity:

nquad passes lambda to the function as required. The lambda is programmed in such a way that it accepts arbitrary number of arguments, so it happily takes it and puts it in the list of powers and sums. Hence you do not get, e.g. 1/u**t+1/v**t -1 but 1/u**t+1/v**t + 1/t**t -1. The function call just does not match the intended function use. If you instead would write us[0]**() + us[1]**() - 1 it works.

Here is the modified code:

from scipy import integrate
import numpy as np

def dblquad_(theta):
    "Double integration"
    c = lambda v, u: ((1+theta)*(u*v)**(-1-theta)) * (u**(-theta)+v**(-theta)-1)**(-1/theta-2)
    return -integrate.dblquad(
        lambda u,v: c(v,u)*np.log(c(v,u)),
        0, 1, lambda u: 0, lambda u: 1
        )[0]

def nquad_(n,theta):
    "Multiple integration"
    c = lambda *us: ((1+theta)*np.prod((us[0], us[1]))**(-1-theta)) * (np.sum(np.power((us[0], us[1]),-theta))-1)**(-1/theta-2)
    return -integrate.nquad(
        func   = lambda *us : c(*us)*np.log(c(*us)),
        ranges = [(0,1) for i in range(n)],
        args=(theta,)
        )[0]

n=2
theta = 1
print(dblquad_(theta))
print(nquad_(n,theta))

Output:

-0.1931471805597395
-0.1931471805597395

Upvotes: 2

Related Questions