Python: интеграл и функция, вложенные в другой интеграл с использованием scipy quad

Question

Мне удалось написать несколько строк кода, используя scipy.integrate.quad для моего класса случайных процессов

У меня есть марковская функция перехода для стандартного броуновского движения

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

Но я хочу вычислить следующее, что я собираюсь написать в коде, который не будет работать. Я пишу это так, чтобы мы могли понять проблему без использования латекса.

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

Вы, вероятно, заметили, что проблема заключается в двумерной вещи, происходящей во внутреннем интеграле. Я сделал следующее, но не уверен в этом:

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]

Я тогда получаю

 0.10806767286289147

Я люблю Python и его лямбда-функции, но, похоже, не уверен в этом. о чем ты думаешь? Спасибо за ваше время.

Answer 1

Для типа интеграла, который вы хотите выполнить, двумерные интегралы, SciPy имеет специальные процедуры.

Преимущество состоит в том, что эти процедуры легче обрабатывают сложные границы (например, если они зависят от другой координаты).

Я переписал ваш пример как:

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])

который печатает тот же результат. Я полагаю, что приведенный выше код легче читать, так как подынтегральная функция написана явно как функция двух переменных.