Python中概率密度函数的更快卷积

Posted 2023-02-23

【中文标题】Python中概率密度函数的更快卷积

【英文标题】：Faster convolution of probability density functions in Python

【发布时间】：2015-05-08 04:50:35

【问题描述】：

假设需要计算一般数量的离散概率密度函数的卷积。对于下面的示例，有四种分布采用指定概率的值 0、1、2：

import numpy as np
pdfs = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1],[0.3,0.7,0.0],[1.0,0.0,0.0]])

卷积可以这样找到：

pdf = pdfs[0]        
for i in range(1,pdfs.shape[0]):
    pdf = np.convolve(pdfs[i], pdf)

看到 0,1,...,8 的概率由下式给出

array([ 0.09 ,  0.327,  0.342,  0.182,  0.052,  0.007,  0.   ,  0.   ,  0.   ])

这部分是我代码中的瓶颈，似乎必须有一些可用的东西来向量化这个操作。有没有人建议让它更快？

或者，您可以使用的解决方案

pdf1 = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1]])
pdf2 = np.array([[0.3,0.7,0.0],[1.0,0.0,0.0]])
convolve(pd1,pd2) 

得到成对卷积

 array([[ 0.18,  0.51,  0.24,  0.07,  0.  ], 
        [ 0.5,  0.4,  0.1,  0. ,  0. ]])

也会有很大帮助。

【问题讨论】：

根据 numpy 文档，np.convolve 的参数只能是一维的。所以我想，这里没有太多要矢量化的东西。但也许它值得使用不同的卷积，比如 scipy 的基于 fft 的卷积？ docs.scipy.org/doc/scipy/reference/generated/…

@SmCaterpillar 我玩了一下，但我对卷积的了解太有限，无法理解那里发生了什么。我理解这里的版本，但我不知道如何为 fft 版本指定权重。

你说的重量是什么意思？我尝试了两种方法，两种卷积都为您的问题提供了相同的结果。但是，fft 的速度要慢得多（由于开销，您的玩具问题太小了，也许当 pdf 本身包含更多值时，您实际上会提高速度）。

@SmCaterpillar 我想您再次将 for 循环用于 scipy 版本并一一进行卷积。我想避免 for 循环并立即将操作应用于所有 pdf 行。

我正在查看这个版本的 convolve 记录docs.scipy.org/doc/scipy/reference/generated/…

【参考方案1】：

您可以使用快速傅立叶变换 (FFT) 有效地计算所有 PDF 的卷积：关键事实是 FFT of the convolution 是各个概率密度函数的 FFT 的乘积。所以转换每个 PDF，将转换后的 PDF 相乘，然后执行逆变换。您需要将每个输入 PDF 用零填充到适当的长度，以避免回绕的影响。

这应该是相当有效的：如果您有m PDF，每个都包含n 条目，那么使用此方法计算卷积的时间应该增长为(m^2)n log(mn)。时间由 FFT 控制，我们有效地计算 m + 1 独立 FFT（m 正向变换和一个逆变换），每个数组的长度不大于 mn。但与往常一样，如果你想要真正的时间，你应该分析一下。

这里有一些代码：

import numpy.fft

def convolve_many(arrays):
    """
    Convolve a list of 1d float arrays together, using FFTs.
    The arrays need not have the same length, but each array should
    have length at least 1.

    """
    result_length = 1 + sum((len(array) - 1) for array in arrays)

    # Copy each array into a 2d array of the appropriate shape.
    rows = numpy.zeros((len(arrays), result_length))
    for i, array in enumerate(arrays):
        rows[i, :len(array)] = array

    # Transform, take the product, and do the inverse transform
    # to get the convolution.
    fft_of_rows = numpy.fft.fft(rows)
    fft_of_convolution = fft_of_rows.prod(axis=0)
    convolution = numpy.fft.ifft(fft_of_convolution)

    # Assuming real inputs, the imaginary part of the output can
    # be ignored.
    return convolution.real

将此应用于您的示例，这就是我得到的：

>>> convolve_many([[0.6, 0.3, 0.1], [0.5, 0.4, 0.1], [0.3, 0.7], [1.0]])
array([ 0.09 ,  0.327,  0.342,  0.182,  0.052,  0.007])

这是基本思想。如果您想对此进行调整，您还可以查看numpy.fft.rfft（及其相反的numpy.fft.irfft），它利用输入是真实的这一事实来生成更紧凑的转换数组。您还可以通过用零填充 rows 数组来获得一些速度，以便总列数对于执行 FFT 是最佳的。此处“最优”的定义取决于 FFT 实现，但例如，2 的幂将是很好的目标。最后，如果所有输入数组的长度相同，则在创建 rows 时可以进行一些明显的简化。但我会将这些潜在的增强功能留给您。

【讨论】：

为什么不使用scipy.signal.fftconvolve() (docs.scipy.org/doc/scipy/reference/generated/…)？

@Dietrich：因为（除非我遗漏了什么）一次只对两个数组进行卷积，并且重复使用它会涉及很多不必要的转换和反转换。

@MarkDickinson 您能否详细说明我们如何将输出（密度概率）与实际结果相匹配？那就是我们如何计算这些概率所属的结果？ result_length 的目的是什么？为什么我们要向每个数组添加一些零，因为我们只填充 rows 数组直到 :len(array)?.

@user2974951 这有点取决于你在做什么。通常，您要对随机变量 X_1、X_2、...、X_n 的 PDF 进行卷积，以获得 X_1 + X_2 + ... + X_n 的 PDF。在这种情况下，每个X_i 需要是离散的，可能的值使用一些间距s 均匀分布，并且该间距需要与所有X_i 匹配。然后结果对应于间距为s 的均匀间隔值，从每个X_i 的最小值之和开始，到每个X_i 的最大值之和。

@user2974951 result_length，顾名思义，就是生成的卷积数组的长度。我们需要将每个输入数组填充到该长度，以便 FFT 兼容。