用于全局变量的洛伦兹拟合的 Matlab 函数
Matlab function for lorentzian fit with global variables
我想用洛伦兹拟合我的数据,所以首先我想测试我对模拟数据的拟合过程:
X = linspace(0,100,200);
Y = 20./((X-30).^2+20)+0.08*randn(size(X));
起始参数
a3 = ((max(X)-min(X))/10)^2;
a2 = (max(X)+min(X))/2;
a1 = max(Y)*a3;
a0 = [a1,a2,a3];
找到适合的最小值
afinal = fminsearch(@devsum,a0);
afinal 是带有适合我的参数的向量。如果我按如下方式测试我的功能
d= devsum(a0)
然后 d= 0,但如果我完全按照函数中的内容进行操作
a=a0;
d = sum((Y - a(1)./((X-a(2)).^2+a(3))).^2)
那么d不等于0。这怎么可能?我的功能非常简单,所以我不知道出了什么问题。
我的功能:
%devsum.m
function d = devsum(a)
global X Y
d = sum((Y - a(1)./((X-a(2)).^2+a(3))).^2);
end
基本上我只是在实现我在这里找到的东西
http://www.home.uni-osnabrueck.de/kbetzler/notes/fitp.pdf
第 7 页
通常最好避免使用全局变量。我通常解决这些问题的方法是首先定义一个函数,它评估你想要拟合的曲线作为 x
和参数的函数:
% lorentz.m
function y = lorentz(param, x)
y = param(1) ./ ((x-param(2)).^2 + param(3))
这样,您可以在以后重用该函数来绘制拟合结果。
然后,您使用要最小化的 属性 定义一个小的匿名函数,只有一个参数作为输入,因为这是 fminsearch
需要的格式。在匿名函数的定义中,测量的 X
和 Y
不是使用全局变量,而是 'captured'(技术术语是对这些变量执行 closure):
fit_error = @(param) sum((y_meas - lorentz(param, x_meas)).^2)
最后,通过使用 fminsearch
:
最小化误差来拟合参数
fitted_param = fminsearch(fit_error, starting_param);
快速演示:
% simulate some data
X = linspace(0,100,200);
Y = 20./((X-30).^2+20)+0.08*randn(size(X));
% rough guess of initial parameters
a3 = ((max(X)-min(X))/10)^2;
a2 = (max(X)+min(X))/2;
a1 = max(Y)*a3;
a0 = [a1,a2,a3];
% define lorentz inline, instead of in a separate file
lorentz = @(param, x) param(1) ./ ((x-param(2)).^2 + param(3));
% define objective function, this captures X and Y
fit_error = @(param) sum((Y - lorentz(param, X)).^2);
% do the fit
a_fit = fminsearch(fit_error, a0);
% quick plot
x_grid = linspace(min(X), max(X), 1000); % fine grid for interpolation
plot(X, Y, '.', x_grid, lorentz(a_fit, x_grid), 'r')
legend('Measurement', 'Fit')
title(sprintf('a1_fit = %g, a2_fit = %g, a3_fit = %g', ...
a_fit(1), a_fit(2), a_fit(3)), 'interpreter', 'none')
结果:
我想用洛伦兹拟合我的数据,所以首先我想测试我对模拟数据的拟合过程:
X = linspace(0,100,200);
Y = 20./((X-30).^2+20)+0.08*randn(size(X));
起始参数
a3 = ((max(X)-min(X))/10)^2;
a2 = (max(X)+min(X))/2;
a1 = max(Y)*a3;
a0 = [a1,a2,a3];
找到适合的最小值
afinal = fminsearch(@devsum,a0);
afinal 是带有适合我的参数的向量。如果我按如下方式测试我的功能
d= devsum(a0)
然后 d= 0,但如果我完全按照函数中的内容进行操作
a=a0;
d = sum((Y - a(1)./((X-a(2)).^2+a(3))).^2)
那么d不等于0。这怎么可能?我的功能非常简单,所以我不知道出了什么问题。 我的功能:
%devsum.m
function d = devsum(a)
global X Y
d = sum((Y - a(1)./((X-a(2)).^2+a(3))).^2);
end
基本上我只是在实现我在这里找到的东西 http://www.home.uni-osnabrueck.de/kbetzler/notes/fitp.pdf 第 7 页
通常最好避免使用全局变量。我通常解决这些问题的方法是首先定义一个函数,它评估你想要拟合的曲线作为 x
和参数的函数:
% lorentz.m
function y = lorentz(param, x)
y = param(1) ./ ((x-param(2)).^2 + param(3))
这样,您可以在以后重用该函数来绘制拟合结果。
然后,您使用要最小化的 属性 定义一个小的匿名函数,只有一个参数作为输入,因为这是 fminsearch
需要的格式。在匿名函数的定义中,测量的 X
和 Y
不是使用全局变量,而是 'captured'(技术术语是对这些变量执行 closure):
fit_error = @(param) sum((y_meas - lorentz(param, x_meas)).^2)
最后,通过使用 fminsearch
:
fitted_param = fminsearch(fit_error, starting_param);
快速演示:
% simulate some data
X = linspace(0,100,200);
Y = 20./((X-30).^2+20)+0.08*randn(size(X));
% rough guess of initial parameters
a3 = ((max(X)-min(X))/10)^2;
a2 = (max(X)+min(X))/2;
a1 = max(Y)*a3;
a0 = [a1,a2,a3];
% define lorentz inline, instead of in a separate file
lorentz = @(param, x) param(1) ./ ((x-param(2)).^2 + param(3));
% define objective function, this captures X and Y
fit_error = @(param) sum((Y - lorentz(param, X)).^2);
% do the fit
a_fit = fminsearch(fit_error, a0);
% quick plot
x_grid = linspace(min(X), max(X), 1000); % fine grid for interpolation
plot(X, Y, '.', x_grid, lorentz(a_fit, x_grid), 'r')
legend('Measurement', 'Fit')
title(sprintf('a1_fit = %g, a2_fit = %g, a3_fit = %g', ...
a_fit(1), a_fit(2), a_fit(3)), 'interpreter', 'none')
结果: