洛伦兹拟合的Matlab函数与全局变量答案

【问题标题】：Matlab function for lorentzian fit with global variables洛伦兹拟合的Matlab函数与全局变量
【发布时间】：2019-04-06 03:32:56
【问题描述】：

我想将洛伦兹拟合到我的数据中，所以首先我想测试我对模拟数据的拟合过程：

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 页

【问题讨论】：

附带一点：不要使用global 变量将数据传递给您的函数（即X 和Y），而是使用匿名函数。请参阅fminsearch docs 中的示例 2
a0 的值是多少，它们与[20,30,20] 有很大不同吗？还有什么不是零，即d 的实际价值是多少？它不应该完全为零，因为您的 Y 包含一些噪音，但它应该非常小
但大多数情况下，您应该检查d = devsum(afinal)
[112.3817 50.0000 100.0000] 是我对a0 的值，与[20, 30 ,20 ] 相差甚远。我最担心的是，无论 a0 是什么向量，我的函数都会以某种方式返回零。当我计算我在函数中声明的内容时，我得到了正确的结果。我还尝试了您的建议并创建了一个函数d=devsum(a,X,Y) 然后我的函数给出了正确的结果，因此这些全局变量不知何故有问题。但是 afinal = fminsearch(@(a0) devsum(a0,X,Y)) 失败了。顺便说一句，我不希望这样：解决我的错误，但我想知道出了什么问题。
所以，首先a0 实际上并没有像你最初的方式那样被改变。 afinal 是更改（优化）的值。所以devsum(afinal) 应该等于0，NOT devsum(a0)

标签： matlab curve-fitting

【解决方案1】：

通常最好避免使用全局变量。我通常解决这些问题的方法是首先定义一个函数，该函数将您想要拟合的曲线评估为x 的函数和参数：

% lorentz.m
function y = lorentz(param, x)
y = param(1) ./ ((x-param(2)).^2 + param(3))

这样，您以后可以重用该函数来绘制拟合结果。

然后，您定义一个带有您想要最小化的属性的小型匿名函数，只有一个参数作为输入，因为这是fminsearch 需要的格式。在匿名函数的定义中，测量的X 和Y 不是使用全局变量，而是“捕获”（技术术语是对这些变量执行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')

结果：

【讨论】：