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]]))