要将曲线拟合到一组点上,我们可以使用普通最小二乘 https://en.wikipedia.org/wiki/Ordinary_least_squares回归。有一个解决方案页面 http://www.mathworks.com/support/solutions/en/data/1-1AVW5/MathWorks 描述了该过程。
作为示例,让我们从一些随机数据开始:
% some 3d points
data = mvnrnd([0 0 0], [1 -0.5 0.8; -0.5 1.1 0; 0.8 0 1], 50);
As @BasSwinckels https://stackoverflow.com/a/18552769/97160表明,通过构建所需的设计矩阵 https://en.wikipedia.org/wiki/Design_matrix, 您可以使用mldivide http://www.mathworks.com/help/matlab/ref/mldivide.html or pinv http://www.mathworks.com/help/matlab/ref/pinv.html to 解决超定系统 http://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html#f4-2064表示为Ax=b
:
% best-fit plane
C = [data(:,1) data(:,2) ones(size(data,1),1)] \ data(:,3); % coefficients
% evaluate it on a regular grid covering the domain of the data
[xx,yy] = meshgrid(-3:.5:3, -3:.5:3);
zz = C(1)*xx + C(2)*yy + C(3);
% or expressed using matrix/vector product
%zz = reshape([xx(:) yy(:) ones(numel(xx),1)] * C, size(xx));
接下来我们可视化结果:
% plot points and surface
figure('Renderer','opengl')
line(data(:,1), data(:,2), data(:,3), 'LineStyle','none', ...
'Marker','.', 'MarkerSize',25, 'Color','r')
surface(xx, yy, zz, ...
'FaceColor','interp', 'EdgeColor','b', 'FaceAlpha',0.2)
grid on; axis tight equal;
view(9,9);
xlabel x; ylabel y; zlabel z;
colormap(cool(64))
如前所述,我们可以通过向自变量矩阵添加更多项来获得高阶多项式拟合(A
in Ax=b
).
假设我们想要拟合一个具有常数项、线性项、交互项和平方项 (1, x, y, xy, x^2, y^2) 的二次模型。我们可以手动执行此操作:
% best-fit quadratic curve
C = [ones(50,1) data(:,1:2) prod(data(:,1:2),2) data(:,1:2).^2] \ data(:,3);
zz = [ones(numel(xx),1) xx(:) yy(:) xx(:).*yy(:) xx(:).^2 yy(:).^2] * C;
zz = reshape(zz, size(xx));
还有一个辅助功能x2fx http://www.mathworks.com/help/stats/x2fx.html在统计工具箱中,有助于构建几个模型阶数的设计矩阵:
C = x2fx(data(:,1:2), 'quadratic') \ data(:,3);
zz = x2fx([xx(:) yy(:)], 'quadratic') * C;
zz = reshape(zz, size(xx));
终于有一个很棒的功能了polyfitn http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitnJohn D'Errico 的文件交换允许您指定所涉及的各种多项式阶数和项:
model = polyfitn(data(:,1:2), data(:,3), 2);
zz = polyvaln(model, [xx(:) yy(:)]);
zz = reshape(zz, size(xx));