Dane ładowane w języku Python i przeprowadzanie wielu dopasowań Gaussa

Question

Szukałem sposobu na wielokrotne dopasowanie Gaussa do moich danych. Większość przykładów, które do tej pory znalazłem, używa rozkładu normalnego do losowania liczb. Ale interesuje mnie spisek moich danych i sprawdzenie, czy są 1-3 piki.

Mogę to zrobić na jeden szczyt, ale nie wiem, jak to zrobić, aby uzyskać więcej.

Na przykład, mam te dane: http://www.filedropper.com/data_11

Próbowałem, używając lmfit i oczywiście scipy, ale bez ładne wyniki.

Dzięki za pomoc!

Answer 1

Po prostu sparametryzuj funkcje modelu sumy pojedynczych Gaussów. Wybierz dobrą wartość dla początkowego odgadnięcia (jest to naprawdę krytyczny krok), a następnie zmień nieco te liczby na scipy.optimize.

Oto jak można to zrobić:

import numpy as np 
import matplotlib.pyplot as plt 
from scipy import optimize 

data = np.genfromtxt('data.txt') 
def gaussian(x, height, center, width, offset): 
    return height*np.exp(-(x - center)**2/(2*width**2)) + offset 
def three_gaussians(x, h1, c1, w1, h2, c2, w2, h3, c3, w3, offset): 
    return (gaussian(x, h1, c1, w1, offset=0) + 
     gaussian(x, h2, c2, w2, offset=0) + 
     gaussian(x, h3, c3, w3, offset=0) + offset) 

def two_gaussians(x, h1, c1, w1, h2, c2, w2, offset): 
    return three_gaussians(x, h1, c1, w1, h2, c2, w2, 0,0,1, offset) 

errfunc3 = lambda p, x, y: (three_gaussians(x, *p) - y)**2 
errfunc2 = lambda p, x, y: (two_gaussians(x, *p) - y)**2 

guess3 = [0.49, 0.55, 0.01, 0.6, 0.61, 0.01, 1, 0.64, 0.01, 0] # I guess there are 3 peaks, 2 are clear, but between them there seems to be another one, based on the change in slope smoothness there 
guess2 = [0.49, 0.55, 0.01, 1, 0.64, 0.01, 0] # I removed the peak I'm not too sure about 
optim3, success = optimize.leastsq(errfunc3, guess3[:], args=(data[:,0], data[:,1])) 
optim2, success = optimize.leastsq(errfunc2, guess2[:], args=(data[:,0], data[:,1])) 
optim3 

plt.plot(data[:,0], data[:,1], lw=5, c='g', label='measurement') 
plt.plot(data[:,0], three_gaussians(data[:,0], *optim3), 
    lw=3, c='b', label='fit of 3 Gaussians') 
plt.plot(data[:,0], two_gaussians(data[:,0], *optim2), 
    lw=1, c='r', ls='--', label='fit of 2 Gaussians') 
plt.legend(loc='best') 
plt.savefig('result.png') 

Jak widać, nie ma prawie żadnej różnicy pomiędzy tymi dwoma atakami (wizualnie). Więc nie można wiedzieć na pewno, gdyby były 3 Gaussians obecne w źródle czy tylko 2. Jednakże, jeśli trzeba było zgadnąć, a następnie sprawdzić na najmniejszą resztę:

err3 = np.sqrt(errfunc3(optim3, data[:,0], data[:,1])).sum() 
err2 = np.sqrt(errfunc2(optim2, data[:,0], data[:,1])).sum() 
print('Residual error when fitting 3 Gaussians: {}\n' 
    'Residual error when fitting 2 Gaussians: {}'.format(err3, err2)) 
# Residual error when fitting 3 Gaussians: 3.52000910965 
# Residual error when fitting 2 Gaussians: 3.82054499044 

W tym przypadku 3 Gaussians daje lepszy wynik, ale moje wstępne domysły były dość trafne.