2013-08-31 342 views
12

我有a,b点的离散规则网格和它们相应的c值,我对它进一步插值得到一条平滑曲线。现在从插值数据,我还想创建一个曲线拟合的多项式方程。如何用多项式拟合三维图?三维曲线拟合

我尝试在MATLAB中做到这一点。我在MATLAB(r2010a)中使用了曲面拟合工具箱来曲线拟合三维数据。但是,如何在MATLAB/MAPLE或任何其他软件中找到适合一组数据的公式,以达到最佳效果。有什么建议?也是最有用的将是一些真正的代码示例看看,PDF文件,在网络上等

这只是我的数据的一小部分。

a = [ 0.001 .. 0.011]; 

b = [1, .. 10]; 

c = [ -.304860225, .. .379710865]; 

在此先感谢。

回答

15

要将曲线拟合到一组点上,我们可以使用ordinary least-squares回归。 MathWorks有一个描述过程的solution page

举个例子,让我们用一些随机数据开始:

% some 3d points 
data = mvnrnd([0 0 0], [1 -0.5 0.8; -0.5 1.1 0; 0.8 0 1], 50); 

由于@BasSwinckels表明,通过构建所需的design matrix,您可以使用mldividepinvsolve the overdetermined system表示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)) 

1st_order_polynomial


正如提到的,我们可以通过添加更多方面的独立变量矩阵(在AAx=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在统计工具箱,有助于构建设计矩阵的一对夫妇模型订单:

C = x2fx(data(:,1:2), 'quadratic') \ data(:,3); 
zz = x2fx([xx(:) yy(:)], 'quadratic') * C; 
zz = reshape(zz, size(xx)); 

最后有一个很好的功能polyfitn由John D'Errico在文件交换中允许您指定所涉及的各种多项式次序和项:

model = polyfitn(data(:,1:2), data(:,3), 2); 
zz = polyvaln(model, [xx(:) yy(:)]); 
zz = reshape(zz, size(xx)); 

2nd_order_polynomial

+0

如何在python中执行同样的操作..!?位帮助将升值... @Amro – diffracteD

+2

@diffracteD:我将代码翻译成Python:https://gist.github.com/amroamroamro/1db8d69b4b65e8bc66a6 – Amro

+0

感谢您的帮助...我一定会尝试它.. ! – diffracteD

3

有可能是对文件交换一些更好的功能,但手工做这件事是这样的:

x = a(:); %make column vectors 
y = b(:); 
z = c(:); 

%first order fit 
M = [ones(size(x)), x, y]; 
k1 = M\z; 
%least square solution of z = M * k1, so z = k1(1) + k1(2) * x + k1(3) * y 

同样,你可以做一个二阶拟合:

%second order fit 
M = [ones(size(x)), x, y, x.^2, x.*y, y.^2]; 
k2 = M\z; 

这似乎对您提供的有限数据集有数值问题。请输入help mldivide了解更多详情。

要对一些规则的网格内插,就可以像这样:

ngrid = 20; 
[A,B] = meshgrid(linspace(min(a), max(a), ngrid), ... 
       linspace(min(b), max(b), ngrid)); 
M = [ones(numel(A),1), A(:), B(:), A(:).^2, A(:).*B(:), B(:).^2]; 
C2_fit = reshape(M * k2, size(A)); % = k2(1) + k2(2)*A + k2(3)*B + k2(4)*A.^2 + ... 

%plot to compare fit with original data 
surfl(A,B,C2_fit);shading flat;colormap gray 
hold on 
plot3(a,b,c, '.r') 

三阶拟合可以使用以下TryHard给出的公式来完成,但是计算公式很快变得单调乏味时,订单增加。如果你不得不这样做,那么最好编写一个函数,该函数可以构造M给定x,yorder

+1

+1三阶拟合与M = [ones(size(x)),x,y,x。^ 2,x。* y,y。^ 2,x。^ 3,x 。^ 2. * y,x。* y。^ 2,y。^ 3]'做得不错... –

+0

对不起,迟到回复! Thankyou回复Bas Swinckels并努力尝试。好的解决方案但是,是否有任何方法可以通过代码或使用MATLAB中的任何工具找到最合适的公式/公式?如何找到最适合给定数据集的公式? – Syeda

+0

我也收到了这个警告。警告:排名不足,排名= 5,tol = 9.9961e-013。你能帮我理解这是什么意思吗?谢谢。 – Syeda

2

这听起来更像是一个哲学问题,而不是具体的实现,特别是位 - “如何找到一个适合一组数据的公式来获得最佳优势?”在我的经验,这是一个选择你必须取决于你想要实现的。

什么为你定义“最好”?对于数据拟合问题,您可以继续添加越来越多的多项式系数,并获得更好的R^2值......但最终会“过度拟合”数据。高阶多项式的缺点是行为超出了用于适应响应表面的样本数据范围之外的行为 - 它可能会以某种狂野的方向快速消失,这可能不适合您尝试建模的任何数据。

您是否了解您正在使用的系统/数据的物理行为?这可以用作用于创建数学模型的一组方程式的基础。我的建议是使用您可以摆脱的最经济(简单)的模型。