高斯擬合工作不使用Python

Question

我試圖用這個Python代碼，以適應列這個矩陣data：

#!/usr/local/bin/env python 
import numpy as np 
import Tkinter #Used for file import 
import tkFileDialog #Used for file import 
import os 
import scipy 
import scipy.optimize as optimize 


root = Tkinter.Tk() 
root.withdraw() #use to hide tkinter window 

filename = os.getcwd() 
background = os.getcwd() 
filename = tkFileDialog.askopenfile(parent=root,mode='rb',title='Choose a file') 
background = tkFileDialog.askopenfile(parent=root,mode='rb',title='Choose a background') 
filename = filename.name 
#filename = r'bb1e03' 
#background = r'bb1e03_background' 


T0 = np.loadtxt(filename, unpack=False) 
bg = np.loadtxt(background, unpack=False) 

T = T0-bg # background subtraction? 
#T = T.clip(min=0) 
T[T<0]=0 
T = np.flipud(T) 
N, M = T.shape 
datax = np.arange(N) 


def gaussian(x, height, center, width, offset): 
    return height*np.exp(-(x - center)**2/(2*width**2)) + offset 

def three_gaussians(x, h1, c1, w1, h2, c2, w2, h3, c3, w3, offset): 
    return (gaussian(x, h1, c1, w1, offset=0) + 
     gaussian(x, h2, c2, w2, offset=0) + 
     gaussian(x, h3, c3, w3, offset=0) + offset) 

def two_gaussians(x, h1, c1, w1, h2, c2, w2, offset): 
    return three_gaussians(x, h1, c1, w1, h2, c2, w2, 0,0,1, offset) 

def one_gaussian(x,h1,c1,w1, offset): 
    return (gaussian(x, h1, c1, w1, offset=0)+offset) 

#errfunc3 = lambda p, x, y: (three_gaussians(x, *p) - y)**2 
#errfunc2 = lambda p, x, y: (two_gaussians(x, *p) - y)**2 
#errfunc1 = lambda p, x, y: (one_gaussian(x, *p) - y)**2 
#output files for fit parameters 
outfile1 = open('results_1gau.txt', 'w') 
outfile2 = open('results_2gau.txt', 'w') 
outfile3 = open('results_3gau.txt', 'w') 
outfile1.write('column\th1\tc1\tw1\toffset\n') 
outfile2.write('column\th1\tc1\tw1\th2\tc2\tw2\toffset\n') 
outfile3.write('column\th1\tc1\tw1\th2\tc2\tw2\th3\tc3\tw3\toffset\n') 

# new matrices for fitted data 
datafit1 = np.empty_like(T) 
datafit2 = np.empty_like(T) 
datafit3 = np.empty_like(T) 

for n in xrange(M): 

    Mmax = T[:,n].max() 
    guess1 = [0.5*Mmax, N/10., 10., 0.] 
    guess2 = [0.5*Mmax, N/10., 10., 0.5*Mmax, N/10., 10., 0.] 
    guess3 = [0.5*Mmax, N/10., 10., 0.5*Mmax, N/10., 10., 
       0.5*Mmax, N/10., 10., 0] 
    #optim3, success = optimize.leastsq(errfunc3, guess3[:], 
    #         args=(datax, data[:,n])) 
    #optim2, success = optimize.leastsq(errfunc2, guess2[:], 
    #         args=(datax, data[:,n])) 

    try: 
     optim1, pcov = optimize.curve_fit(one_gaussian, datax, T[:,n], guess1) 
    except: 
     optim1 = [0, 0, 1, 0] 

    try: 
     optim2, pcov = optimize.curve_fit(two_gaussians, datax, T[:,n], guess2) 
    except: 
     optim2 = [0, 0, 1, 0, 0, 1, 0] 

    try: 
     optim3, pcov = optimize.curve_fit(three_gaussians, datax, T[:,n], guess3) 
    except: 
     optim3 = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] 

    # write parameters to file (1 gau) 
    s = '{}'.format(n) 
    for x in guess1: 
     s += '\t{:g}'.format(x) 
    outfile1.write(s + '\n') 

    # write parameters to file (2 gau) 
    s = '{}'.format(n) 
    for x in guess2: 
     s += '\t{:g}'.format(x) 
    outfile2.write(s + '\n') 

    # write parameters to file (3 gau) 
    s = '{}'.format(n) 
    for x in guess3: 
     s += '\t{:g}'.format(x) 
    outfile3.write(s + '\n') 
    # fill new matrices with fitted data 
    datafit1[:,n] = one_gaussian(datax, *optim1) 
    datafit2[:,n] = two_gaussians(datax, *optim2) 
    datafit3[:,n] = three_gaussians(datax, *optim3) 

T = datafit1 

我看過最相關配件相關的職位，但我找不到什麼是錯的用我的代碼。它應該可以工作，但是最終的矩陣「T」只顯示具有常數的列，而不是平滑的高斯形狀曲線。請看看並告訴我我做錯了什麼。我曾在其他程序中嘗試過，例如OriginLab和配件很好。

謝謝。

Answer 1

您遇到了向曲線擬合算法提供錯誤猜測的經典問題。這完全是由於你不必要的顛倒矩陣T，然後不考慮高斯的新位置（參數center，傳遞給gaussian()--我記得this code）。

你看，這裏的時候我做配件上自己的原始數據會發生什麼：

T = T0-bg # background subtraction? 
fitparams_me, fitparams_you = [], [] 

for colind in xrange(16,19): 
    column = T[:,colind] 
    guess = column.max(), column.argmax(), 3, 0 # Good guess for a SINGLE gaussian 

    popt, pcov = optimize.curve_fit(one_gaussian, datax, column, p0=guess) 
    fitparams_me.append(popt) 
print(fitparams_me) 

其中顯示：

[array([ 365.40098996, 91.24095009, 1.11390434, -0.99632476]), 
array([ 348.4327168 , 92.0262556 , 1.26650618, -1.08018819]), 
array([ 413.21526868, 90.8569241 , 1.0445618 , -1.0565371 ])] 

而這些導致非常好的擬合。

現在，您正在做的是：您首先倒置矩陣，但您一直假設峯位於第一行。然而，這不再是這種情況，此代碼亮點：

T = np.flipud(T) 

for colind in xrange(16,19): 
    column = T[:,colind] 
    guess = column.max(), column.argmax(), 3, 0 # Good guess for a SINGLE gaussian 
    your_guess = [0.5*Mmax, N/10., 10., 0.] 
    print guess[1], your_guess[1] 
    popt, pcov = optimize.curve_fit(one_gaussian, datax, column, p0=your_guess) 
    fitparams_you.append(popt) 

# printed results: 
932 102.4 
931 102.4 
932 102.4 

所以，每次我還在猜測正確，其中最大的出現，但你假設它始終圍繞你的數據的102行（形狀爲1024, 1024）。僅通過翻轉你的專欄容易

>>> print(fitparams_you) 

[array([ -1.640e-07, 1.024e+02, 1.000e+01, 2.046e-10]), 
array([ -1.640e-07, 1.024e+02, 1.000e+01, 2.046e-10]), 
array([ -1.640e-07, 1.024e+02, 1.000e+01, 2.046e-10])] 

你可以解決這個問題：

popt, pcov = optimize.curve_fit(one_gaussian, datax, column[::-1], p0=your_guess) 

或者你

它應是毫不奇怪，從曲線擬合的結果差別很大，從礦井可以嘗試使用像argmax這樣的技巧使您的算法更加健壯。