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/6097EBE81EE0C29EE95A497128C1C2E420220704110110130/86BF2D

这是我下面的代码。有人可以告诉我问题是什么，我该如何解决？ /I.sstatic.net/nzkim.png“ alt =”在此处输入图像说明”> 在图的图片上，蓝色情节是我的实验性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