使用含数值积分的sinh函数拟合时curve_fit报错求助

错误原因分析 #

你遇到的错误核心是把函数对象直接传给了np.sinh，同时代码里还有几处语法和逻辑错误：

• funcY中直接使用np.sinh(funcX)，但funcX是函数，不是计算后的数值数组，必须调用funcX(x, O_m)得到结果后再传入双曲正弦函数。
• integr函数里的积分被积函数写错：1/x(np.sqrt(...))是函数调用语法，实际应该是1/(x * np.sqrt(...))；另外intg.quad的args=(0.02)是多余的，因为lambda已经捕获了外部的O_m参数。
• funcY的return语句没有缩进，违反Python语法。
• 最后参数解包时ans_b, ans_c和定义的参数Hu, O_m不对应，应该改为ans_Hu, ans_O_m。

修复后的代码 #

分步修复点 #

1. 修正integr函数的被积函数和积分参数：

def integr(x, O_m):
    # 修复被积函数的运算语法，移除多余的args参数，lambda变量名改为t避免冲突
    return intg.quad(lambda t: 1/(t * np.sqrt((O_m / t) + (1 - O_m))), x, 1)[0]

2. 修正funcY的调用逻辑和缩进：

def funcY(x, Hu, O_m):
    # 先调用funcX计算得到数值数组，再传入np.sinh
    fx = funcX(x, O_m)
    return (litsped / (x * Hu * np.sqrt(abs(1 - O_m)))) * np.sinh(fx)

3. 修正参数解包的变量名：

ans_Hu, ans_O_m = params

完整修复代码 #

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数据文件
    rawdata = list(csv.reader(i,delimiter=",")) #转换为数据列表

exmdata = np.array(rawdata[1:],dtype=float) #转换为数值数组
xdata = exmdata[:,1]
ydata = exmdata[:,2]
error = exmdata[:,3]

# 绘制导入的数据
plt.title("观测SNe Ia数据")
plt.figure(1,dpi=120)
plt.xlabel("膨胀因子")
plt.ylabel("距离 (Mpc)")
plt.plot(xdata,ydata,label = "观测SNe Ia数据")
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和Hu是需要通过模型和观测数据估计的未知参数
def integr(x, O_m):
    return intg.quad(lambda t: 1/(t * np.sqrt((O_m / t) + (1 - O_m))), x, 1)[0]

O_m_guess = 0.02 # O_m的初始猜测值，范围0.01到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 # 真空中的光速 (m/s)
def funcY(x, Hu, O_m):
    fx = funcX(x, O_m)
    return (litsped / (x * Hu * np.sqrt(abs(1 - O_m)))) * np.sinh(fx)

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

TotalInt = intg.trapz(ydata,xdata) # 计算数值积分验证数据导入

print("总面积为: ", TotalInt)
print("拟合得到的Hu: ", ans_Hu, "误差: ", perr[0])
print("拟合得到的O_m: ", ans_O_m, "误差: ", perr[1])

补充说明 #

• 把integr里的lambda变量名改为t是为了避免和外部传入的积分下限x混淆，减少不必要的bug。
• 如果拟合过程中出现收敛问题，可以尝试调整初始猜测值init_guess，放宽边界范围bnds，或者给curve_fit添加maxfev参数增加迭代次数。

内容的提问来源于stack exchange，提问作者Mike Smith