最小二乘拟合平面和直线matlab

利用Matlab实现直线和平面的拟合1、直线拟合的matlab代码% Fitting a best-fit line to data, both noisy and non-noisyx = rand(1,10);n = rand(size(x)); % Noisey = 2*x + 3; % x and y satisfy y = 2*x + 3yn = y + n; % x and yn roughly satisfy yn = 2*x + 3 due to the noise % Determine coefficients for non-noisy line y=m1*x+b1Xcolv = x(:); % Make X a column vectorYcolv = y(:); % Make Y a column vectorConst = ones(size(Xcolv)); % Vector of ones for constant term Coeffs = [Xcolv Const]\Ycolv; % Find the coefficientsm1 = 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 curvefigureplot(x,y,'ro')hold onx2 = 0:0.01:1;y2 = m1*x2+b1; % Evaluate fitted curve at many pointsplot(x2, y2, 'g-')title(sprintf('Non-noisy data: y=%f*x+%f',m1,b1))% Determine coefficients for noisy line yn=m2*x+b2Xcolv = x(:); % Make X a column vectorYncolv = yn(:); % Make Yn a column vectorConst = ones(size(Xcolv)); % Vector of ones for constant term NoisyCoeffs = [Xcolv Const]\Yncolv; % Find the coefficientsm2 = NoisyCoeffs(1);b2 = NoisyCoeffs(2);% Plot the original points and the fitted curvefigureplot(x,yn,'ro')hold onx2 = 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=3Xcolv = x(:); % Make X a column vectorYcolv = y(:); % Make Y a column vectorZcolv = z(:); % Make Z a column vectorConst = ones(size(Xcolv)); % Vector of ones for constant termCoefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficientsXCoeff = Coefficients(1); % X coefficientYCoeff = Coefficients(2); % X coefficientCCoeff = Coefficients(3); % constant term% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeffL=plot3(x,y,z,'ro'); % Plot the original data pointsset(L,'Markersize',2*get(L,'Markersize')) % Making the circle markers larger set(L,'Markerfacecolor','r') % Filling in the markershold 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 surfacetitle(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,% 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 = a1x n+ … + a n x + a n+1的系数a = [a1, …, a n, a n+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命令进行平面拟合以少量数据为例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)这样，图就出来啦%r p q就是你的x,y,zr = randi(10,20,1);p = randi(10,20,1);q = randi(10,20,1);b = regress(r,[p q]);scatter3(r,p,q,'filled');hold onrfit = min(r):1:max(r);pfit = min(p):1:max(p);[RFIT PFIT] = meshgrid(rfit,pfit);QFIT = b(1) * RFIT + b(2) * PFIT;mesh(RFIT,PFIT,QFIT);view(60,10);这个程序有问题，只能拟合得到Z＝AX+BY，得不到Z＝AX+BY+D的形式%点X,Y,Z到平面Ax+By+Cz+D=0的距离为%d(ABCD,XYZ)=|AX+BY+CZ+D|/sqrt(A^2+B^2+C^2)%ABCD四个变量只有三个是互相独立的%设A=cos(a);B=sin(a)*cos(b);C=sin(a)*sin(b)%那么A^2+B^2+C^2=1,距离公式化简为%d(abc,XYZ)=|cos(a)*X+sin(a)*cos(b)*Y+sin(a)*sin(b)*Z+c|%现在有已知点序列X，Y，Z，求参数a b c%先构造一个函数fun(p) 输入参数为p，其中p(1)=a,p(2)=b,p(3)=c%使用lsqnonlin求得p,使得sum((fun(p))^2)最小fun=@(p) cos(p(1))*X+sin(p(1))*cos(p(2))*Y+sin(p(1))*sin(p(2))*Z+p(3);p = lsqnonlin(fun,[0 0 0]);A=cos(p(1));B=sin(p(1))*cos(p(2));C=sin(p(1))*sin(p(2));D=p(3);。