通过集成创建一个新的 lambda 函数

Question

我在 Python 中使用lambda 函数定义了一个现有函数。该函数表示概率分布的 PDF。我想构建另一个代表 CDF 的 lambda 函数单线。

我不希望使用带有 def 关键字的单独函数定义。

以下是我一直在使用的代码部分：

import numpy as np
import scipy.integrate as integrate

#define range
dx=0.01
X  = np.arange(0,13,dx)

#define a piecewise function for the spline
ul = 1.0
f_pdf = lambda x: np.piecewise(x, [x < ul, x >= ul], [x[x<ul],0])
f_cdf = lambda x: integrate.quad(f_pdf,0,x)

#print the function evaluations
print(f_pdf(X))
print(f_cdf(X))

请注意，我最近发现，对于分段定义，我需要在 x<ul 的情况下限制返回数组的范围，例如 x[x<ul]，以便它能够正确处理不同大小的数组。

我从最后一个命令得到的错误包括：

    ---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-ff6550e6092c> in <module>()
     13 #print the function evaluations
     14 print(f_pdf(X))
---> 15 print(f_cdf(X))

<ipython-input-4-ff6550e6092c> in <lambda>(x)
      9 ul = 1.0
     10 f_pdf = lambda x: np.piecewise(x, [x < ul, x >= ul], [x[x<ul],0])
---> 11 f_cdf = lambda x: integrate.quad(f_pdf,0,x)
     12 
     13 #print the function evaluations

/path/anaconda3/lib/python3.6/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    321     if (weight is None):
    322         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323                        points)
    324     else:
    325         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

/path/anaconda3/lib/python3.6/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    370 def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):
    371     infbounds = 0
--> 372     if (b != Inf and a != -Inf):
    373         pass   # standard integration
    374     elif (b == Inf and a != -Inf):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

当我尝试让它返回与集成端点数组对应的值时。我认为函数integrate.quad不适合将数组作为集成端点，但是有没有合适的替代方案？

我正在使用 Python 3.6.1。

Answer 1

piecewise的第三个参数应该是

可调用对象、f(x,*args,**kw) 或标量列表

而x[x < ul] 两者都不是。正如所写，即使是单个评估 f_cdf(3) 也因此失败。

如果在那里使用可调用对象，

f_pdf = lambda x: np.piecewise(x, [x < ul, x >= ul], [lambda x: x, 0])

评估f_cdf(3) 有效，返回(0.49999999999999994, 5.551115123125782e-16)。您可能只想要这个元组的第一部分，而不是第二部分，这是对集成误差的估计。所以用

f_cdf = lambda x: integrate.quad(f_pdf, 0, x)[0]

由于quad 不接受数组作为其第三个参数（它必须是浮点数），因此尝试将 X 传递给f_cdf 将失败。解决这个问题的一个懒惰方法是在函数上拍vectorize：

f_cdf = np.vectorize(lambda x: integrate.quad(f_pdf,0,x)[0])

将使使用f_cdf(X) 成为可能。然而，这很慢。

与可变上限积分

一种非常快速但不是很准确的方法是使用cumtrapz 方法：

print(integrate.cumtrapz(f_pdf(X), X))

比f_cdf(X) 快得多，但输出有明显的集成错误。

要获得准确的结果，可以编写一个 for 循环，在从 X[i-1] 到 X[i] 的非常短的时间间隔内与 quad 进行积分并添加结果。然而，这并不容易放在一行lambda 中。

精度问题

您的被积函数 f_cdf 在 ul 处是不连续的，这会损害积分的精度。为提高精度，在积分范围内提供points=[ul]参数。像这样：

integrate.quad(f_pdf, 0, x)[0] if x < ul else integrate.quad(f_pdf, 0, x, points=[ul])[0]

Answer 2

请始终提供您的任何错误的完整回溯，以便我们真正知道它在哪里不起作用。

但是从错误中我已经可以说，当您编写 np.piecewise(x, [x < ul, x >= ul], [x[x<ul],0]) 时，您在阵列上运行的逻辑测试就是问题所在。也就是说，当你写x < ul时，你的意思是x中的所有值都必须是x[i] < ul[i]吗？你的意思是至少有一个值满足这个条件？

该问题的修复基本上由错误给出：在np.all() 或np.any() 函数内的数组之间编写任何逻辑测试，例如np.all(x < ul)，如果您希望x 中的所有值小于ul 中的值在同一位置。