Python: doğrusal olmayan en küçük kareler ile iki eğri Gauss uydurma

Question

Matematik bilgim sınırlıdır, bu yüzden muhtemelen sıkışmış durumdayım. İki Gauss zirvesine uymaya çalıştığım bir spektrumum var. En büyük zirveye uyabilirim, ama en küçük zirveye sığamam. Anladığım kadarıyla, iki zirve için Gauss işlevini toplamalıyım ama nerede yanlış gittiğimi bilmiyorum. benim şimdiki çıkışı bir görüntüsü gösterilmiştir:

mavi çizgi verilerim ve yeşil hat benim şimdiki seçimdir. Şu anda aşağıdaki kodu kullanarak, sığdırıyorum benim verilerindeki ana zirvenin solunda bir omuz vardır: Bu kod bana sadece bir fonksiyonu uydurma gösterilmesi koşulu çalıştı

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

iki Gauss dağılımının kombinasyonu.

Sadece iki Gauss işlevini ekleyen ve sonra bunları gerçek verilerden çıkartan bir artık işlevi yapıyorum.

Numpy'nin en küçük kareler fonksiyonuna geçtiğim parametreler (p) şunları içerir: ilk Gauss fonksiyonunun ortalaması (m), birinci ve ikinci Gauss fonksiyonlarından ortalama farkı (dm, yani yatay kayma) , birinci (sd1) standart sapması ve ikinci (sd2) standart sapması. Aralarında karar verebilir verilerde kesinlikle hiçbir şey yoktur -

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 ve 4 dejenere bulunmaktadır. İki yerine tek bir sıfır seviyesi parametresi kullanmalısınız (yani kodunuzdan birini kaldırabilirsiniz). bu muhtemelen uyumunuzu durduran şeydir (burada bunun mümkün olmadığını söyleyen yorumları göz ardı edin - bu verilerde en az iki zirve vardır ve kesinlikle buna uymalısınız).

(Neden bunu önerdiğim açık olmayabilir, fakat ne oluyor? 0 ve 4 katsayıları birbirini iptal edebilir, her ikisi de sıfır olabilir veya biri 100 ve diğeri 100 olabilir. yol, bu kadar iyi bir şeydir: bu, "ne zaman bir değer, ne olursa olsun, bir tek doğru cevap olmadığı zaman, ne olması gerektiğine çalışmak için harcayan zaman harcıyor" uydurma rutini karıştırır, diğeri sadece bunun negatif ve fit aynı olacaktır. Hiç bir sıfır düzeyine gerek olabilir gibi

aslında, arsa gelen görünüyor. İkisini de düşürmeyi ve zevkin nasıl göründüğünü görmeyi deneyebilirim. Ayrıca, en küçük karelerde 1 ve 5 katsayısını (veya sıfır noktasını) sığdırmaya gerek yoktur. Bunun yerine, model lineer olduğundan, her bir döngü değerini hesaplayabilirsiniz. Bu, işleri daha hızlı hale getirecek, ancak kritik değil. Sadece matematiklerinin çok iyi olmadığını söylediğini fark ettim, muhtemelen bunu görmezden gel.

Answer 3

Sen scikit-learn gelen Gauss karışım modelleri kullanabilirsiniz:

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

Ayrıca ncomp parametresi ile istediğiniz Gauss sayısını sığdırmak için aşağıdaki fonksiyonu kullanabilirsiniz:

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