使用 lmfit 的二维高斯拟合
2D Gaussian fit using lmfit
我需要将二维高斯拟合到我读入的数据集。我选择的拟合例程是 lmfit,因为它可以轻松实现边界条件和固定参数。因为我不是最有效率的程序员,所以我在实现我的需求时遇到了问题。这是我所做的:
from numpy import *
from math import *
from lmfit import Parameters,minimize,report_fit
## fails to run
# from https://www.w3resource.com/python-exercises/numpy/python-numpy-exercise-79.php
x,y = meshgrid(linspace(-1,1,10),linspace(-1,1,10))
#d = sqrt(x*x+y*y)
#sigma, mu = 1.0, 0.0
#g = exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )
def gaussian2D(p,x,y):
height = p["height"].value
centroid_x = p["centroid_x"].value
centroid_y = p["centroid_y"].value
sigma_x = p["sigma_x"].value
sigma_y = p["sigma_y"].value
background = p["background"].value
return height*exp(-(((centroid_x-x)/sigma_x)**2+((centroid_y-y)/sigma_y)**2)/2.0)+background
def residuals(p,x,y,z):
return z - gaussian2D(p,x,y)
initial = Parameters()
initial.add("height",value=1.)
initial.add("centroid_x",value=0.)
initial.add("centroid_y",value=0.)
initial.add("sigma_x",value=1.)
initial.add("sigma_y",value=3.)
initial.add("background",value=0.)
xx,yy = meshgrid(x,y)
fit = minimize(residuals,initial,args=(array(xx).flatten(),array(yy).flatten(),array(g).flatten()))
popt = fit.params
print report_fit(fit)
首先,生成 2D 高斯的示例代码无法 运行 并给出 "TypeError: only size-1 arrays can be converted to Python scalars" for d = sqrt(xx+yy)。因为无论如何我都在使用文件中的数据,所以我正在处理网站 here.
上提供的示例数据
一些 research told me 将二维数组转换为一维数据,以便 lmfit 能够处理它们。我尝试在我的数组上使用 flatten 方法来实现它,但没有成功,给出了同样的错误(类型错误:只有大小为 1 的数组可以转换为 Python 标量)。我不够精通,无法完全理解 link.
中的代码
如果有任何帮助,我将不胜感激,尤其是。因为我更喜欢定义自己的函数来拟合数据,而不是依赖内置模型。
我认为你很接近,只是混淆了何时(或多久)调用 meshgrid
。修改后的版本是
import numpy as np
from lmfit import Parameters, minimize, report_fit
x, y = np.meshgrid(np.linspace(-1, 1, 10), np.linspace(-1, 1, 10))
def gaussian2D(x, y, cen_x, cen_y, sig_x, sig_y, offset):
return np.exp(-(((cen_x-x)/sig_x)**2 + ((cen_y-y)/sig_y)**2)/2.0) + offset
def residuals(p, x, y, z):
height = p["height"].value
cen_x = p["centroid_x"].value
cen_y = p["centroid_y"].value
sigma_x = p["sigma_x"].value
sigma_y = p["sigma_y"].value
offset = p["background"].value
return (z - height*gaussian2D(x,y, cen_x, cen_y, sigma_x, sigma_y, offset))
# test data
g = gaussian2D(x, y, 1.2, 2.1, 0.5, 0.7, 1.1)
initial = Parameters()
initial.add("height",value=1.)
initial.add("centroid_x",value=0.)
initial.add("centroid_y",value=0.)
initial.add("sigma_x",value=1.)
initial.add("sigma_y",value=3.)
initial.add("background",value=0.)
fit = minimize(residuals, initial, args=(x, y, g))
print(report_fit(fit))
也就是说,定义一个您可以更好地使用和测试的 gaussian2D()
函数,然后有一个简单的 objective 函数来调用它。
我需要将二维高斯拟合到我读入的数据集。我选择的拟合例程是 lmfit,因为它可以轻松实现边界条件和固定参数。因为我不是最有效率的程序员,所以我在实现我的需求时遇到了问题。这是我所做的:
from numpy import *
from math import *
from lmfit import Parameters,minimize,report_fit
## fails to run
# from https://www.w3resource.com/python-exercises/numpy/python-numpy-exercise-79.php
x,y = meshgrid(linspace(-1,1,10),linspace(-1,1,10))
#d = sqrt(x*x+y*y)
#sigma, mu = 1.0, 0.0
#g = exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )
def gaussian2D(p,x,y):
height = p["height"].value
centroid_x = p["centroid_x"].value
centroid_y = p["centroid_y"].value
sigma_x = p["sigma_x"].value
sigma_y = p["sigma_y"].value
background = p["background"].value
return height*exp(-(((centroid_x-x)/sigma_x)**2+((centroid_y-y)/sigma_y)**2)/2.0)+background
def residuals(p,x,y,z):
return z - gaussian2D(p,x,y)
initial = Parameters()
initial.add("height",value=1.)
initial.add("centroid_x",value=0.)
initial.add("centroid_y",value=0.)
initial.add("sigma_x",value=1.)
initial.add("sigma_y",value=3.)
initial.add("background",value=0.)
xx,yy = meshgrid(x,y)
fit = minimize(residuals,initial,args=(array(xx).flatten(),array(yy).flatten(),array(g).flatten()))
popt = fit.params
print report_fit(fit)
首先,生成 2D 高斯的示例代码无法 运行 并给出 "TypeError: only size-1 arrays can be converted to Python scalars" for d = sqrt(xx+yy)。因为无论如何我都在使用文件中的数据,所以我正在处理网站 here.
上提供的示例数据一些 research told me 将二维数组转换为一维数据,以便 lmfit 能够处理它们。我尝试在我的数组上使用 flatten 方法来实现它,但没有成功,给出了同样的错误(类型错误:只有大小为 1 的数组可以转换为 Python 标量)。我不够精通,无法完全理解 link.
中的代码如果有任何帮助,我将不胜感激,尤其是。因为我更喜欢定义自己的函数来拟合数据,而不是依赖内置模型。
我认为你很接近,只是混淆了何时(或多久)调用 meshgrid
。修改后的版本是
import numpy as np
from lmfit import Parameters, minimize, report_fit
x, y = np.meshgrid(np.linspace(-1, 1, 10), np.linspace(-1, 1, 10))
def gaussian2D(x, y, cen_x, cen_y, sig_x, sig_y, offset):
return np.exp(-(((cen_x-x)/sig_x)**2 + ((cen_y-y)/sig_y)**2)/2.0) + offset
def residuals(p, x, y, z):
height = p["height"].value
cen_x = p["centroid_x"].value
cen_y = p["centroid_y"].value
sigma_x = p["sigma_x"].value
sigma_y = p["sigma_y"].value
offset = p["background"].value
return (z - height*gaussian2D(x,y, cen_x, cen_y, sigma_x, sigma_y, offset))
# test data
g = gaussian2D(x, y, 1.2, 2.1, 0.5, 0.7, 1.1)
initial = Parameters()
initial.add("height",value=1.)
initial.add("centroid_x",value=0.)
initial.add("centroid_y",value=0.)
initial.add("sigma_x",value=1.)
initial.add("sigma_y",value=3.)
initial.add("background",value=0.)
fit = minimize(residuals, initial, args=(x, y, g))
print(report_fit(fit))
也就是说,定义一个您可以更好地使用和测试的 gaussian2D()
函数,然后有一个简单的 objective 函数来调用它。