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

利用Matlab实现直线和平面的拟合1、直线拟合的matlab代码%Fittingabest-fitlinetodata,bothnoisyandnon-noisyx=rand(1,10);n=rand(size(x));%Noisey=2*x+3;%xandysatisfyy=2*x+3yn=y+n;%xandynroughlysatisfyyn=2*x+3duetothenoise%Determinecoefficientsfornon-noisyliney=m1*x+b1Xcolv=x(:);%MakeXacolumnvectorYcolv=y(:);%MakeYacolumnvectorConst=ones(size(Xcolv));%VectorofonesforconstanttermCoeffs=[XcolvConst]\Ycolv;%Findthecoefficientsm1=Coeffs(1);b1=Coeffs(2);%Tofitanotherfunctiontothisdata,simplychangethefirst%matrixonthelinedefiningCoeffs%Forexample,thiscodewouldfitaquadratic%y=Coeffs(1)*x^2+Coeffs(2)*x+Coeffs(3)%Coeffs=[Xcolv.^2XcolvConst]\Ycolv;%Notethe.^beforetheexponentofthefirstterm%Plottheoriginalpointsandthefittedcurvefigureplot(x,y,'ro')holdonx2=0:0.01:1;y2=m1*x2+b1;%Evaluatefittedcurveatmanypointsplot(x2,y2,'g-')title(sprintf('Non-noisydata:y=%f*x+%f',m1,b1))%Determinecoefficientsfornoisylineyn=m2*x+b2Xcolv=x(:);%MakeXacolumnvectorYncolv=yn(:);%MakeYnacolumnvectorConst=ones(size(Xcolv));%VectorofonesforconstanttermNoisyCoeffs=[XcolvConst]\Yncolv;%Findthecoefficientsm2=NoisyCoeffs(1);b2=NoisyCoeffs(2);%Plottheoriginalpointsandthefittedcurvefigureplot(x,yn,'ro')holdonx2=0:0.01:1;yn2=m2*x2+b2;plot(x2,yn2,'g-')title(sprintf('Noisydata:y=%f*x+%f',m2,b2))2、平面拟合matlab代码x=rand(1,10);y=rand(1,10);z=(3-2*x-5*y)/4;%Equationoftheplanecontaining%(x,y,z)pointsis2*x+5*y+4*z=3Xcolv=x(:);%MakeXacolumnvectorYcolv=y(:);%MakeYacolumnvectorZcolv=z(:);%MakeZacolumnvectorConst=ones(size(Xcolv));%VectorofonesforconstanttermCoefficients=[XcolvYcolvConst]\Zcolv;%FindthecoefficientsXCoeff=Coefficients(1);%XcoefficientYCoeff=Coefficients(2);%XcoefficientCCoeff=Coefficients(3);%constantterm%Usingtheabovevariables,z=XCoeff*x+YCoeff*y+CCoeffL=plot3(x,y,z,'ro');%Plottheoriginaldatapointsset(L,'Markersize',2*get(L,'Markersize'))%Makingthecirclemarkerslargerset(L,'Markerfacecolor','r')%Fillinginthemarkersholdon[xx,yy]=meshgrid(0:0.1:1,0:0.1:1);%Generatingaregulargridforplottingzz=XCoeff*xx+YCoeff*yy+CCoeff;surf(xx,yy,zz)%Plottingthesurfacetitle(sprintf('Plottingplanez=(%f)*x+(%f)*y+(%f)',XCoeff,YCoeff,CCoeff))%Byrotatingthesurface,youcanseethatthepointslieontheplane%Also,ifyoumultiplybothsidesoftheequationinthetitleby4,%yougettheequationinthecommentonthethirdlineofthisexample如何用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=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];plot(t,y,'o')%画散点图p=polyfit(t,y,2)(二次多项式拟合)计算结果：p=-0.04451.07114.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的NAGFoundationToolbox中也有一些曲面拟合函数，如e02daf，e02cf，e02def可分别求出矩形网格点数据、散点数据的最小平方误差双三次样条曲面拟合，e02def等可求出曲面拟合的函数值。用matlab的regress命令进行平面拟合以少量数据为例x=[15637]';y=[29358]';z=[435116]';scatter3(x,y,z,'filled')holdon即可将散点绘制出来我们继续X=[ones(5,1)xy];//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)这样，图就出来啦%rpq就是你的x,y,zr=randi(10,20,1);p=randi(10,20,1);q=randi(10,20,1);b=regress(r,[pq]);scatter3(r,p,q,'filled');holdonrfit=min(r):1:max(r);pfit=min(p):1:max(p);[RFITPFIT]=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，求参数abc%先构造一个函数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,[000]);A=cos(p(1));B=sin(p(1))*cos(p(2));C=sin(p(1))*sin(p(2));D=p(3);