पायथन: गैर-रैखिक कम-वर्गों के साथ दो-वक्र गाऊशियन फिटिंग

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

इस कोड को मुझे प्रदान है कि आप केवल एक समारोह है कि एक फिटिंग कर रहे हैं के लिए काम किया दो गॉसियन वितरण का संयोजन।

मैंने अभी एक अवशिष्ट कार्य किया है जो दो गॉसियन कार्यों को जोड़ता है और फिर उन्हें वास्तविक डेटा से घटा देता है।

पैरामीटर (पी) जो मैंने न्यूम्पी के कम से कम वर्ग समारोह में पारित किया है उनमें शामिल हैं: पहले गॉसियन फ़ंक्शन (एम) का अर्थ, पहले और दूसरे गॉसियन कार्यों (डीएम, यानी क्षैतिज शिफ्ट) से माध्य में अंतर , पहले (एसडी 1) का मानक विचलन, और दूसरे (एसडी 2) के मानक विचलन।

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 पतित कर रहे हैं - वहाँ बिल्कुल डेटा है कि उन दोनों के बीच तय कर सकते हैं में कुछ नहीं है। आपको दो के बजाय एक शून्य स्तर पैरामीटर का उपयोग करना चाहिए (यानी अपने कोड से उनमें से एक को हटा दें)। यह शायद आपके फिट को रोक रहा है (यहां टिप्पणियों को अनदेखा करें कि यह संभव नहीं है - उस डेटा में कम से कम दो चोटियों को स्पष्ट रूप से स्पष्ट किया गया है और आपको निश्चित रूप से उसमें फिट होना चाहिए)।

(यह स्पष्ट नहीं हो सकता है कि मैं इसका सुझाव क्यों दे रहा हूं, लेकिन क्या हो रहा है कि coeffs 0 और 4 एक-दूसरे को रद्द कर सकते हैं। वे दोनों शून्य हो सकते हैं, या कोई 100 हो सकता है और दूसरा -100 - या तो वैसे, फिट उतना ही अच्छा है। यह फिटिंग दिनचर्या "भ्रमित" करता है, जो अपना समय व्यतीत करने की कोशिश करता है कि उन्हें क्या करना चाहिए, जब कोई भी सही उत्तर न हो, क्योंकि जो भी मूल्य है, दूसरा हो सकता है उस के नकारात्मक, और फिट वही होगा)।

वास्तव में, साजिश से, ऐसा लगता है कि शून्य स्तर की आवश्यकता नहीं हो सकती है। मैं उन दोनों को छोड़ने और फिट दिखने के तरीके को देखने की कोशिश करता हूं।

भी, कम से कम वर्गों में coeffs 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