2017-10-21 201 views
1

我相信我成功實現了使用來自scipy的曲線擬合的高斯擬合。但是我遇到的問題是...適合性不是很好,因爲優化的參數正在改變質心。使用*** curve_fit ***從scipy和python 3.x改進高斯擬合

data =np.loadtxt('mock.txt') 
    my_x=data[:,0] 
    my_y=data[:,1] 

    def gauss(x,mu,sigma,A): 
     return A*np.exp(-(x-mu)**2/2/sigma**2) 
    def trimodal_gauss(x,mu1,sigma1,A1,mu2,sigma2,A2,mu3,sigma3,A3): 
     return gauss(x,mu1,sigma1,A1)+gauss(x,mu2,sigma2,A2)+gauss(x,mu3,sigma3,A3) 



    """"" 
    Gaussian fitting parameters recognized in each file 
    """"" 
    first_centroid=(10180.4*2+9)/9 
    second_centroid=(10180.4*2+(58.6934*1)+7)/9 
    third_centroid=(10180.4*2+(58.6934*2)+5)/9 
    centroid=[] 
    centroid+=(first_centroid,second_centroid,third_centroid) 

    apparent_resolving_power=1200 
    sigma=[] 
    for i in range(len(centroid)): 
     sigma.append(centroid[i]/((apparent_resolving_power)*2.355)) 

    height=[1,1,1] 

    p=[]  

    p = np.array([list(t) for t in zip(centroid, sigma, height)]).flatten() 


    popt, pcov = curve_fit(trimodal_gauss,my_x,my_y,p0=p) 

輸出:enter image description here

據我所知,在這裏有很多山峯,但我真的需要它,以適應只有三個高斯,但在正確的質心(在我最初的猜測給出)。換句話說,我真的不希望我給出的質心不變。有沒有人遇到過這樣的挑戰?並且可以請我幫助我做些什麼來實現目標?

+0

一般情況下,似乎你最好配合正確的號碼(至少5個,也許6個),並且只採用你實際關心的三個結果。你目前的方法會做不好的工作,因爲你不適合的高峯會影響你關心的三個高峯的結果。它會「認爲」額外的東西是三個峯值的一部分。 – NichtJens

回答

0

如果使用lmfit模塊(https://github.com/lmfit/lmfit-py),你可以很容易地把邊界上的高斯函數的重心,甚至解決這些問題。 Lmfit還可以輕鬆建立多峯模型。

你沒有給一個完整的例子或鏈接到你的數據,但lmfit到您的數據的擬合可能是這樣的:

import numpy as np 
from lmfit import GaussianModel 
data =np.loadtxt('mock.txt') 
my_x=data[:,0] 
my_y=data[:,1] 

model = (GaussianModel(prefix='p1_') + 
      GaussianModel(prefix='p2_') + 
      GaussianModel(prefix='p3_')) 

params = model.make_params(p1_amplitude=100, p1_sigma=2, p1_center=2262, 
          p2_amplitude=100, p2_sigma=2, p2_center=2269, 
          p3_amplitude=100, p3_sigma=2, p3_center=2276, 
         ) 

# set boundaries on the Gaussian Centers: 
params['p1_center'].min = 2260 
params['p1_center'].max = 2264 

params['p2_center'].min = 2267 
params['p2_center'].max = 2273 

params['p3_center'].min = 2274 
params['p3_center'].max = 2279 

# or you could just fix one of the centroids like this: 
params['p3_center'].vary = False 

# if needed, you could force all the sigmas to be the same value 
# or related by simple mathematical expressions 
params['p2_sigma'].expr = 'p1_sigma' 
params['p3_sigma'].expr = '2*p1_sigma' 

# fit this model to data: 
result = model.fit(my_y, params, x=my_x) 

# print results 
print(result.fit_report()) 

# evaluate individual gaussian components: 
peaks = model.eval_components(params=result.params, x=my_x) 

# plot results: 
plt.plot(my_x, my_y, label='data') 
plt.plot(my_x, result.best_fit, label='best fit') 
plt.plot(my_x, peaks['p1_']) 
plt.plot(my_x, peaks['p2_']) 
plt.plot(my_x, peaks['p3_']) 
plt.show() 
+0

'scipy.curve_fit'知道當前版本的界限。看到這裏:https://stackoverflow.com/questions/16760788/python-curve-fit-library-that-allows-me-to-assign-bounds-to-parameters – NichtJens

0

您應該爲中心定義三個固定值的單獨函數。然後你只適合剩餘參數的這些函數的求和函數。

簡而言之,您的trimodal_gauss()不應採取mu s,但只有A s和sigma s。 mu應該是常數。

這樣做的一個微不足道的(但不是很普遍)的方法是:

def trimodal_gauss(x, sigma1, A1, sigma2, A2, sigma3, A3): 
    mu1 = 1234 # put your mu's here 
    mu2 = 2345 
    mu3 = 3456 
    g1 = gauss(x, mu1, sigma1, A1) 
    g2 = gauss(x, mu2, sigma2, A2) 
    g3 = gauss(x, mu3, sigma3, A3) 
    return g1 + g2 + g3 

從這個人們可以通過一個「發電機」爲trimodal_gauss函數需要三個概括的想法mu(或n?) s並創建其他參數的功能。像這樣:

def make_trimodal_gauss(mu1, mu2, mu3): 
    def result_function(x, sigma1, A1, sigma2, A2, sigma3, A3): 
     g1 = gauss(x, mu1, sigma1, A1) 
     g2 = gauss(x, mu2, sigma2, A2) 
     g3 = gauss(x, mu3, sigma3, A3) 
     return g1 + g2 + g3 
    return result_function 


mu1 = 1234 # put your mu's here 
mu2 = 2345 
mu3 = 3456 

trimodal_gauss = make_trimodal_gauss(mu1, mu2, mu3) 

#usage like this: trimodal_gauss(x, sigma1, A1, sigma2, A2, sigma3, A3) 
+0

感謝您的評論。如果mu不是函數的一部分,我們如何確保高斯函數適合於感興趣的質心? – user7852656

+0

我添加了一個方法來做到我的答案。 – NichtJens