Выполнить двойной интеграл по массиву - PullRequest
Новая
       36

Выполнить двойной интеграл по массиву

0 голосов
/ 15 ноября 2018

Я выполняю одномерный интеграл по массиву, как описано здесь и здесь . Как указано в этих ответах, я не могу использовать scipy.integrate.quad для векторизации интегралов по массиву, потому что он использует адаптивный алгоритм , поэтому я использую numpy.trapz ниже (я также мог бы использовать scipy.integrate.simps или scipy.integrate.romb) 

import numpy as np

# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)

# 1D function to integrate
def distFunc(x, c1=1., c2=.1):
    B1 = ((data1 - (1. / x)) / data2)**2
    B2 = ((x - c1) / c2)**2
    f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
    return f

# Values in x to evaluate the integral.
x = np.linspace(.1, 10, 100).reshape(-1, 1)

# Integral in x for each of the Ndata values defined above.
int_exp = np.trapz(distFunc(x), x, axis=0)

Это прекрасно работает для одного измерения, но теперь я хотел бы выполнить двойной интеграл, заменив константу c2 на переменную: 

# 2D function to integrate
def distFunc(x, y, c1=1.):
    B1 = ((data1 - (1. / x)) / data2)**2
    B2 = ((x - c1) / y)**2
    f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / y
    return f

Для того, что я мог собрать, единственная доступная функция - scipy.integrate.dblquad bu, что означает, что я больше не мог применять интеграл ко всему массиву за один проход, и мне пришлось бы использовать for петля, которая значительно медленнее.

Есть ли какое-нибудь решение для этого? Я открыт практически для всего, если это разумно с точки зрения производительности (я включаю этот двойной интеграл в MCMC, и его нужно оценивать миллионы раз)

ADD

Это моя попытка с 1-мерным интегралом с использованием scipy.integrate.quad внутри цикла for (то есть: одно значение данных в массиве за раз). Этот процесс более чем в 50 раз медленнее, чем использование np.trapz по всему массиву. 

import numpy as np
from scipy.integrate import quad

# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)

# Function to integrate
def distFunc2(x, data1_i, data2_i, c1=1., c2=.1):
    B1 = ((data1_i - (1. / x)) / data2_i)**2
    B2 = ((x - c1) / c2)**2
    f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
    return f

s = t.time()
int_exp = np.zeros(Ndata)
for i in range(Ndata):
    int_exp[i] = quad(distFunc2, .1, 10., args=(data1[i], data2[i]))[0]
print(t.time() - s)

ADD 2

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

# Define some random data
Ndata = 10
data1 = np.random.uniform(.1, 10., Ndata)
data2 = np.random.uniform(.1, .2, Ndata)

c1 = .1
print(integ_dblquad(c1, data1, data2))
print(integ_trapz(c1, data1, data2))

def integ_dblquad(c1, data1, data2):
    def distFunc(y, x, d1_i, d2_i, c1):
        B1 = ((d1_i - (1. / x)) / d2_i)**2
        B2 = ((x - c1) / y)**2
        return (np.exp(-.5 * B1) / d2_i) * np.exp(-.5 * B2) / y

    int_exp = np.zeros(data1.size)
    for i in range(data1.size):
        int_exp[i] = dblquad(
            distFunc, .1, 10., lambda x: 0, lambda x: 5.,
            args=(data1[i], data2[i], c1))[0]

    return np.sum(np.log(int_exp))

def integ_trapz(c1, data1, data2):
    def distFunc2d(x, y):
        B1 = ((data1 - (1. / x)) / data2)**2
        B2 = ((x - c1) / y)**2
        return (np.exp(-.5 * B1) / data2) * np.exp(-.5 * B2) / y

    # Values in x to evaluate the integral.
    x = np.linspace(.1, 10, 1000)
    y = np.linspace(.1, 5., 1000)
    # Integral in x for each of the Ndata values defined above.
    int_exp2d = np.trapz(np.trapz(distFunc2d(x[:, np.newaxis], y[:, np.newaxis, np.newaxis]), y, axis=0), x, axis=0)

    return np.sum(np.log(int_exp2d))

1 Ответ

0 голосов
/ 15 ноября 2018

Если я правильно понял вашу проблему, вы можете просто позвонить trapz дважды: 

import numpy as np

# Define some random data
Ndata = 500
data1 = np.random.uniform(.1, .8, Ndata)
data2 = np.random.uniform(.01, .2, Ndata)

# 1D function to integrate
def distFunc(x, c1=1., c2=.1):
    B1 = ((data1 - (1. / x)) / data2)**2
    B2 = ((x - c1) / c2)**2
    f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / c2
    return f

def distFunc2d(x, y, c1=1.):
    B1 = ((data1 - (1. / x)) / data2)**2
    B2 = ((x - c1) / y)**2
    f = np.exp(-.5 * B1) * np.exp(-.5 * B2) / y
    return f

# Values in x to evaluate the integral.
x = np.linspace(.1, 10, 100)
y = np.linspace(.1, 10, 100)

# Integral in x for each of the Ndata values defined above.
int_exp = np.trapz(distFunc(x[:,np.newaxis]), x, axis=0)
int_exp2d = np.trapz(np.trapz(distFunc2d(x[:,np.newaxis],y[:,np.newaxis,np.newaxis]), y, axis=0), x, axis=0)
...