使用uncertainties线性拟合时遇ValueError: numpy.object_数据类型错误

问题背景 #

我是Python新手，处理带误差的Excel数据集时遇到以下需求与问题：

• 数据集包含周期、pdot、流量等数据列及对应误差列
• 需要完成：绘制数据、带误差的线性拟合（获取标准差、p值等指标）、基于拟合结果预测缺失参数
• 已实现无误差的拟合代码，但无法获取p值，也未考虑误差传播；尝试用unumpy处理误差后，拟合时抛出类型错误

已实现的无误差拟合代码 #

dist_array1= np.multiply(3.08567758128*10**21,dist_array)
dist_array2 = np.multiply(dist_array1,dist_array1)
e1=np.multiply(4*math.pi,dist_array2)
L_gamma = np.multiply(e1,flux_array)
Gamma_Eff = np.divide(L_gamma,edot_array)
Tau = np.divide(period_array,pdot_array)
constant = 2.94*10**8
t1=np.power(period_array,-5)
t2=np.multiply(t1,pdot_array)
t3=np.power(t2,1/2)
B_LC = np.multiply(constant,t3)
c1=np.multiply(10**15,pdot_array)
c2=np.log(c1)
c3=np.log(period_array)
c4=1-np.multiply(11/7,c3)+np.multiply(4/7,c2)
c5=3.56-c3-c2
Zeta1=1+np.divide(c4,c5)
c6=0.8-np.multiply(2/7,c3)+np.multiply(2/7,c2)
Zeta2=1+np.divide(c6,1.3)
c8=0.6-np.multiply(11/14,c3)+np.multiply(2/7,c2)
Zeta3=1+np.divide(c8,1.3)

拟合部分代码 #

x1 = np.log(period_array)
y1 = np.log(Gamma_Eff)
coef1, V1 = np.polyfit(x1,y1,1, cov=True)
poly1d_fn1 = np.poly1d(coef1) 
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(30,10))
fig.suptitle('Figure 1')
ax1.plot(x1,y1, 'yo', x1, poly1d_fn1(x1), '-k')

x2 = np.log(Tau)
coef2, V2 = np.polyfit(x2,y1,1, cov=True)
poly1d_fn2 = np.poly1d(coef2) 
ax2.plot(x2,y1, 'yo', x2, poly1d_fn2(x2), '-k')

x3= np.log(B_LC)
coef3, V3 = np.polyfit(x3,y1,1, cov=True)
poly1d_fn3 = np.poly1d(coef3) 
ax3.plot(x3,y1, 'yo', x3, poly1d_fn3(x3), '-k')

ax1.set(xlabel='log P (s)', ylabel='log η')
ax2.set(xlabel='log τ (yr)', ylabel='log η')
ax3.set(xlabel='log B_LC (G)', ylabel='log η')

# 获取不确定度
sigma_period_1=np.sqrt(V1[0][0])
sigma_period_2=np.sqrt(V1[1][1])
sigma_Tau_1=np.sqrt(V2[0][0])
sigma_Tau_2=np.sqrt(V2[1][1])
sigma_B_LC_1=np.sqrt(V3[0][0])
sigma_B_LC_2=np.sqrt(V3[1][1])

尝试误差传播后的错误情况 #

为实现误差传播，修改代码用unumpy包装带误差的数组：

period_array= unumpy.uarray(period_array,perioderr_array) # 合并数值与误差
pdot_array=unumpy.uarray(pdot_array,pdoterr_array)
flux_array=unumpy.uarray(flux_array,flux_err_array)
c2=unumpy.log(c1) # 使用unumpy避免log函数报错
c3=unumpy.log(period_array)

执行拟合时：

x1 = unumpy.log(period_array)
y1 = unumpy.log(Gamma_Eff)
coef1, V1 = np.polyfit(x1,y1,1, cov=True)

抛出错误：

ValueError: data type <class 'numpy.object_'> not inexact

问题原因与解决方案 #

核心原因 #

numpy.polyfit仅支持数值型数组（如float64），而unumpy.uarray是object类型数组，每个元素是带误差的自定义对象，numpy无法对其进行数值拟合计算，因此抛出类型错误。

解决方案1：用scipy.optimize.curve_fit处理带误差的拟合 #

curve_fit支持传入y轴误差作为权重，能输出参数协方差，还可通过统计方法计算p值：

import numpy as np
from scipy.optimize import curve_fit
import scipy.stats as stats
import matplotlib.pyplot as plt

# 1. 用unumpy完成误差传播后，提取名义值和标准差
x_nominal = unumpy.nominal_values(x1)
x_err = unumpy.std_devs(x1)
y_nominal = unumpy.nominal_values(y1)
y_err = unumpy.std_devs(y1)

# 2. 定义线性拟合模型
def linear_model(x, slope, intercept):
    return slope * x + intercept

# 3. 带权重的拟合（sigma传入y的误差，absolute_sigma=True表示直接用误差作为权重）
coef, cov = curve_fit(linear_model, x_nominal, y_nominal, sigma=y_err, absolute_sigma=True)
slope, intercept = coef

# 4. 计算p值（基于F检验）
y_pred = linear_model(x_nominal, slope, intercept)
ss_res = np.sum((y_nominal - y_pred)**2)  # 残差平方和
ss_tot = np.sum((y_nominal - np.mean(y_nominal))**2)  # 总平方和
r_squared = 1 - (ss_res / ss_tot)
n = len(x_nominal)
param_count = 2  # 斜率+截距
f_stat = (r_squared / (param_count-1)) / ((1 - r_squared) / (n - param_count))
p_value = stats.f.sf(f_stat, param_count-1, n-param_count)

# 5. 绘制带误差棒的拟合图
plt.errorbar(x_nominal, y_nominal, yerr=y_err, fmt='yo', label='原始数据')
plt.plot(x_nominal, y_pred, '-k', label=f'拟合线: y={slope:.2f}x + {intercept:.2f}\nR²={r_squared:.3f}, p值={p_value:.3e}')
plt.xlabel('log P (s)')
plt.ylabel('log η')
plt.legend()
plt.show()

# 6. 获取参数的标准差（从协方差矩阵提取）
slope_std = np.sqrt(cov[0][0])
intercept_std = np.sqrt(cov[1][1])

解决方案2：仅用unumpy做误差传播，提取数值后用传统方法拟合 #

如果只需要unumpy处理变量的误差传播，得到x1和y1的名义值与标准差后，可选择scipy.stats.linregress（适合简单线性拟合，若需权重仍推荐curve_fit）：

from scipy.stats import linregress

# 提取unumpy数组的名义值和标准差
x_vals = unumpy.nominal_values(x1)
y_vals = unumpy.nominal_values(y1)

# 简单线性拟合（无权重，若需权重用curve_fit）
result = linregress(x_vals, y_vals)
print(f"斜率: {result.slope:.2f}, 截距: {result.intercept:.2f}")
print(f"R²: {result.rvalue**2:.3f}, p值: {result.pvalue:.3e}")

内容的提问来源于stack exchange，提问作者Ertuğrul Karamanlı