python - 使用curve_fit对噪声数据进行高斯拟合
问题描述
我在将高斯拟合到我的数据时遇到问题。目前我的代码的输出如下 所示。橙色是数据,蓝色是高斯拟合,绿色是内置的高斯拟合器,但是我不希望使用它,因为它从不完全从零开始,而且我无权访问代码。我希望我的输出看起来像这样,其中红色绘制的是高斯拟合。
我已经尝试阅读有关curve_fit文档的信息,但是充其量我得到的拟合看起来像这样,它适合所有数据,但是,这是不可取的,因为我只对中心峰感兴趣,这是我的主要问题-我不知道如何让curve_fit像第二张图像一样在中心峰上拟合高斯。
我考虑过使用像 np.random.choice() 这样的权重函数,或者查看数据文件的最大值,然后查看中心峰两侧的二阶导数,以查看拐点变化的位置,但不确定如何最好来实现这一点。
我最好怎么做?我已经做了很多谷歌搜索,但不能完全改变curve_fit以满足我的需要。
为任何指针干杯!
这是一个数据文件。
https://drive.google.com/file/d/1qrAkD74U6L46GoGnvMiUHdPuLEToS6Pv/view?usp=sharing
这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from matplotlib.pyplot import figure
plt.close('all')
fpathB4 = 'E:\.1. Work - Current Projects + Old Projects\Current Projects\PF 4MHz Laser System\.8. 1050 SYSTEM\AC traces'
fpath = fpathB4.replace('\\','/') + ('/')
filename = '300'
with open(fpath+filename) as f:
dataraw = f.readlines()
FWHM = dataraw[8].split(':')[1].split()[0]
FWHM = np.float(FWHM)
print("##### For AC file -", filename, "#####")
print("Auto-co guess -", FWHM, "ps")
pulseduration = FWHM/np.sqrt(2)
pulseduration = str(pulseduration)
dataraw = dataraw[15:]
print("Pulse duration -", pulseduration, "ps" + "\n")
time = np.array([])
acf1 = np.array([]) #### DATA
fit = np.array([]) #### Gaussian fit
for k in dataraw:
data = k.split()
time = np.append(time, np.float(data[0]))
acf1= np.append(acf1, np.float(data[1]))
fit = np.append(fit, np.float(data[2]))
n = len(time)
y = acf1.copy()
x = time.copy()
mean = sum(x*y)/n
sigma = sum(y*(x-mean)**2)/n
def gaus(x,a,x0,sigma):
return a*np.exp(-(x-x0)**2/(2*sigma**2))
popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])
plt.plot(x,gaus(x,*popt)/np.max(gaus(x,*popt)))
figure(num=1, figsize=(8, 3), dpi=96, facecolor='w', edgecolor='k') # figsize = (length, height)
plt.plot(time, acf1/np.max(acf1), label = 'Data - ' + filename, linewidth = 1)
plt.plot(time, fit/np.max(fit), label = '$FWHM_{{\Delta t}}$ (ps) = ' + pulseduration)
plt.autoscale(enable = True, axis = 'x', tight = True)
plt.title("Auto-Correlation Data")
plt.xlabel("Time (ps)")
plt.ylabel("Intensity (a.u.)")
plt.legend()
解决方案
我认为问题可能在于数据并不完全类似于高斯。由于采集仪器的时间分辨率,您似乎具有某种 Airy/sinc 函数。不过,如果您只对中心感兴趣,您仍然可以使用单个高斯拟合它:
import fitwrap as fw
import pandas as pd
df = pd.read_csv('300', skip_blank_lines=True, skiprows=13, sep='\s+')
def gaussian_no_offset(x, x0=2, sigma=1, amp=300):
return amp*np.exp(-(x-x0)**2/sigma**2)
fw.fit(gaussian_no_offset, df.time, df.acf1)
x0: 2.59158 +/- 0.00828 (0.3%) initial:2
sigma: 0.373 +/- 0.0117 (3.1%) initial:1
amp: 355.02 +/- 9.65 (2.7%) initial:300
如果您想稍微精确一点,我可以想到峰值的 sinc 平方函数和宽高斯偏移。拟合似乎更好,但它实际上取决于数据实际代表什么......
def sinc(x, x0=2.5, amp=300, width=1, amp_g=20, sigma=3):
return amp*(np.sinc((x-x0)/width))**2 + amp_g*np.exp(-(x-x0)**2/sigma**2)
fw.fit(sinc, df.time, df.acf1)
x0: 2.58884 +/- 0.0021 (0.1%) initial:2.5
amp: 303.84 +/- 3.7 (1.2%) initial:300
width: 0.49211 +/- 0.00565 (1.1%) initial:1
amp_g: 81.32 +/- 2.11 (2.6%) initial:20
sigma: 1.512 +/- 0.0351 (2.3%) initial:3
推荐阅读
- flutter - 如何解决“RangeError (index): Invalid value: Only valid value is 0: 1”错误
- c# - 由于单元测试而面试失败 - 我需要了解原因
- sql - 用于避免重复列的 SQL Server 查询
- reactjs - 在其他组件中或在单独的文件中创建反应功能组件
- c++ - 是否有一个 C++ 函数可以将向量分成三个单独的向量?
- google-cloud-platform - Flex 模板发行说明
- javascript - 处理字符串原型中的空类型
- laravel - Laravel 按日期分组
- c# - Range.FitTextWidth 始终返回零
- google-analytics - 是否可以使用 gtag.js 通过自定义参数对 Google Analytics 中的用户进行分组?