Гауссовский пример не работает с использованием Python

Я пытаюсь разместить в столбцах эту матрицу data с помощью этого Python код:

#!/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, и примерка работает хорошо.

Спасибо.


person ticofiz    schedule 08.02.2015    source источник


Ответы (1)


Вы столкнулись с классической проблемой предоставления неверного предположения для алгоритма подбора кривой. Это полностью из-за вашего ненужного переворачивания матрицы T вверх ногами, а затем без учета нового местоположения гауссианов (параметр с именем center, переданный в gaussian() - я помню этот код).

Видите ли, вот что происходит, когда я подхожу к вашим исходным данным:

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.

person Oliver W.    schedule 08.02.2015
comment
Большое Вам спасибо. Очень четкое и красивое объяснение! - person ticofiz; 09.02.2015