scipy.optimize curve_fit 拟合以上限为自变量的积分函数时出错

Question

我正在尝试使用 scipy.optimize curve_fit 将积分函数拟合到我的数据中，如下所示：

import numpy as np
from scipy.integrate import quad, nquad, odeint
from scipy.special import gammainc, gamma
import math
import os
import sys
import contextlib
from scipy.optimize import curve_fit
import string

#Global model constants, all units are CGI
day=86400.     #seconds in a day
year=3.15436e7 #seconds in a year
Msun=1.99e33   #solar mass in grams
c=2.99792458e10#speed of light
sb=5.67051e-5  #Stefan-Boltzmann constant
kms2cms=1.e5   #km/s in cm/s
r15=1.e15      #radii in units of 10^15cm
tni=8.8        #Ni-56 decay time-scale
tco=111.3      #Co-56 decay time-scale
eni=3.9e10     #Ni-56 decay specific energy generation rate
eco=6.8e9      #Co-56 decay specific energy generation rate
TH=5500.       #Hydrogen ionization temperature in K
A0=1.e14       #Radioactive decay model gamma-ray leakage parameter 
E51=1.e51      #Energy in units of 1 F.O.E. (10^51 erg)
L45=1.e45      #Luminosity in units of 10^45 erg/s
nmax=1000000   #Grid resolution for fallback accretion model 

def rad_decay_dep(t,td,r0,vej):
    return ((r0*r15/(vej*td*day*kms2cms))+t/td)*np.exp((t/td)**2  
   (2.*r0*r15*t/(vej*kms2cms*(td**2)*day)))

# Integrants for Ni-56 and Co-56 decay energy depositions

def rad_decay_int1(t,td,r0,vej):
    return rad_decay_dep(t,td,r0,vej)*np.exp(-t/tni)

def rad_decay_int2(t,td,r0,vej):
    return rad_decay_dep(t,td,r0,vej)*np.exp(-t/tco)

# Final radioactive decay luminosity integral function

def Lum_rad(x,Mni,td,r0,vej,A):
    return (2.*Mni*Msun/td)*np.exp(-((x/td)**2+(2.*r0*r15*x
    /(vej*kms2cms*(td**2)*day))))* \
    ((eni-eco)*quad(rad_decay_int1,0,x,args=(td,r0,vej))[0] +  
    eco*quad(rad_decay_int2,0,x,args=(td,r0,vej))[0])* \
    (1.-np.exp(-A*A0/(x*day)**2))

xdata, ydata = np.loadtxt(sys.argv[1], usecols=(0,1), unpack=True)

popt, pcov = curve_fit(Lum_rad, xdata, ydata)

我收到以下错误：

Traceback (most recent call last):
  File "SLSNFit.py", line 11, in <module>
    popt, pcov = curve_fit(LCmods.Lum_rad, xdata, ydata)
  File "/usr/local/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 676, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 377, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/usr/local/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/usr/local/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 455, in func_wrapped
    return func(xdata, *params) - ydata
  File "/Users/macbro/Desktop/LCmods.py", line 86, in Lum_rad
    ((eni-eco)*quad(rad_decay_int1,0,x,args=(td,r0,vej))[0] + eco*quad(rad_decay_int2,0,x,args=(td,r0,vej))[0])* \
  File "/usr/local/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 315, in quad
    points)
  File "/usr/local/lib/python2.7/site-packages/scipy/integrate/quadpack.py", line 364, 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()

我尝试使用 np.vectorize 对函数进行矢量化，但这也不起作用（我得到 TypeError: is not a Python function）。

输入文件 (sn2006gy) 的 link 在这里：https://www.dropbox.com/s/ng4lknaja4igwhm/sn2006gy?dl=0。我运行使用： >python test.py sn2006gy 并且我得到了上述错误。任何帮助将不胜感激。

Answer 1

scipy.integral.quad 无法将数组处理为积分限制。您可以定义一个包装函数，它接受非标量限制，例如

def integral(func, a, b, args):
    ret = np.asarray([quad(func, a, q, args)[0] for q in b])
    return ret

并像这样修改 Lum_rad：

def Lum_rad(x,Mni,td,r0,vej,A):
    return (2.*Mni*Msun/td)*np.exp(-((x/td)**2+(2.*r0*r15*x
    /(vej*kms2cms*(td**2)*day))))* \
    ((eni-eco)*integral(rad_decay_int1,0,x,args=(td,r0,vej)) +  
    eco*integral(rad_decay_int2,0,x,args=(td,r0,vej)))* \
    (1.-np.exp(-A*A0/(x*day)**2))

程序在一段时间后结束说

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1200.

现在轮到你来解决数学问题了。

scipy.optimize curve_fit 拟合以上限为自变量的积分函数时出错

scipy.optimize curve_fit error while fitting an integral function with upper limit as independent variable

curve-fitting

scipy