多项式拟合 matlab，对系数有一些约束

Question

我有数据，我应该使用必须属于以下类型的函数进行插值：

f(x) = ax<sup>4</sup> + bx<sup>2</sup> + c

与a > 0 和b ≤ 0。不幸的是，MATLAB 的polyfit 不允许对多项式的系数进行任何约束。有人知道是否有 MATLAB 函数可以做到这一点吗？不然怎么实现？

非常感谢您，

伊丽莎白

Answer 1

您可以尝试使用fminsearch、fminunc 手动定义您的目标函数。

或者，您可以稍微不同地定义您的问题：

f(x) = a<sup>2</sup>x<sup>4</sup> - b<sup>2</sup>x<sup>2</sup> + c

现在，新的a 和b 可以针对无约束进行优化，同时确保您正在寻找的最终a 和b 是正数（负数分别为负数）。 )。

Answer 2

没有约束，这个问题可以写成一个简单的线性系统来解决：

% Your design matrix ([4 2 0] are the powers of the polynomial)
A = bsxfun(@power, your_X_data(:), [4 2 0]);

% Best estimate for the coefficients, [a b c], found by 
% solving A*[a b c]' = y in a least-squares sense
abc = A\your_Y_data(:)

这些约束当然会自动满足iff约束模型确实是您的数据的基础。例如，

% some example factors
a = +23.9;
b = -15.75;
c = 4;

% Your model
f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3);

% generate some noisy XY data
x = -1:0.01:1;
y = f(x, [a b c]) + randn(size(x));

% Best unconstrained estimate a, b and c from the data
A = bsxfun(@power, x(:), [4 2 0]);
abc = A\y(:);

% Plot results
plot(x,y, 'b'), hold on
plot(x, f(x, abc), 'r')
xlabel('x (nodes)'), ylabel('y (data)')

但是，如果您对受约束模型未准确描述的数据施加约束，事情可能会出错：

% Note: same data, but flipped signs 
a = -23.9;
b = +15.75;
c = 4;

f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3);

% generate some noisy XY data
x = -1:0.01:1;
y = f(x, [a b c]) + randn(size(x));

% Estimate a, b and c from the data, Forcing a>0 and b<0
abc = fmincon(@(Y) sum((f(x,Y)-y).^2), [0 0 0], [-1 0 0; 0 +1 0; 0 0 0], zeros(3,1));

% Plot results
plot(x,y, 'b'), hold on
plot(x, f(x, abc), 'r')
xlabel('x (nodes)'), ylabel('y (data)')

（此解决方案有a == 0，表示模型选择不正确）。

如果a == 0 的完全相等是一个问题：如果你设置a == eps(0) 当然没有区别。从数值上讲，这对于真实世界的数据不会很明显，但它仍然是非零的。

无论如何，我怀疑您的模型选择不当，并且约束是使一切正常工作的“修复”，或者您的数据实际上应该在尝试进行任何拟合之前是无偏/重新调整的，或者类似的先决条件适用（我经常看到人们做这种事情，所以是的，我在这方面有点偏见:)。

那么...这些限制背后的真正原因是什么？

Answer 3

如果您有曲线拟合工具箱，那么fit 确实允许使用“上”和“下”选项设置约束。你会想要类似的东西。

M=fit(x, f, 'poly4', 'upper', [-inf, 0, -inf, 0, -inf], 'lower', [0, 0, 0, 0, -inf]);

注意使用 -inf 将特定系数设置为不受约束。

这将给出一个具有相关系数的 cfit 对象。您可以使用例如 M.p1 来访问这些 x^4 术语。或者，您可以使用 feval 在任何您想要的点评估函数。

我认为您也可以使用优化工具箱中的 lsqcurvefit 来做类似的事情。