Python: dopasowanie krzywej gaussowskiej w dwóch krzywiznach z nieliniowymi najmniejszymi kwadratami

Question

Moja wiedza na temat matematyki jest ograniczona, dlatego prawdopodobnie utknąłem. Mam widma, do których próbuję dopasować dwa szczyty Gaussa. Mogę zmieścić się na największym piku, ale nie mogę zmieścić się na najmniejszym szczycie. Rozumiem, że muszę podsumować funkcję Gaussa dla dwóch szczytów, ale nie wiem, gdzie popełniłem błąd. Obraz mojego prąd wyjściowy jest pokazany:

Niebieska linia to moje dane i zielona linia jest mój obecny dopasowanie. Jest barku na lewo od głównego piku w moim danych, które mam obecnie stara się dopasować, stosując następujący kod:

import matplotlib.pyplot as pt 
import numpy as np 
from scipy.optimize import leastsq 
from pylab import * 

time = [] 
counts = [] 


for i in open('/some/folder/to/file.txt', 'r'): 
    segs = i.split() 
    time.append(float(segs[0])) 
    counts.append(segs[1]) 

time_array = arange(len(time), dtype=float) 
counts_array = arange(len(counts)) 
time_array[0:] = time 
counts_array[0:] = counts 


def model(time_array0, coeffs0): 
    a = coeffs0[0] + coeffs0[1] * np.exp(- ((time_array0-coeffs0[2])/coeffs0[3])**2) 
    b = coeffs0[4] + coeffs0[5] * np.exp(- ((time_array0-coeffs0[6])/coeffs0[7])**2) 
    c = a+b 
    return c 


def residuals(coeffs, counts_array, time_array): 
    return counts_array - model(time_array, coeffs) 

# 0 = baseline, 1 = amplitude, 2 = centre, 3 = width 
peak1 = np.array([0,6337,16.2,4.47,0,2300,13.5,2], dtype=float) 
#peak2 = np.array([0,2300,13.5,2], dtype=float) 

x, flag = leastsq(residuals, peak1, args=(counts_array, time_array)) 
#z, flag = leastsq(residuals, peak2, args=(counts_array, time_array)) 

plt.plot(time_array, counts_array) 
plt.plot(time_array, model(time_array, x), color = 'g') 
#plt.plot(time_array, model(time_array, z), color = 'r') 
plt.show() 

Answer 1

Ten kod pracował dla mnie zapewnienie, że jesteś sylwetkę tylko funkcję, która jest połączenie dwóch rozkładów Gaussa.

Właśnie stworzyłem funkcję residuals, która dodaje dwie funkcje gaussowskie, a następnie odejmuje je od rzeczywistych danych.

Parametry (p) przekazywane do funkcji najmniejszych kwadratów Numpy'ego obejmują: średnią pierwszej funkcji Gaussa (m), różnicę średniej z pierwszej i drugiej funkcji Gaussa (dm, czyli przesunięcie w poziomie) , standardowe odchylenie pierwszego (sd1) i standardowe odchylenie drugiego (sd2).

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

###################################### 
# Setting up test data 
def norm(x, mean, sd): 
    norm = [] 
    for i in range(x.size): 
    norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))] 
    return np.array(norm) 

mean1, mean2 = 0, -2 
std1, std2 = 0.5, 1 

x = np.linspace(-20, 20, 500) 
y_real = norm(x, mean1, std1) + norm(x, mean2, std2) 

###################################### 
# Solving 
m, dm, sd1, sd2 = [5, 10, 1, 1] 
p = [m, dm, sd1, sd2] # Initial guesses for leastsq 
y_init = norm(x, m, sd1) + norm(x, m + dm, sd2) # For final comparison plot 

def res(p, y, x): 
    m, dm, sd1, sd2 = p 
    m1 = m 
    m2 = m1 + dm 
    y_fit = norm(x, m1, sd1) + norm(x, m2, sd2) 
    err = y - y_fit 
    return err 

plsq = leastsq(res, p, args = (y_real, x)) 

y_est = norm(x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3]) 

plt.plot(x, y_real, label='Real Data') 
plt.plot(x, y_init, 'r.', label='Starting Guess') 
plt.plot(x, y_est, 'g.', label='Fitted') 
plt.legend() 
plt.show() 

Answer 2

coeffs 0 i 4 są zdegenerowane - nie ma absolutnie nic w danych, które mogą zdecydować się między nimi. powinieneś użyć jednego parametru poziomu zerowego zamiast dwóch (tj. usuń jeden z nich z kodu). to prawdopodobnie powstrzymuje twój atak (zignoruj ​​komentarze, mówiąc, że nie jest to możliwe - są wyraźnie co najmniej dwa szczyty w tych danych i na pewno będziesz w stanie się do tego dopasować).

(może nie być jasne, dlaczego sugeruję to, ale co się dzieje, to, że coeffs 0 i 4 mogą anulować siebie nawzajem, mogą być równe zero lub jeden może być 100, a drugi -100 - albo sposób, dopasowanie jest równie dobre, to "myli" rutynę dopasowania, która spędza czas próbując ustalić, czym powinny być, kiedy nie ma jednej prawidłowej odpowiedzi, ponieważ bez względu na to, kim jest, druga może być po prostu ujemny, a dopasowanie będzie takie samo).

w rzeczywistości, z działki wygląda na to, że w ogóle może nie być potrzeby stosowania poziomu zerowego. chciałbym spróbować rzucić oba i zobaczyć, jak wygląda dopasowanie.

również, nie ma potrzeby, aby dopasować coeffs 1 i 5 (lub punkt zerowy) w najmniejszych kwadratach. zamiast tego, ponieważ model jest liniowy, można obliczyć ich wartości w każdej pętli. to sprawi, że rzeczy będą szybsze, ale nie krytyczne. Właśnie zauważyłem, że mówisz, że twoje matematyki nie jest tak dobre, więc prawdopodobnie zignorować to.

Answer 3

Można użyć Gaussa modele Mieszaninę z scikit-learn:

from sklearn import mixture 
import matplotlib.pyplot 
import matplotlib.mlab 
import numpy as np 
clf = mixture.GMM(n_components=2, covariance_type='full') 
clf.fit(yourdata) 
m1, m2 = clf.means_ 
w1, w2 = clf.weights_ 
c1, c2 = clf.covars_ 
histdist = matplotlib.pyplot.hist(yourdata, 100, normed=True) 
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=3) 
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=3) 
plotgauss1(histdist[1]) 
plotgauss2(histdist[1]) 

Można również skorzystać z funkcji poniżej, aby dopasować liczbę Gaussa chcesz z ncomp parametru:

from sklearn import mixture 
%pylab 

def fit_mixture(data, ncomp=2, doplot=False): 
    clf = mixture.GMM(n_components=ncomp, covariance_type='full') 
    clf.fit(data) 
    ml = clf.means_ 
    wl = clf.weights_ 
    cl = clf.covars_ 
    ms = [m[0] for m in ml] 
    cs = [numpy.sqrt(c[0][0]) for c in cl] 
    ws = [w for w in wl] 
    if doplot == True: 
     histo = hist(data, 200, normed=True) 
     for w, m, c in zip(ws, ms, cs): 
      plot(histo[1],w*matplotlib.mlab.normpdf(histo[1],m,np.sqrt(c)), linewidth=3) 
    return ms, cs, ws