0

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.

Yes
  • 3
  • 3
  • Is there a question here? It seems you're just asking for opinions. – B. Mehta Jul 29 '18 at 21:19
  • You aren't sure about the answer or about the most efficient way? To address the former, try using some known benchmark analytical integral. – Sheldore Jul 29 '18 at 21:22
  • I guess that I'm ashamed to be so direct as to ask, "Teacher, is my answer correct?" . I guess that would pretty much sum it up. Do you think the value I got is correct? – Yes Jul 29 '18 at 21:23
  • As far as efficiency goes, I'm happy with this. I actually tried Sympy at first. It seemed to not be able to do this. – Yes Jul 29 '18 at 21:25

1 Answers1

0

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.

Pierre de Buyl
  • 7,074
  • 2
  • 16
  • 22
  • Thank you so much for this valuable learning experience. I like this explanation and additional example a lot! – Yes Jul 30 '18 at 17:13