python - scipy curve_fi 返回初始参数估计

Question

我正在尝试使用 scipy curve_fit 函数将高斯函数拟合到我的数据以估计理论功率谱密度。这样做时，curve_fit 函数总是返回初始参数 (p0=[1,1,1]) ，从而告诉我拟合不起作用。我不知道问题出在哪里。我正在使用 Windows 11 上的 anaconda 发行版中的 python 3.9 (spyder 5.1.5)。这是指向数据文件的 Wetransfer 链接 https://wetransfer.com/downloads/6097ebe81ee0c29ee95a497128c1c2e420220704110130/86bf2d

下面是我的代码。有人可以告诉我问题是什么，我该如何解决？ 在图的图片上，蓝色图是我的实验PSD，橙色图是拟合的结果。

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.constants as cst

File = np.loadtxt('test5.dat')
X = File[:, 1]
Y = File[:, 2]


f_sample = 50000
time=[]
for i in range(1,len(X)+1):
    t=i*(1/f_sample)
    time= np.append(time,t)
N = X.shape[0]  # number of observation
N1=int(N/2)
delta_t = time[2] - time[1]
T_mes = N * delta_t
freq = np.arange(1/T_mes, (N+1)/T_mes, 1/T_mes)
freq=freq[0:N1]
fNyq = f_sample/2  # Nyquist frequency
nb = 350
freq_block = []

#  discrete fourier tansform
X_ft = delta_t*np.fft.fft(X, n=N)
X_ft=X_ft[0:N1]

plt.figure()
plt.plot(time, X)
plt.xlabel('t [s]')
plt.ylabel('x [micro m]')

# Experimental power spectrum on both raw and blocked data
PSD_X_exp = (np.abs(X_ft)**2/T_mes)
PSD_X_exp_b = []

STD_PSD_X_exp_b = []
for i in range(0, N1+2, nb):
    freq_b = np.array(freq[i:i+nb])  # i-nb:i
    psd_b = np.array(PSD_X_exp[i:i+nb])
    freq_block = np.append(freq_block, (1/nb)*np.sum(freq_b))
 
    PSD_X_exp_b = np.append(PSD_X_exp_b, (1/nb)*np.sum(psd_b))
STD_PSD_X_exp_b = np.append(STD_PSD_X_exp_b, PSD_X_exp_b/np.sqrt(nb))

plt.figure()
plt.loglog(freq, PSD_X_exp)
plt.legend(['Raw Experimental PSD'])
plt.xlabel('f [Hz]')
plt.ylabel('PSD')
plt.figure()
plt.loglog(freq_block, PSD_X_exp_b)
plt.legend(['Experimental PSD after blocking'])
plt.xlabel('f [Hz]')
plt.ylabel('PSD')

kB = cst.k  # Boltzmann constant [m^2kg/s^2K]
T = 273.15 + 25  # Temperature [K]
r = (2.8 / 2) * 1e-6  # Particle radius [m]
v = 0.00002414 * 10 ** (247.8 / (-140 + T))  # Water viscosity [Pa*s]
gamma = np.pi * 6 * r * v  # [m*Pa*s]
Do = kB*T/gamma  # expected value for D
f3db_o = 50000  # expected value for f3db
fc_o = 300  # expected value pour fc

n = np.arange(-10,11)


def theo_spectrum_lorentzian_filter(x, D_, fc_, f3db_):
    PSD_theo=[]
    for i in range(0,len(x)):
        # print(i)
        psd_theo=np.sum((((D_*Do)/2*math.pi**2)/((fc_*fc_o)**2+(x[i]+n*f_sample)
                    ** 2))*(1/(1+((x[i]+n*f_sample)/(f3db_*f3db_o))**2)))
       
        PSD_theo= np.append(PSD_theo,psd_theo)
            

    return PSD_theo

popt, pcov = curve_fit(theo_spectrum_lorentzian_filter, freq_block, PSD_X_exp_b, p0=[1, 1, 1], sigma=STD_PSD_X_exp_b, absolute_sigma=True, check_finite=True,bounds=(0.1, 10),  method='trf', jac=None)

D_, fc_, f3db_ = popt
D1 = D_*Do
fc1 = fc_*fc_o
f3db1 = f3db_*f3db_o

print('Diffusion constant D = ', D1, ' Corner frequency fc= ',fc1, 'f3db(diode,eff)= ', f3db1)

Answer 1

我相信我已经成功地拟合了您的数据。这是我采用的方法。

首先，我绘制了您的模型(使用 popt=[1, 1, 1])和您拥有的数据。我注意到您的数据明显低于模型。然后我开始摆弄参数。我想向上推模型。我通过将 popt[0] 乘以越来越大的值来做到这一点。我最终以 1E13 作为大概值。请注意，我不知道您的模型在物理上是否可行。然后我临时操纵了你的拟合函数，将 D_ 乘以 1E13 并运行了你的代码。我很合身:

所以我认为这是 1) 不合适的起始值和 2) 不合适的界限的问题。在你的位置上，我会修改这个模型，检查单位是否有任何问题等等。

这是我用来尝试拟合您的模型的内容:

plt.figure()
plt.loglog(freq_block[:170], PSD_X_exp_b[:170], label='Exp')
plt.loglog(freq_block[:170], 
           theo_spectrum_lorentzian_filter(
               freq_block[:170],
               1E13*popt[0], popt[1], popt[2]),
           label='model'
          )
plt.xlabel('f [Hz]')
plt.ylabel('PSD')
plt.legend()

我将数据限制在 170 点，因为有一些奇怪的反向值让我感到不舒服。如果我是你，我会重新检查它们。

这是我使用的模型代码。我没有更改 curve_fit 调用(除了将 x 限制为 :170。

def theo_spectrum_lorentzian_filter(x, D_, fc_, f3db_):
    PSD_theo=[]
    D_ = 1E13*D_  # I only changed here
    for i in range(0,len(x)):
        psd_theo=np.sum((((D_*Do)/2*math.pi**2)/((fc_*fc_o)**2+(x[i]+n*f_sample)
                    ** 2))*(1/(1+((x[i]+n*f_sample)/(f3db_*f3db_o))**2)))
       
        PSD_theo= np.append(PSD_theo,psd_theo)
            
    return PSD_theo