Вопрос по scipy, gaussian, python, least-squares – Python: двухкривая гауссова аппроксимация с нелинейными наименьшими квадратами

17

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

Current Output

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

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

Ваш Ответ

3   ответа
16

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

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

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

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

Вы можете использовать модели гауссовой смеси из Scikit учиться:

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
Это будет соответствовать гистограмме данных, а не самим данным. Rob
4

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

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

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

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

Неприятность, несмотря на это, мне это кажется правдоподобным. Если вы можете приспособить всю свою модель за один раз, это имеет бесчисленные преимущества. Upvoted. nes1983
errr. Благодарность? :) andrew cooke

Похожие вопросы