使用curve_fit拟合复合高斯函数时np.exp出现溢出错误

问题描述 #

我用以下代码拟合复合高斯函数：

def curve_f(x,a,b,c,d,e,f):
  return (a*np.exp(((x-b)**2)/(2*-c**2)))+(d*np.exp(((x-e)**2)/(2*f**2)))

将该函数传入scipy.optimize.curve_fit处理数据时，出现以下警告：

<ipython-input-27-1baee2b89bf1>:2: RuntimeWarning: overflow encountered in exp
  return (a*np.exp(((x-b)**2)/(2*-c**2)))+(d*np.exp(((x-e)**2)/(2*f**2)))
/usr/local/lib/python3.10/dist-packages/scipy/optimize/_minpack_py.py:862: RuntimeWarning: overflow encountered in square
cost = np.sum(infodict['fvec'] ** 2)

我尝试把数据转为float128，但问题没解决。Y轴数据样例如下：

[   0    0    0    0    0    0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0   -1 -110 1290   33    0
0    0    0    0   -1    1    1    2    2    0   -1   -2   -1    1
1    1    2    2    2    4   -3    5    5    7    3   10    3   11
10    4   16   13   32   23   21   23   37   42   42   43   56   62
41   56   71   59   61   84   82   76   73   78   73   69   84   85
74   65   75   59   74   73   76   74   60   95   63   91   85   82
73   72   76   82   71   85   71   60   53   71   84   95   83   50
88   63   76   68   81   73   80   60   89   78   76   86   73   99
91   94  104   98  114  117  110  122  104  106  121  135  130  161
146  133  172  145  186  202  222  216  205  200  234  219  231  265
269  266  296  286  304  303  325  319  381  328  362  370  339  360
423  405  368  400  411  416  448  454  453  468  480  427  433  490
471  505  500  442  483  482  462  436  457  483  435  431  432  416
410  438  402  422  387  403  366  371  369  339  341  306  327  331
310  296  277  259  268  228  214  225  234  204  208  205  205  176
168  190  192  142  161  126  107  122  112   86  107   92   90   89
78   82   59   69   67   57   34   42   46   44   35   36   37   29
28   29   25   20   20   18   20   17    8   22   12    9    9    7
13    8   11    5    5   14    5    3    6    4    6    9    2    5
5    2    4    5    8    7    5    1    4    3    2    6   -1    3
3    3    2    3    2    3    1    1    2    0    1    2    3    0
1    0    5    3    2    1    1    2    2    0    4    3    1    2
2    4    2    4    0    4    1    2    1    1    1    1    0    2
-1    2    1    1    1    4    0    4    2    3    0   -1    1    3
1    2    0    2    1    1    1    1    0    1    0    1    1    1
1    0    0    1    2    2    0    3    2    1    1    1    0    0
1    1    1    1    1    1    2    2    1    2    0    1   -1    2
1    1    2    1    0    1    1    0]

问题原因 #

溢出警告本质是拟合过程中，优化算法尝试的参数导致np.exp()的输入值过大，超出浮点数表示范围。另外数据里存在负值（比如-110），高斯函数本身非负，这些异常点会增大残差，进一步导致优化过程中参数失控。

解决步骤 #

1. 预处理数据，移除无效负值
   高斯函数输出不可能为负，数据里的负值属于异常值，建议直接替换为0：

   y_data[y_data < 0] = 0
   

2. 提供合理的初始参数猜测
   curve_fit默认初始参数全为1，极易引发exp计算溢出。根据数据分布手动估计初始值：

   • 峰值幅度a、d：参考数据峰值500左右，设为400、100；
   • 峰值位置b、e：数据峰值集中在中间区域，设为120、150；
   • 标准差c、f：根据数据宽度设为30、40。
     示例代码：

   initial_guess = [400, 120, 30, 100, 150, 40]
   popt, pcov = curve_fit(curve_f, x_data, y_data, p0=initial_guess)
   

3. 约束参数范围
   给参数添加边界，避免优化时出现不合理取值（比如标准差必须为正）：

   # 每个参数的(下限, 上限)，对应a,b,c,d,e,f
   bounds = (0, [np.inf, np.inf, np.inf, np.inf, np.inf, np.inf])
   popt, pcov = curve_fit(curve_f, x_data, y_data, p0=initial_guess, bounds=bounds)
   

4. 优化函数计算稳定性
   重构高斯函数写法，把负号移到分子，避免分母出现负数，提升计算可读性和稳定性：

   def curve_f(x,a,b,c,d,e,f):
       term1 = a * np.exp(-(x - b)**2 / (2 * c**2))
       term2 = d * np.exp(-(x - e)**2 / (2 * f**2))
       return term1 + term2
   

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