NOAA CO₂测量数据linear fit参数误差问题求助

问题背景 #

处理NOAA的CO₂浓度数据时，线性拟合曲线绘制正常，但拟合参数、参数不确定度和约化卡方值均显示无穷大（inf）。已将负的uncertainty替换为0.1，调整拟合初始猜测值无效，需确认该问题是线性拟合特性还是代码/拟合错误。数据截至2023年9月。

代码实现 #

import numpy as np
import matplotlib.pyplot as plt
from numpy import cos, sqrt, exp, sin, pi, array
from matplotlib import cm
data=np.loadtxt('co2_data.txt', skiprows=2)
from scipy.optimize import curve_fit

def linfit(x, p0, p1):
    return (p0*x)+p1
def quadfit(x, p0, p1, p2):
    return p0*(x**2)+(p1*x)+p2

x=data[:,2]
t=array([i-x[0] for i in x])
y=data[:,3]
yerr=data[:,7]
for i in range(len(yerr)):
    if yerr[i]<0:
       yerr[i]=.1
guesses = (1.5, 310)
(p0, p1), cc = curve_fit(linfit, t, y, p0 = guesses, sigma = yerr)
(up0, up1) = np.sqrt(np.diag(cc))
print(f'p0 = {p0:.4} +/- {up0:.4f} ppm/year')
print(f'p1 = {p1:.4f} +/- {up1:.4f} ppm')
xm = np.linspace(t[0], t[-1], 201)
ym1 = linfit(xm, p0, p1)
yfit1=linfit(t, p0, p1)
n = len(t)
p=2
rcsq = sum(((y - yfit1)/yerr)**2)/(n-p)
print(f'Reduced chisquared for linear fit = {rcsq:.5f}')

报错信息 #

/usr/local/lib/python3.10/dist-packages/scipy/optimize/_minpack_py.py:968: RuntimeWarning: divide by zero encountered in divide
  transform = 1.0 / sigma
<ipython-input-44-c1614db438c8>:9: OptimizeWarning: Covariance of the parameters could not be estimated
  (p0, p1), cc = curve_fit(linfit, t, y, p0 = guesses, sigma = yerr)
<ipython-input-44-c1614db438c8>:18: RuntimeWarning: divide by zero encountered in divide
  rcsq = sum(((y - yfit1)/yerr)**2)/(n-p)
p0 = 1.5 +/- inf ppm/year
p1 = 310.0000 +/- inf ppm
Reduced chisquared for linear fit = inf

原因分析与解决办法 #

• 核心原因：仅处理了负的uncertainty，但未处理值为0的uncertainty。curve_fit计算时会用到1/sigma，约化卡方计算也依赖1/yerr，一旦存在yerr=0的情况，就会触发除零错误，导致结果为inf。
• 修复方案：
  修改yerr的处理逻辑，将所有小于等于0的值替换为合理的极小值（比如0.1）：

  # 用numpy批量处理更高效
  yerr[yerr <= 0] = 0.1
  

  或者保留循环写法：

  for i in range(len(yerr)):
      if yerr[i] <= 0:
          yerr[i] = 0.1
  

• 结论：这不是线性拟合的特性，完全是数据中存在0不确定性值引发的计算错误，和初始猜测值无关。修改后curve_fit可正常计算协方差矩阵，约化卡方也会得到有效数值。

内容的提问来源于stack exchange，提问作者Giau Diep