利用Matlab实现直线和平面的拟合 
2011-04-14 10:45:43|  分类： 算法思想 |举报|字号 订阅
直线和平面拟合是很常用的两个算法，原理非常简单。但如果matlab不太熟的话，写起来也不是那么容易。搜了很久才到这两个代码，保存之，免得日后麻烦。
1、直线拟合的matlab代码
% Fitting a best-fit line to data, both noisy and non-noisy
x = rand(1,10);
n = rand(size(x)); % Noise
y = 2*x + 3; % x and y satisfy y = 2*x + 3
yn = y + n; % x and yn roughly satisfy yn = 2*x + 3 due to the noise
% Determine coefficients for non-noisy line y=m1*x+b1
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
Coeffs = [Xcolv Const]\Ycolv; % Find the coefficients
m1 = Coeffs(1);
b1 = Coeffs(2);
% To fit another function to this data, simply change the first
% matrix on the line defining Coeffs
% For example, this code would fit a quadratic
% y = Coeffs(1)*x^2+Coeffs(2)*x+Coeffs(3)
% Coeffs = [Xcolv.^2 Xcolv Const]\Ycolv;
% Note the .^ before the exponent of the first term
% Plot the original points and the fitted curve
figure
plot(x,y,'ro')
hold on
x2 = 0:0.01:1;
y2 = m1*x2+b1; % Evaluate fitted curve at many points
plot(x2, y2, 'g-')
title(sprintf('Non-noisy data: y=%f*x+%f',m1,b1))
% Determine coefficients for noisy line yn=m2*x+b2
Xcolv = x(:); % Make X a column vector
Yncolv = yn(:); % Make Yn a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
NoisyCoeffs = [Xcolv Const]\Yncolv; % Find the coefficients
m2 = NoisyCoeffs(1);
b2 = NoisyCoeffs(2);
% Plot the original points and the fitted curve
figure
plot(x,yn,'ro')
hold on
x2 = 0:0.01:1;
yn2 = m2*x2+b2;
plot(x2, yn2, 'g-')
title(sprintf('Noisy data: y=%f*x+%f',m2,b2))
2、平面拟合matlab代码
x = rand(1,10);
y = rand(1,10);
z = (3-2*x-5*y)/4; % Equation of the plane containing
% (x,y,z) points is 2*x+5*y+4*z=3
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
L=plot3(x,y,z,'ro'); % Plot the original data points
set(L,'Markersize',2*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','r') % Filling in the markers
hold on
[xx, yy]=meshgrid(0:0.1:1,0:0.1:1); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
% By rotating the surface, you can see that the points lie on the plane
% Also, if you multiply both sides of the equation in the title by 4,
matlab拟合数据
% you get the equation in the comment on the third line of this example
如何用matlab最小二乘法进行平面拟合
MATLAB软件提供了基本的曲线拟合函数的命令：
多项式函数拟合： a = polyfit（xdata，ydata，n）
其中n表示多项式的最高阶数，xdata，ydata 为要拟合的数据，它是用数组的方式输入。输出参数a为拟合多项式y = a1xn + … + anx + an+1的系数a = [a1, …, an, an+1]。
多项式在x处的值y可用下面程序计算。
y = polyval (a, x)
一般的曲线拟合： p = curvefit（‘Fun’p0，xdata，ydata）
其中Fun表示函数Fun (p, xdata)的M-文件，p0表示函数的初值。curvefit命令的求解问题形式是：
min{p} sum {(Fun (p, xdata)-ydata).^2}
若要求解点x处的函数值可用程序f = Fun(p, x) 计算。
例如已知函数形式 y = ae - bx + ce – dx ，并且已知数据点(xi, yi), i = 1,2,…, n，要确定四个未知参数a, b, c, d。
使用curvefit命令，数据输入xdata = [x1,x2, …, xn]; ydata = [y1,y2, …, yn];初值输入p0 = [a0,b0,c0,d0]; 并且建立函数y = ae - bx + ce – dx的M-文件（Fun.m）。若定义p1 = a, p2 = b, p3 = c, p4 = d , 则输出p = [p1, p2, p3, p4]。
引例求解：
  t=[1:16]; %数据输入
  y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6];
  plot(t,y,'o') %画散点图
  p=polyfit(t,y,2) (二次多项式拟合)
计算结果：
  p = -0.0445    1.0711    4.3252  %二次多项式的系数
从而得到某化合物的浓度y与时间t的拟合函数：
y = 4.3252+1.0711t – 0.0445t2
对函数的精度如何检测呢？仍然以图形来检测，将散点与拟合曲线画在一个画面上。
  xi=linspace(0,16,160);
  yi=polyval(p,xi);
  plot(x,y,'o',xi,yi)
  在MATLAB的NAG Foundation Toolbox中也有一些曲面拟合函数，如e02daf，e02cf，e02def可分别求出矩形网格点数据、散点数据的最小平方误差双三次样条曲面拟合，e02def等可求出曲面拟合的函数值。
用matlab的regress命令进行平面拟合
(2011-08-16 22:00:38)
转载▼
标签：
教育
分类： 数学软件
以少量数据为例
x = [1 5 6 3 7]';
y = [2 9 3 5 8]';
z = [4 3 5 11 6]';
scatter3(x,y,z,'filled')
hold on
即可将散点绘制出来


我们继续
 X = [ones(5,1) x y]; //5为size(x)
 b = regress(z,X) //拟合，其实是线性回归，但可以用来拟合平面。regress命令还有其它用
法，但一般这样就可以满足要求了。
 于是显示出
 b =

    6.5642
  -0.1269
  -0.0381
这就表示 z = 6.5643 - 0.1269 * x - 0.0381 * y 是拟合出来的平面的方程

下面把它绘制出来
xfit = min(x):0.1:max(x);  //注 0.1表示数据的间隔
yfit = min(y):0.1:max(y);
[XFIT,YFIT]= meshgrid (xfit,yfit); //制成网格数据
ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT;
mesh (XFIT,YFIT,ZFIT)
这样，图就出来啦

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。