慢 scipy 双积分

Question

我正在尝试获取作为集成结果的函数 expected_W 或 H：

地点：

• theta 是一个包含两个元素的向量：theta_0 和 theta_1
• f(beta | theta) 是 beta 的正态密度，均值为 theta_0，方差为 theta_1
• q(epsilon) 是 epsilon 的正态密度，均值为 0，方差为 sigma_epsilon（默认设置为 1）。
• w(p, theta, eps, beta) 是我作为输入的函数，因此我无法准确预测它的外观。它可能是非线性的，但不是特别讨厌。

这是我实现问题的方式。我确信我制作的包装函数很混乱，所以我也很乐意在这方面获得任何帮助。

from __future__ import division
from scipy import integrate
from scipy.stats import norm
import math
import numpy as np


def exp_w(w_B, sigma_eps = 1, **kwargs):
    '''
    Integrates the w_B function

    Input:
    + w_B : the function to be integrated. 
    + sigma_eps : variance of the epsilon term. Set to 1 by default
    '''

    #The integrand function gives everything under the integral:
    # w(B(p, \theta, \epsilon, \beta)) f(\beta | \theta ) q(\epsilon)
    def integrand(eps, beta, p, theta_0, theta_1, sigma_eps=sigma_eps):
        q_e = norm.pdf(eps, loc=0, scale=math.sqrt(sigma_eps))
        f_beta = norm.pdf(beta, loc=theta_0, scale=math.sqrt(theta_1))

        return w_B(p = p, 
                   theta_0 = theta_0, theta_1 = theta_1,
                   eps = eps, beta=beta)* q_e *f_beta

    #limits of integration. Using limited support for now.
    eps_inf = lambda beta : -10 # otherwise: -np.inf
    eps_sup = lambda beta : 10  # otherwise: np.inf
    beta_inf = -10
    beta_sup = 10

    def integrated_f(p, theta_0, theta_1):
        return integrate.dblquad(integrand, beta_inf, beta_sup,
            eps_inf, eps_sup,
            args = (p, theta_0, theta_1))
    # this integrated_f is the H referenced at the top of the question
    return integrated_f

我用一个简单的 w 函数测试了这个函数，我知道它的解析解（通常情况并非如此）。

def test_exp_w():
    def w_B(p, theta_0, theta_1, eps, beta):
        return 3*(p*eps + p*(theta_0 + theta_1) - beta)

    # Function that I get
    integrated = exp_w(w_B, sigma_eps = 1)

    # Function that I should get
    def exp_result(p, theta_0, theta_1):
        return 3*p*(theta_0 + theta_1) - 3*theta_0

    args = np.random.rand(3)
    d_args = {'p' : args[0], 'theta_0' : args[1], 'theta_1' : args[2]}

    if not (np.allclose(
    integrated(**d_args)[0], exp_result(**d_args)) ):
        raise Exception("Integration procedure isn't working!")    

因此，我的实现似乎正在运行，但就我的目的而言它非常慢。我需要将这个过程重复数万或数十万次（这是价值函数迭代中的一个步骤。如果人们认为它相关，我可以提供更多信息）。

对于scipy 0.14.0 版和numpy 1.8.1 版，计算这个积分需要15 秒。

有人对如何解决这个问题有任何建议吗？ 首先， tt 可能有助于获得有界的积分域，但我还没有弄清楚如何做到这一点，或者 SciPy 中的高斯正交是否能很好地处理它（它使用 Gauss-Hermite 吗？） .

感谢您的宝贵时间。

---- 编辑：添加分析时间 -----

%lprun 结果表明大部分时间都花在 _distn_infraestructure.py:1529(pdf) 和 _continuous_distns.py:97(_norm_pdf) 每个都有高达 83244 个电话号码。

Answer 1

如果函数不是很糟糕的话，集成函数所花费的时间听起来很长。

我建议你做的第一件事是分析时间花在哪里。是花在dblquad 还是其他地方？在集成过程中向w_B 拨打了多少电话？如果时间花在dblquad，调用次数非常多，可以在集成中使用更宽松的容差吗？

似乎高斯乘法实际上使您能够极大地限制积分，因为高斯的大部分能量都在一个非常小的区域内。您可能想尝试计算合理的更严格的界限。您已经将区域限制为 -10..10； -100..100、-10..10 和 -1..1 之间是否有任何显着的性能变化？

---

如果你知道你的功能比较流畅，那么有米老鼠版的集成：

• 确定两个轴的合理上限和下限（通过高斯）
• 计算合理的网格密度（例如，每个方向 100 个点）
• 计算每个点的 w_B（如果可能需要 w_B 的矢量化版本，这会更快）
• 总结起来

这是非常低技术但也非常快。它是否为您提供足够好的外部迭代结果是一个有趣的问题。它只是可能。