Python: двухкристальный гауссовский фитинг с нелинейными наименьшими квадратами

Question

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

Синяя линия мои данные и зеленая линия мой текущий нужным. Существует плечо слева от основного пика в моих данных, которые я в настоящее время пытается вписаться, используя следующий код:

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

Этот код работает для меня, при условии, что вы только фитинг функцию, которая является сочетание двух гауссовских распределений.

Я только что сделал функцию остатков, которая добавляет две функции Гаусса, а затем вычитает их из реальных данных.

Параметры (p), которые я передал функции наименьших квадратов Numpy, включают в себя: среднее значение первой гауссовой функции (m), разность средних от первой и второй гауссовых функций (dm, т.е. горизонтальный сдвиг) , стандартное отклонение первого (sd1) и стандартное отклонение второго (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 и 4 являются вырожденными - нет абсолютно ничего в данных, которые могут решить между ними. вы должны использовать один параметр нулевого уровня вместо двух (т. е. удалить один из них из вашего кода). это, вероятно, то, что останавливает вашу посадку (игнорируйте комментарии здесь, говоря, что это невозможно), в этих данных явно есть как минимум два пика, и вы, безусловно, сможете соответствовать этому).

(может быть неясно, почему я предлагаю это, но происходит то, что коэффициенты 0 и 4 могут отменить друг друга. Оба они могут быть равны нулю, или один может быть 100, а другой -100 - либо так как это подходит, это «сбивает с толку» подгоняющую процедуру, которая тратит свое время на то, чтобы выяснить, какими они должны быть, когда нет единого правильного ответа, потому что независимо от того, какое значение оно имеет, другое может быть просто отрицательный результат, а подгонка будет одинаковой).

Фактически, из сюжета, похоже, не может быть необходимости в нулевом уровне вообще. я бы постарался сбросить оба из них и посмотреть, как выглядит подгонка.

также нет необходимости устанавливать коэффициенты 1 и 5 (или нулевую точку) в наименьших квадратах. вместо этого, поскольку модель линейна в тех, которые вы могли бы рассчитать их значения для каждого цикла. это сделает вещи быстрее, но не критично. я просто заметил, что вы говорите, что ваша математика не так хороша, поэтому, вероятно, проигнорируйте это.

Answer 3

Вы можете использовать гауссовые смеси из 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]) 

Вы также можете использовать функцию ниже, чтобы соответствовать числу гауссовых вы хотите с ncomp параметром:

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