Python：使用非線性最小二乘法的雙曲線高斯擬合

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

此代碼爲我工作提供，你只擬合函數是一個兩個高斯分佈的組合。

我只是做了一個殘差函數，它添加了兩個高斯函數，然後從真實數據中減去它們。

我傳遞給Numpy最小二乘函數的參數（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() 

Answer 2

coeffs 0和4退化 - 是絕對沒有的，可以決定它們之間的數據。你應該使用一個零水平參數而不是兩個（即從你的代碼中刪除其中的一個）。這可能是阻止你的合適（忽略這裏的評論，說這是不可能的 - 這些數據中至少有兩個高峯，你當然應該能夠適應這一點）。

（可能不太清楚爲什麼我提出這個建議，但是發生的事情是係數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