Python - 如何将曲线拟合到包含数值计算积分的函数？

Question

我有以下代码：

import numpy as np
import scipy.integrate as spi
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as mh

def GUFunction(z, Omega_Lambda):
    integral = spi.quad(lambda zvar: AuxIntegrandum(zvar, Omega_Lambda), 0.0, z)[0]
    DL = (1+z) * c/H0 * integral *1000000
    return (5*(mh.log(DL,10)-1))

def AuxIntegrandum(z, Omega_Lambda):
    Omega_m = 1 - Omega_Lambda

    return 1 / mh.sqrt(Omega_m*(1+z)**3 + Omega_Lambda)

def DataFit(filename):
    print curve_fit(GUFunction, ComputeData(filename)[0], ComputeData(filename)[1])

DataFit("data.dat")

data.dat 在第一列有 z 值，在第二列有 GUF(z) 值。

执行此代码时，编译器告诉我将数组与值（+inf 或 -inf）进行比较是不明确的。
我认为这是指集成边界，它看起来是否要集成到无穷大。出于某种原因，它显然会将数据文件中的所有 z 值放入集成边界。
是否有一些我不知道的技巧可以让您将曲线拟合到数值积分函数？

这是确切的错误：

Traceback (most recent call last):
  File "plot.py", line 83, in <module>
    DataFit("data.dat")
  File "plot.py", line 67, in DataFit
    print curve_fit(GUFunction, ComputeData(filename)[0], ComputeData(filename)[1])
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 736, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 377, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 454, in func_wrapped
    return func(xdata, *params) - ydata
  File "plot.py", line 57, in GUFunction
    integral = spi.quad(lambda zvar: AuxIntegrandum(zvar, Omega_Lambda), 0.0, z)[0]
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 323, in quad
    points)
  File "/home/joshua/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 372, in _quad
    if (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()

Answer 1

简答：curve_fit 尝试在 xdata 数组上计算目标函数，但 quad 不能接受向量参数。您需要通过例如定义目标函数对输入数组的列表推导。

让我们做一个最小可重现的例子：

In [33]: xdata = np.linspace(0, 3, 11)

In [34]: ydata = xdata**3

In [35]: def integr(x):
    ...:     return quad(lambda t: t**2, 0, x)[0]
    ...: 

In [36]: def func(x, a):
    ...:     return integr(x) * a
    ...: 

In [37]: curve_fit(func, xdata, ydata)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-37-4660c65f85a2> in <module>()
----> 1 curve_fit(func, xdata, ydata)

 [... removed for clarity ...]

~/virtualenvs/py35/lib/python3.5/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()

这正是您看到的错误。好的，错误来自quad，它试图评估func(xdata, a)，归结为integr(xdata)，这不起作用。 （我是怎么发现的？我将import pdb; pdf.set_trace() 放在func 函数中，然后在调试器中四处寻找）。

然后，让我们让目标函数处理数组参数：

In [38]: def func2(x, a):
    ...:     return np.asarray([integr(xx) for xx in x]) * a
    ...: 

In [39]: curve_fit(func2, xdata, ydata)
Out[39]: (array([ 3.]), array([[  3.44663413e-32]]))