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

Question

尝试将模型拟合到观测数据。该代码使用 0.5 到 1.0 范围内的数据作为自变量，并使用 scipy curve_fit 和数值积分。要积分的函数还包含未知参数，然后使用三角函数 sinh(被积函数) 对被积函数进行评估。

应用 curve_fit 后，我收到一条错误消息“ufunc 的循环不支持没有可调用 sinh 方法的函数类型的参数 0”。我对 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)

您仍然需要在 intgr 中放入积分的分析形式，并用您的 CSV 文件数据替换我的随机数组。您之前提到的错误确实是由于您传递了整个函数而不是在某个时刻评估的函数。请首先尝试实现这些步骤，并确保您可以独立调用这三个函数而不会出现错误。如果你立即解决整个程序，那么寻找错误是相当困难的。首先尝试让各个部分正常工作；）。如果您在实施这些更改后仍然需要帮助，例如实际的验配程序，请再次询问我;)。