scipy.optimize.curve_fit 不合适

Question

数据集的 scipy.optimize.curve_fit 拟合不够好，从前两个峰值和第三个峰值之间可以看出。

我的代码如下

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
    
x = [58.1]
for j in x:
    if j < 71.9:
        x.append(0.05 + j)
    
y = np.array([12,10,16,18,13,14,15,21,26,16,16,16,15,15,31,22,23,22,23,31,31,36,29,30,44,43,45,39,67,67,62,83,79,92,93,128,177,187,251,289,334,414,509,549,623,680,751,769,816,742,701,721,599,540,498,425,339,293,212,172,157,142,115,95,84,86,81,78,58,47,52,62,47,44,39,44,50,34,28,36,36,28,34,24,28,28,37,40,33,45,44,32,32,47,41,46,41,56,46,50,59,71,72,63,62,77,89,62,93,96,114,97,123,140,137,146,138,187,213,219,239,255,287,314,336,386,401,434,467,499,534,635,597,588,597,671,725,693,751,759,727,758,763,720,706,655,650,595,608,565,522,511,433,499,448,387,371,365,324,312,326,306,292,266,282,238,265,278,242,209,226,254,236,233,232,225,220,205,209,234,214,234,234,205,194,223,201,226,214,234,241,226,269,272,303,261,272,277,284,308,324,286,307,329,281,280,257,227,239,234,221,181,195,190,164,158,131,118,118,106,108,70,104,74,76,70,52,57,56,36,44,40,48,40,33,23,28,25,24,20,20,24,31,19,23,24,15,18,12,14,11,14,14,21,20,14,12,10,11,13,5,14,11,10,11,8,8,6,10,5,11,13,8,14,6,8,3,7])

def _1_PseudoVoigt(x, alpha,amp,cen,sigma,b):
    return (1-alpha)*amp*(1/(sigma*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x-cen)/sigma)**2))) + \
            alpha*amp/np.pi*(sigma/((x-cen)**2+sigma**2)) + b

def _PseudoVoigt(x, alpha1,amp1,cen1,sigma1,alpha2,amp2,cen2,sigma2,alpha3,amp3,cen3,sigma3,b):
    return _1_PseudoVoigt(x,alpha1,amp1,cen1,sigma1,b) + _1_PseudoVoigt(x,alpha2,amp2,cen2,sigma2,b) + \
            _1_PseudoVoigt(x,alpha3,amp3,cen3,sigma3,b)

bounds_min = [0,0,60,0, 0,0,64,0, 0.5,0,68,0, 1]
bounds_max = [1,np.inf,62,1, 1,np.inf,66,1, 1,np.inf,69,1, np.inf]

alpha1,amp1,cen1,sigma1 = 0.5,600,60.4,1
alpha2,amp2,cen2,sigma2 = 0.5,1000,64.9,1
alpha3,amp3,cen3,sigma3 = 0.5,350,68.3,0.7
b = 1
    
popt_PseudoVoigt, pcov_PseudoVoigt = curve_fit(lambda x,alpha1,amp1,cen1,sigma1,alpha2,amp2,cen2,sigma2,alpha3,amp3,cen3,sigma3,b: \
                    _PseudoVoigt(x,alpha1,amp1,cen1,sigma1,alpha2,amp2,cen2,sigma2,alpha3,amp3,cen3,sigma3,b), \
                    x, y, \
                    p0=[alpha1,amp1,cen1,sigma1,alpha2,amp2,cen2,sigma2,alpha3,amp3,cen3,sigma3,b], bounds = (bounds_min,bounds_max))

pars_1 = np.append(popt_PseudoVoigt[0:4], popt_PseudoVoigt[12])
pars_2 = np.append(popt_PseudoVoigt[4:8], popt_PseudoVoigt[12])
pars_3 = np.append(popt_PseudoVoigt[8:12], popt_PseudoVoigt[12])

PseudoVoigt_peak_1 = _1_PseudoVoigt(x, *pars_1)
PseudoVoigt_peak_2 = _1_PseudoVoigt(x, *pars_2)
PseudoVoigt_peak_3 = _1_PseudoVoigt(x, *pars_3)

