curve_fit 的 np.exp 函数中出现溢出错误

Question

我正在尝试使用以下代码来拟合组合高斯：

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]

Answer 1

工作正常。如 StackOverflow 高斯拟合问题（一遍又一遍，一遍又一遍）所示，添加良好的边界和初始猜测。

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


def curve_f(x: np.ndarray, a: float, b: float, c: float, d: float, e: float, f: float, g: float) -> np.ndarray:
    return (
        1/np.sqrt(2*np.pi) *
        (
            + a/c * np.exp(-0.5*((x - b)/c)**2)
            + d/f * np.exp(-0.5*((x - e)/f)**2)
        )
        + g
    )


y = 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, -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])

# Ignore massive spike at beginning
y = y[30:]

x = np.arange(len(y))
x0 = [5000, 50, 25, 25000, 150, 25, 0]
popt, *_ = curve_fit(
    f=curve_f,
    xdata=x,
    ydata=y,
    p0=x0,
    bounds=Bounds(
        lb=[1e-2,   0, 1e-2, 1e-2,   0, 1e-2, -1000],
        ub=[1e+6, 400, 1000, 1e+6, 400, 1000,  1000],
    )
)
print(popt)

plt.scatter(x, y, label='experiment', marker='+', color='green')
plt.plot(x, curve_f(x, *x0), label='guess')
plt.plot(x, curve_f(x, *popt), label='fit')
plt.legend()
plt.show()