高斯拟合无法使用 Python
gaussian fitting not working using Python
我正在尝试使用此 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 并且拟合效果很好。
谢谢。
您遇到了向曲线拟合算法提供 错误猜测 的经典问题。这完全是由于您不必要地颠倒矩阵 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.
等技巧使您的算法更健壮
我正在尝试使用此 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 并且拟合效果很好。
谢谢。
您遇到了向曲线拟合算法提供 错误猜测 的经典问题。这完全是由于您不必要地颠倒矩阵 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.