Python: Integral and funcion nested in another integral using scipy quad

Question

I have managed to write a few lines of code using scipy.integrate.quad for my stochastic process class

I have the Markov transition function for standard Brownian motion

import numpy as np
def p(x,t):
    return (1/np.sqrt(2*np.pi*t))*np.exp(-x**2/(2*t))

But I want to compute the following that I am going to write in code that would not work. I write it like this so we can understand the problem without the use of latex.

 from scipy.integrate import quad
 integral = quad(quad(p(y-x),1,np.inf)*p(x,1),1,np.inf) 

You probably noticed that the problem is the bivariate thing going on in the inner integral. I did the following but am unsure of it:

p_xy = lambda y,x: p(y-x,1)
inner = lambda x : quad(p_xy,1,np.inf,args = (x,))[0]
outer = lambda x: inner(x)*p(x,1)
integral = quad(outer,1,np.inf)[0]

I then get

 0.10806767286289147

I love Python and its lambda functions but seem to not be sure about this. What are your thoughts? Thank you for your time.

Answer 1

For the type of integral you wish to perform, bivariate integrals, SciPy has dedicated routines.

The advantage is that these routines handle complex boundaries more easily (were the bounds depend on the other coordinate, for instance).

I rewrote your example as:

import numpy as np
from scipy.integrate import nquad

def p(x,t):
    return (1/np.sqrt(2*np.pi*t))*np.exp(-x**2/(2*t))

def integrand(x, y):
    return p(y-x, 1)*p(x, 1)

integral = nquad(integrand, ((1, np.inf), (1, np.inf)))

print(integral[0])

which prints out the same result. I believe that the code above is easier to read as the integrand is written explicitly as a function of the two variables.