plt.figure(figsize=[10,6])
plt.plot(x,y, marker='.', markersize='5', linestyle='none')
plt.plot(x, PseudoVoigt_peak_1)
plt.plot(x, PseudoVoigt_peak_2)
plt.plot(x, PseudoVoigt_peak_3)
plt.plot(x,_PseudoVoigt(x, *popt_PseudoVoigt), 'k--')
plt.show()

这导致了这张图

可以看出前两个峰值和第三个峰值之间的拟合不太好。有办法改善吗？或者我应该尝试另一个像 lmfit 这样的软件包？

非常感谢您的帮助！

Answer 1

TL；博士

问题更多在于数据而不是拟合程序。主要是第三个峰左侧没有足够的数据（共流出太多）。

此外，通过放松对第三个峰值（中心和 alpha 参数）的一些限制，可以进一步改善拟合。

最后，第三个峰可能不遵循 Voigt 模型，因此即使正确洗脱也根本不适合。

MCVE

让我们展示您能够将一些合成数据集与您的程序相匹配：

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from sklearn.metrics import r2_score

x = np.arange(58.1, 71.96, 0.05)
y = np.array([12, 10, 16, 18, 13, 14, 15, 21, 26, 16, 16, 16, 15, 15, 31, 22, 23, 22, 23, 31, 31, 36, 29, 30, 44, 43, 45, 39, 67, 67, 62, 83, 79, 92, 93, 128, 177, 187, 251, 289, 334, 414, 509, 549, 623, 680, 751, 769, 816, 742, 701, 721, 599, 540, 498, 425, 339, 293, 212, 172, 157, 142, 115, 95, 84, 86, 81, 78, 58, 47, 52, 62, 47, 44, 39, 44, 50, 34, 28, 36, 36, 28, 34, 24, 28, 28, 37, 40, 33, 45, 44, 32, 32, 47, 41, 46, 41, 56, 46, 50, 59, 71, 72, 63, 62, 77, 89, 62, 93, 96, 114, 97, 123, 140, 137, 146, 138, 187, 213, 219, 239, 255, 287, 314, 336, 386, 401, 434, 467, 499, 534, 635, 597, 588, 597, 671, 725, 693, 751, 759, 727, 758, 763, 720, 706, 655, 650, 595, 608, 565, 522, 511, 433, 499, 448, 387, 371, 365, 324, 312, 326, 306, 292, 266, 282, 238, 265, 278, 242, 209, 226, 254, 236, 233, 232, 225, 220, 205, 209, 234, 214, 234, 234, 205, 194, 223, 201, 226, 214, 234, 241, 226, 269, 272, 303, 261, 272, 277, 284, 308, 324, 286, 307, 329, 281, 280, 257, 227, 239, 234, 221, 181, 195, 190, 164, 158, 131, 118, 118, 106, 108, 70, 104, 74, 76, 70, 52, 57, 56, 36, 44, 40, 48, 40, 33, 23, 28, 25, 24, 20, 20, 24, 31, 19, 23, 24, 15, 18, 12, 14, 11, 14, 14, 21, 20, 14, 12, 10, 11, 13, 5, 14, 11, 10, 11, 8, 8, 6, 10, 5, 11, 13, 8, 14, 6, 8, 3, 7])

def peak(x, alpha, A, x0, s, b):
    return A * (
        (1 - alpha)*(1 / (s * (np.sqrt(2 * np.pi)))) * np.exp(-0.5 * np.power((x - x0) / s, 2)) + \
        alpha / np.pi * (s / (np.power((x - x0), 2) + np.power(s, 2)))
    ) + b

def model(x, alpha1, A1, x01, s1, alpha2, A2, x02, s2, alpha3, A3, x03, s3, b):
    return (
          peak(x, alpha1, A1, x01, s1, b)
        + peak(x, alpha2, A2, x02, s2, b)
        + peak(x, alpha3, A3, x03, s3, b)
    )

bounds_min = [
    0, 0, 60, 0,
    0, 0, 64, 0,
    0, 0, 67, 0,  # Enlarge alpha and x0
    -10,
]
bounds_max = [
    1, np.inf, 62, 2,
    1, np.inf, 66, 2,
    1, np.inf, 70, 2,
    10,
]

