scipy.optimize curve_fit錯誤

我試圖將一個整體的功能使用scipy.optimize curve_fit我的數據：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功能，但這個沒有任何工作（我得到的類型錯誤：不是一個Python函數）。

輸入文件（sn2006gy）的鏈接位於：https://www.dropbox.com/s/ng4lknaja4igwhm/sn2006gy?dl=0。我運行：> python test.py sn2006gy和我得到上述錯誤。任何幫助將不勝感激。

來源

2016-09-30 manolis07gr

示例代碼缺少幾個變量定義，所以除了不顯示數據外，我無法運行它。該錯誤可能與函數的返回值的不合適形狀有關。有功能打印結果的形狀，看看你得到什麼。 –

我的壞基督徒！我從兩個文件複製代碼，忘記粘貼變量定義，它們在下面。我也無法弄清楚如何上傳輸入文件。天= 86400。年= 3.15436e7 Msun = 1.99e33 C = 2.99792458e10 SB = 5.67051e-5 kms2cms = 1.e5 R15 = 1.e15 TNI = 8.8 TCO = 111.3 ENI = 3.9e10 生態= 6.8e9 TH = 5500。 A0 = 1.e14 E51 = 1.e51 L45 = 1.e45 nmax = 1000000 – manolis07gr

請將這些定義添加到問題中的代碼中，而不是將它們轉儲到註釋中。您可以將輸入文件上傳到外部網站（例如文件提取器）並添加鏈接到您的問題。 –

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.

後，現在它是由你來解決數學結束。

來源

2016-10-03 16:56:12

scipy.optimize curve_fit錯誤

回答

相關問題