我可以强制 lmfit 高斯模型都具有正振幅吗？

如何解决我可以强制 lmfit 高斯模型都具有正振幅吗？

我正在对一些数据拟合高斯曲线，虽然“拟合”非常令人印象深刻，但有时会使用负振幅来实现它。下面的代码生成数据、拟合和组成高斯曲线的图。您会看到一条幅度为 577 的压倒性高斯曲线和另一条幅度为 -570 的高斯曲线，这有助于拟合曲线。

我更喜欢会牺牲一点“精度”来强制所有振幅为正的拟合。有什么办法可以强制模型的参数保持振幅为正？

import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import signal
import math
from lmfit import models

perc = [0,0.000019,0.000294,0.001946,0.007003,0.017276,0.033345,0.055974,0.086622,0.126049,0.173105,0.226678,0.286038,0.351443,0.42241,0.498498,0.578561,0.661742,0.746826,0.83285,0.919105,1.00482,1.08944,1.17218,1.25306,1.33216,1.40968,1.48517,1.55839,1.62919,1.69746,1.7625,1.82327,1.87912,1.93017,1.97694,2.01944,2.05721,2.08996,2.11797,2.14111,2.15795,2.16612,2.16457,2.15517,2.14131,2.125,2.10617,2.08647,2.07406,2.08323,2.12419,2.19053,2.25702,2.29272,2.28612,2.25781,2.24839,2.29069,2.38262,2.48227,2.52351,2.44913,2.24458,1.94567,1.61752,1.32053,1.08224,0.895022,0.730309,0.559126,0.375057,0.205875,0.099411,0.06067,0.063161,0.079791,0.089018,0.073056,0.041268,0.014748,0.002668,0.000197,0]
gsize = [0.04,0.04391,0.048203,0.052916,0.058089,0.063768,0.070002,0.076845,0.084358,0.092605,0.101658,0.111597,0.122507,0.134483,0.147631,0.162064,0.177907,0.1953,0.214393,0.235353,0.258361,0.283619,0.311346,0.341784,0.375198,0.411878,0.452145,0.496347,0.544872,0.59814,0.656615,0.720807,0.791275,0.868632,0.953552,1.04677,1.14911,1.26145,1.38477,1.52015,1.66876,1.8319,2.011,2.2076,2.42342,2.66033,2.92042,3.20592,3.51934,3.8634,4.2411,4.65572,5.11087,5.61052,6.15902,6.76114,7.42212,8.14773,8.94427,9.81869,10.7786,11.8323,12.9891,14.2589,15.6529,17.1832,18.863,20.7071,22.7315,24.9538,27.3934,30.0714,33.0113,36.2385,39.7813,43.6704,47.9397,52.6264,57.7713,63.4192,69.6192,76.4253,83.8969,92.0988,101.103,110.987,121.837,133.748,146.824,161.177,176.935,194.232,213.221,234.066,256.948,282.068,309.644,339.916,373.147,409.626,449.672,493.633,541.892,594.869,653.025,716.866,786.949,863.883,948.338,1041.05,1142.83,1254.55,1377.2,1511.84,1659.64,1821.89]

xinit = gsize
xlog = [math.log(xval) for xval in xinit]
x = np.array(xlog)
y = np.array(perc)

peaks = signal.find_peaks_cwt(y,(2.5,25))             
xstep = x.ptp() / len(x)

model,params = None,None

for i,peak_index in enumerate(peaks):
    this_model = models.GaussianModel(prefix=f'p{1+i:d}_')
    this_params = this_model.make_params(amplitude=y[peak_index],center=x[peak_index],sigma=2*xstep)
    if model is None:
        model = this_model
        params = this_params
    else:
        model += this_model
        params.update(this_params)

result = model.fit(y,params,x=x)
print(result.fit_report())

params_list = ['p1_','p2_','p3_','p4_','p5_','p6_']
plt.plot(figsize=(60,50))
plt.ylim(-4.0,4.0)
for param in params_list:
   try:
                
       center = result.params[param + 'center'].value
       sigma = result.params[param + 'sigma'].value
       mean = center
       standard_deviation = sigma
       amplitude = result.params[param + 'amplitude'].value
       x_values = np.arange(-2,10,0.1)
       y_values = scipy.stats.norm(mean,standard_deviation)
       plt.plot(x_values,(amplitude*y_values.pdf(x_values)))

   except KeyError:
       continue

plt.plot(x,y,label='data')
plt.plot(x,result.best_fit,label='fit')
plt.legend()
plt.show()

解决方法

这里的问题是你有太多的自由度，它会过拟合。

如果你拟合一个只有三个高斯分布的模型

for i,peak_index in enumerate(peaks[:3]):
    this_model = models.GaussianModel(prefix=f'p{1+i:d}_')
    this_params = this_model.make_params(amplitude=y[peak_index],center=x[peak_index],sigma=2*xstep)
    if model is None:
        model = this_model
        params = this_params
    else:
        model += this_model
        params.update(this_params)

它不会使用负振幅

,

是的，您可以将每个“振幅”参数设置为正值。类似的东西：

for i,peak_index in enumerate(peaks):
    this_model = models.GaussianModel(prefix=f'p{1+i:d}_')
    this_params = this_model.make_params(amplitude=y[peak_index],sigma=2*xstep)
    this_params[f'p{1+i:d}_amplitude'].min = 0.  # <--- set min value here!

    if model is None:
        model = this_model
        params = this_params
    else:
        model += this_model
        params.update(this_params)

应该这样做。