我们使用一些可以模仿您的数据的参数：

p0 = [
    5.47454610e-01, 7.11445399e+02, 6.04956408e+01, 3.22783279e-01,
    6.07553375e-01, 1.48345124e+03, 6.51134457e+01, 7.05897429e-01,
    2.90916120e-02, 5.87380342e+02, 6.78878208e+01, 9.19652025e-01,
    1.40132572e-01,
]

我们创建一个带有光噪声的合成数据集：

np.random.seed(12345)
ys = model(x, *p0)
ys += 0.1 * np.random.normal(size=y.size)

该程序可以正确回归所有参数：

popt, pcov = optimize.curve_fit(
    model, x, ys,
    bounds=(bounds_min, bounds_max)
)
# array([5.47458889e-01, 7.11360396e+02, 6.04956820e+01, 3.22733587e-01,
#        6.07783229e-01, 1.48357521e+03, 6.51134710e+01, 7.05874792e-01,
#        2.77973904e-02, 5.87092002e+02, 6.78878046e+01, 9.19619386e-01,
#        1.42916997e-01])

spopt = np.sqrt(np.diag(pcov))
rsd = spopt / popt
rsd

# array([5.61128239e-04, 1.32138508e-04, 2.96521977e-07, 7.28917004e-05,
#        4.33596797e-04, 1.21671442e-04, 5.40000252e-07, 6.52427558e-05,
#        3.84556129e-02, 4.05365365e-04, 1.70150662e-06, 1.43308769e-04,
#        4.98309017e-02])

但是您可以看到

alpha3

和

b

具有更大的 RSD。

如果我们添加更多噪声，其级别与您的数据集相当，则这两个参数都变得毫无意义：

np.random.seed(12345)
ys = model(x, *p0)
ys += 10 * np.random.normal(size=y.size)

popt
array([ 5.58589007e-01,  7.07259660e+02,  6.04996493e+01,  3.18390870e-01,
        6.36346065e-01,  1.49709562e+03,  6.51149588e+01,  7.03686126e-01,
        1.65684280e-11,  5.79736506e+02,  6.78855558e+01,  9.21699598e-01,
       -5.14052072e-02])

与第三个峰值相比，它强烈暗示要关闭第二个峰值，并在其左侧丢失太多数据而无法正确回归，并且当噪声增加时，这种效果会变得更糟。在这种情况下，第三个峰值确实遵循 Voigt 模型。

如果我们对您的数据集执行相同的过程，我们会得到：

popt, pcov = optimize.curve_fit(
    model, x, y,
    bounds=(bounds_min, bounds_max)
)
# array([5.47454610e-01, 7.11445399e+02, 6.04956408e+01, 3.22783279e-01,
#        6.07553375e-01, 1.48345124e+03, 6.51134457e+01, 7.05897429e-01,
#        2.90916120e-02, 5.87380342e+02, 6.78878208e+01, 9.19652025e-01,
#        1.40132572e-01])

rsd
# array([1.10523592e-01, 2.60282421e-02, 5.84082859e-05, 1.43564134e-02,
#        8.54478232e-02, 2.39650775e-02, 1.06361209e-04, 1.28483646e-02,
#        7.23149551e+00, 7.98235431e-02, 3.35089472e-04, 2.82361968e-02,
#        1.00113097e+01])

我们可以看到上面提到的对

alpha3

的影响。我们还看到

alpha3

和

x03

有点超出您提供的初始边界，因此在您的拟合中，它们达到了约束边界并停止优化。

score = r2_score(y, model(x, *popt))
# 0.9918069181541809

即使缺乏第三个峰值的适应性，分数也还算可以接受。

结论

拟合 13 个参数并非易事。我们证明，比

alpha3

是更敏感的参数，主要是因为共流出太强。

你有几个选择：

• 减少数据集上的噪声，以将其对

  alpha3

  的影响降低到可接受的水平；
• 提高峰分辨率以减少共洗脱（在第三个峰左侧获得更多点）；
• 更改第三峰的模型。