Python uncertainties包结合curve_fit高斯拟合误差传播报错排查

我希望使用Python的uncertainties包实现高斯拟合流程中的误差传播。持有一组带测量不确定度的实验测量数据x与y，y与x近似满足高斯函数关系，采用高斯曲线对数据拟合，期望得到各拟合参数的不确定度。
初始实现代码如下：

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

def gaus_fit(x,a,x0,sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2*sigma/abs(sigma)))

def curve_fit_float(x_arr,y_arr,p0):
    popt_X,pcov_X = curve_fit(gaus_fit,x_arr,y_arr,p0)
    return popt_X[2]

arr_x=np.array([    0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,   263.46478873,   380.        ,   614.74647887,
         946.07042254,  1838.18309859,  3226.18309859,  5361.07042254,
        7910.49295775, 10306.29577465, 12224.07042254, 12463.97183099,
       11558.87323944,  9225.8028169 ,  6181.09859155,  3764.05633803,
        2075.78873239,   983.43661972,   479.70422535,   218.84507042,
         118.09859155,    73.53521127,    61.29577465,    51.85915493,
          18.63380282,    22.77464789,    18.76056338,    20.84507042,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ,     0.        ,
           0.        ,     0.        ,     0.        ])

arr_x_coor=np.array([-18.64222222, -18.08666667, -17.53111111, -16.97555556,
       -16.42      , -15.86444444, -15.30888889, -14.75333333,
       -14.19777778, -13.64222222, -13.08666667, -12.53111111,
       -11.97555556, -11.42      , -10.86444444, -10.30888889,
        -9.75333333,  -9.19777778,  -8.64222222,  -8.08666667,
        -7.53111111,  -6.97555556,  -6.42      ,  -5.86444444,
        -5.30888889,  -4.75333333,  -4.19777778,  -3.64222222,
        -3.08666667,  -2.53111111,  -1.97555556,  -1.42      ,
        -0.86444444,  -0.30888889,   0.24666667,   0.80222222,
         1.35777778,   1.91333333,   2.46888889,   3.02444444,
         3.58      ,   4.13555556,   4.69111111,   5.24666667,
         5.80222222,   6.35777778,   6.91333333,   7.46888889,
         8.02444444,   8.58      ,   9.13555556,   9.69111111,
        10.24666667,  10.80222222,  11.35777778,  11.91333333,
        12.46888889,  13.02444444,  13.58      ,  14.13555556,
        14.69111111,  15.24666667,  15.80222222,  16.35777778,
        16.91333333,  17.46888889,  18.02444444,  18.58      ,
        19.13555556,  19.69111111,  20.24666667])

arr_x_coor_uncertainty=[0.]*71
arr_x_uncertainty=[0.]*71

for ii in np.arange(0,71):
    arr_x_coor_uncertainty[ii]=un.ufloat(arr_x_coor[ii],0.555/2)
    arr_x_uncertainty[ii]=un.ufloat(arr_x[ii],np.sqrt(arr_x[ii]))

wrapped_fit=un.wrap(curve_fit_float)

# fit curve with a Gaussian
mean_X=(arr_x_coor[np.size(arr_x_coor)-1]+arr_x_coor[0])/2
sigma_X=1.0
wrapped_fit(arr_x_coor_uncertainty,arr_x_uncertainty,[max(arr_x),mean_X,sigma_X])

运行时抛出如下报错：

File "/opt/anaconda3/lib/python3.8/site-packages/uncertainties/core.py", line 669, in f_with_affine_output
    return f(*args, **kwargs)

  File "/Users/huangze/PhD/academics/post_PSI_data_analysis/debug_odr_uncertainty_package.py", line 18, in curve_fit_float
    popt_X,pcov_X = curve_fit(gaus_fit,x_arr,y_arr,p0)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 734, in curve_fit
    ydata = np.asarray_chkfinite(ydata, float)

  File "/opt/anaconda3/lib/python3.8/site-packages/numpy/lib/function_base.py", line 486, in asarray_chkfinite
    a = asarray(a, dtype=dtype, order=order)

  File "/opt/anaconda3/lib/python3.8/site-packages/numpy/core/_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)

  File "/opt/anaconda3/lib/python3.8/site-packages/uncertainties/core.py", line 2712, in raise_error
    raise TypeError("can't convert an affine function (<class 'uncertainties.core.Variable'>) to float; use x.nominal_value")

TypeError: can't convert an affine function (<class 'uncertainties.core.Variable'>) to float; use x.nominal_value

1. nominal_value的调用位置 #

报错核心原因是scipy.optimize.curve_fit只接受纯浮点型数值输入，无法直接处理uncertainties生成的带不确定度的UFloat对象。需要在传入curve_fit之前，在curve_fit_float函数内部提取输入数组中所有UFloat的标称值，转换为纯浮点数组即可：

def curve_fit_float(x_arr,y_arr,p0):
    # 提取UFloat标称值转为纯浮点数组
    x_nom = np.array([v.nominal_value for v in x_arr])
    y_nom = np.array([v.nominal_value for v in y_arr])
    popt_X,pcov_X = curve_fit(gaus_fit,x_nom,y_nom,p0)
    return popt_X[2]

补充：高斯函数定义里2*sigma**2*sigma/abs(sigma)可以简化为2*sigma**2，sigma/abs(sigma)的结果恒为±1，不影响高斯函数的标准形式，还可以避免拟合迭代中sigma取负值时指数分母为负的异常。

2. uncertainties与curve_fit的适配性 #

二者不能直接配合实现可靠的不确定度传播：

• un.wrap的设计目标是对纯解析浮点运算的函数做自动微分，实现输入不确定度到输出的线性传播。但curve_fit是迭代优化算法，内部包含条件判断、迭代终止等非解析逻辑，不属于纯数值计算函数，即使修复了类型错误，un.wrap也无法正确穿透迭代过程传播不确定度，得到的结果没有参考价值。
• 如果只考虑y方向的测量不确定度，不需要额外使用uncertainties：curve_fit原生支持传入sigma参数指定y值的不确定度，返回的pcov协方差矩阵对角线元素开平方，就是各拟合参数的标准不确定度。
• 如果需要同时考虑x、y两个方向的测量不确定度，不要使用curve_fit，应改用scipy.optimize.ODR（正交距离回归）做拟合，拟合完成后直接从输出的协方差矩阵提取参数不确定度即可。

内容的提问来源于stack exchange，提问作者Ze Huang