用数值积分的三角函数求解curve_fit问题

Question

试图将模型拟合到观测数据中。代码使用0.5到1.0范围内的数据作为自变量，具有scipy curve_fit和数值积分。要积分的函数还包括一个未知参数，然后使用三角函数sinh(integrand)对被积函数进行求值。

应用curve_fit后，我得到一个错误消息"循环的ufunc不支持参数0类型的函数没有可调用的sinh方法"。我在Python 3上遇到死胡同了吗?希望不会。

这个求值代码是

#O_m, Hu are unknown parameters to be estimated with model, data
def integr(x,O_m):
return intg.quad(lambda x: 1/x(np.sqrt((O_m/x) + (1-O_m))) , x, 1, args=(0.02))[0]
O_m = 0.02 #Guess for value of O_m, which shall lie between 0.01 and 1.0
def funcX(x,O_m):
result = np.asarray([integr(xx,O_m) for xx in x]) * np.sqrt(abs(1-O_m))
return result
litsped=299793 #the constant speed of light in a vacuum (m/s)
def funcY(x,Hu,O_m):
return (litsped/(x * Hu * np.sqrt(abs(1-O_m))))*np.sinh(funcX)
init_guess = [65,0.02]
bnds=([50,0.001],[80,1.0])
params, pcov = curve_fit(funcY, xdata, ydata, p0 = init_guess, bounds = bnds, sigma = error, absolute_sigma = True)
ans_Hu, ans_O_m = params
perr = np.sqrt(np.diag(pcov))

---

完整的代码如下-到目前为止，我已经得到了这个curve_fit。

import numpy as np
import csv
import matplotlib.pylab as plt
from scipy.optimize import curve_fit
from scipy import integrate as intg
with open("Riess_1998_D_L.csv",'r') as i:  #SNe Ia data file
rawdata = list(csv.reader(i,delimiter=",")) #make a data list
exmdata = np.array(rawdata[1:],dtype=float) #convert to data array
xdata = exmdata[:,1]
ydata = exmdata[:,2]
error = exmdata[:,3]
#plot of imported data
plt.title("Observed SNe Ia Data")
plt.figure(1,dpi=120)
plt.xlabel("Expansion factor")
plt.ylabel("Distance (Mpc)")
plt.plot(xdata,ydata,label = "Observed SNe Ia data")
plt.xlim(0.5,1)
plt.ylim(0.0,9000)
plt.xscale("linear")
plt.yscale("linear")
plt.errorbar(xdata, ydata, yerr=error, fmt='.k', capsize = 4)
# O_m and Hu are the unknown parameters which shall be estimated using the model and observational data
def integr(x,O_m):
return intg.quad(lambda x: 1/x(np.sqrt((O_m/x) + (1-O_m))) , x, 1, args=(0.02))[0]
O_m = 0.02 # Guess for value of O_m, which are between 0.01 and 1.0
def funcX(x,O_m):
result = np.asarray([integr(xx,O_m) for xx in x])* np.sqrt(abs(1-O_m))
return result
litsped=299793 #the constant speed of light in a vacuum (m/s)
def funcY(x,Hu,O_m):
return (litsped/(x*Hu*np.sqrt(abs(1-O_m))))*np.sinh(funcX)
init_guess = [65,0.02]
bnds=([50,0.001],[80,1.0])
params, pcov = curve_fit(funcY, xdata, ydata, p0 = init_guess, bounds = bnds, sigma = error, absolute_sigma = True)
ans_b, ans_c = params
perr = np.sqrt(np.diag(pcov))
TotalInt = intg.trapz(ydata,xdata) #Compute numerical integral to check data import
print("The total area is: ", TotalInt)

Answer 1

更多的信息将是有用的，例如，你的xdata/ydata是什么?你能把你的代码重写成一个最小可复制的例子吗?

注:你可以将stackoverflow上的内容格式化为代码，在代码前后写' ' '，以提高可读性;)

Answer 2

您的问题与拟合程序无关。对我来说，理解这些代码有点困难。如果我理解正确的话，我建议这样做:

import numpy as np
import csv
import matplotlib.pylab as plt
from scipy.optimize import curve_fit
from scipy import integrate as intg
exmdata = np.array(np.random.random((10,4)),dtype=float) #convert to data array
xdata = exmdata[:,1]
ydata = exmdata[:,2]
error = exmdata[:,3]
#plot of imported data
plt.title("Observed SNe Ia Data")
plt.figure(1,dpi=120)
plt.xlabel("Expansion factor")
plt.ylabel("Distance (Mpc)")
plt.plot(xdata,ydata,label = "Observed SNe Ia data")
plt.xlim(0.5,1)
plt.ylim(0.0,9000)
plt.xscale("linear")
plt.yscale("linear")
plt.errorbar(xdata, ydata, yerr=error, fmt='.k', capsize = 4)
# O_m and Hu are the unknown parameters which shall be estimated using the model and observational data
def integr(x,O_m):
return 5*x+3*O_m#Some analytical form
O_m = 0.02 # Guess for value of O_m, which are between 0.01 and 1.0
def funcX(x,O_m):
result = integr(x,O_m)* np.sqrt(abs(1-O_m))
return result
litsped=299793 #the constant speed of light in a vacuum (m/s)
def funcY(x,Hu,O_m):
return (litsped/(x*Hu*np.sqrt(abs(1-O_m))))*np.sinh(funcX(x,O_m))
init_guess = np.array([65,0.02])
bnds=([50,0.001],[80,1.0])

params, pcov = curve_fit(funcY, xdata, ydata, p0 = init_guess, bounds = bnds, sigma = error, absolute_sigma = True)

您仍然需要在intr中放入积分的分析形式，并将我的随机数组替换为您的CSV文件数据。你之前提到的错误实际上是由于你传递了整个函数，而不是函数在某个点的值。请首先尝试实现这些步骤，并确保您可以独立调用您的三个函数而不会出现错误。如果您立即处理整个程序，那么查找bug是相当困难的。试着让各个部分先工作;)。如果你在完成这些修改后仍然需要帮助，比如实际的试穿过程，请再问我;)。