Гауссовская подгонка к шумным данным с использованием Curve_fit

Question

У меня проблемы с подгонкой гауссовости к моим данным. В настоящее время вывод для моего кода выглядит как this . Где оранжевый - это данные, синий - это гауссовская подгонка, а зеленый - встроенный гауссовский установщик, однако я не буду его использовать, поскольку он никогда не начинается с нуля, и у меня нет доступа к коду. Я бы хотел, чтобы мой вывод выглядел примерно так: this , где нарисованное красным цветом соответствует подгонка по Гауссу.

Я попытался прочитать документацию по Curve_Fit, однако в лучшем случае получаю подгонку, которая выглядит как это , которое подходит для всех данных, однако, это нежелательно, поскольку меня интересует только центральный пик, который является моей главной проблемой - я не знаю, как настроить кривую в соответствии с гауссианой на центральном пике как на втором изображении.

Я рассмотрел использование весовой функции, например np.random.choice (), или просмотр максимального значения файла данных, а затем просмотр второй производной по обе стороны от центрального пика, чтобы увидеть, где происходят изменения перегиб, но я не уверен, как лучше это реализовать.

Как бы я лучше go об этом? Я много гуглил, но не могу понять, как изменить кривую в соответствии с моими потребностями.

Спасибо за любые указатели!

Это файл данных.

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()

Answer 1

Я думаю, что проблема может заключаться в том, что данные не полностью гауссовоподобны. Кажется, у вас есть какая-то функция Airy / sin c из-за временного разрешения вашего инструмента сбора данных. Тем не менее, если вас интересует только центр, вы все равно можете установить его, используя один гауссиан:

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

Если вы хотите чтобы быть немного более точным, я могу представить квадратную функцию sin c для пика и широкий гауссовский сдвиг. Подгонка выглядит лучше, но это действительно зависит от того, что на самом деле представляют данные ...

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