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

person Harpal    schedule 13.04.2012    source источник
comment
В данном случае это было бы довольно сложно, так как два пика довольно близки друг к другу - для меньшего «гауссова» нет определенного пика. Обычно (я думаю) нужно идентифицировать все интересующие пики, а затем перебирать каждый пик, маскируя все остальные пики и подгоняя каждый пик. Полная подгонка тогда является суммой всех этих подгонок. Похоже, вам нужно определить большой пик и его протяженность, а затем замаскировать его из данных перед подгонкой к меньшему пику.   -  person Chris    schedule 13.04.2012


Ответы (3)


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

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

Параметры (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()

Результаты кода.

person Usagi    schedule 13.04.2012
comment
Итак, я предполагаю, что для n гауссов мне нужно будет сложить n гауссовских функций вместе и вычесть их из данных? - person Harpal; 14.04.2012
comment
@Харпал - Да. Вы можете изменить код, чтобы использовать n кривых. Я бы просто убедился, что закодировал алгоритм таким образом, чтобы никакие две кривые не имели одинаковое среднее значение. - person Usagi; 17.04.2012
comment
Строка y_est = норма(х, plsq[0][0], plsq[0][2]) + норма(х, plsq[0][1], plsq[0][3]) должна быть y_est = норма (x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3]); в вашем примере это не очевидно, потому что одно из средств равно нулю. Отредактировал это. В остальном отличное решение :) - person Kyle; 21.06.2013

Вы можете использовать смешанные модели Гаусса из 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
person bougui    schedule 04.10.2013
comment
Это будет соответствовать гистограмме данных, а не самим данным. - person Rob; 11.01.2016

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

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

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

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

person andrew cooke    schedule 14.04.2012
comment
Несмотря на колючесть, это действительно кажется мне правдоподобным. Если вы можете подогнать всю модель за один раз, это дает бесчисленные преимущества. Проголосовал. - person nes1983; 14.04.2012