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

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

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

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

Ответы

Ответ 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()

Ответ 2

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

Ответ 3

coeffs 0 и 4 вырождены - в данных, которые могут решить между ними, нет абсолютно ничего. вы должны использовать один параметр нулевого уровня вместо двух (т.е. удалить один из них из вашего кода). это, вероятно, то, что останавливает вашу пригодность (игнорируйте комментарии здесь, говоря, что это невозможно), в этих данных явно есть как минимум два пика, и вы, безусловно, должны быть в состоянии соответствовать этому).

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

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

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