最小二乘法圆拟合及matlab程序
最小二乘法曲线拟合-原理及matlab实现
曲线拟合（curve-fitting ）：工程实践中，用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ϕ来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势（即使)(x ϕ最好地逼近()x f ，而不必满足插值原则。
因此没必要取)(i x ϕ=i y ，只要使i i i y x -=)(ϕδ尽可能地小）。
原理：给定数据点},...2,1,0),,{(m i y x i i =。
求近似曲线)(x ϕ。
并且使得近似曲线与()x f 的偏差最小。
近似曲线在该点处的偏差i i i y x -=)(ϕδ，i=1,2,...,m 。
常见的曲线拟合方法:1.使偏差绝对值之和最小2.使偏差绝对值最大的最小3.使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程：1. 设拟合多项式为：kk x a x a a x +++=...)(10ϕ2. 各点到这条曲线的距离之和，即偏差平方和如下：3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数，因而我们得到了：.......4、 把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5. 将这个范德蒙得矩阵化简后可得到:6. 也就是说X*A=Y ，那么A = (X'*X)-1*X'*Y ，便得到了系数矩阵A ，同时，我们也就得到了拟合曲线。
MATLAB实现：MATLAB提供了polyfit（）函数命令进行最小二乘曲线拟合。
调用格式：p=polyfit(x,y,n)[p,s]= polyfit(x,y,n)[p,s,mu]=polyfit(x,y,n)x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。
x 必须是单调的。
矩阵s包括R（对x进行QR分解的三角元素）、df(自由度）、normr(残差）用于生成预测值的误差估计。
最小二乘法曲线拟合的Matlab程序
方便大家使用的最小二乘法曲线拟合的Matlab程序非常方便用户使用,直接按提示操作即可;这里我演示一个例子:(红色部分为用户输入部分,其余为程序运行的结果,结果图为Untitled.fig,Untitled2.fig) 请以向量的形式输入x,y.x=[1,2,3,4]y=[3,4,5,6]通过下面的交互式图形，你可以事先估计一下你要拟合的多项式的阶数，方便下面的计算.polytool()是交互式函数,在图形上方[Degree]框中输入阶数,右击左下角的[Export]输出图形回车打开polytool交互式界面回车继续进行拟合输入多项式拟合的阶数m = 4Warning: Polynomial is not unique; degree >= number of data points. > In polyfit at 72In zxecf at 64输出多项式的各项系数a = 0.0200000000000001a = -0.2000000000000008a = 0.7000000000000022a = 0.0000000000000000a = 2.4799999999999973输出多项式的有关信息 SR: [4x5 double]df: 0normr: 2.3915e-015Warning: Zero degrees of freedom implies infinite error bounds.> In polyval at 104In polyconf at 92In zxecf at 69观测数据拟合数据x y yh1.0000 3.0000 3.00002.0000 4.0000 4.00003 5 54.0000 6.0000 6.0000剩余平方和 Q = 0.000000标准误差 Sigma = 0.000000相关指数 RR = 1.000000请输入你所需要拟合的数据点，若没有请按回车键结束程序.输入插值点x0 = 3输出插值点拟合函数值 y0 = 5.0000>>结果：untitled.figuntitled2.fig一些matlab优化算法代码的分享代码的目录如下：欢迎讨论1.约束优化问题：minRosen(Rosen梯度法求解约束多维函数的极值)(算法还有bug) minPF(外点罚函数法解线性等式约束)minGeneralPF(外点罚函数法解一般等式约束)minNF(内点罚函数法)minMixFun(混合罚函数法)minJSMixFun(混合罚函数加速法)minFactor(乘子法)minconPS(坐标轮换法)(算法还有bug)minconSimpSearch(复合形法)2.非线性最小二乘优化问题minMGN(修正G-N法)3.线性规划：CmpSimpleMthd(完整单纯形法)4.整数规划(含0-1规划)DividePlane(割平面法)ZeroOneprog(枚举法)5.二次规划QuadLagR(拉格朗日法)ActivedeSet(起作用集法)6.辅助函数(在一些函数中会调用)minNT(牛顿法求多元函数的极值)Funval(求目标函数的值)minMNT(修正的牛顿法求多元函数极值)minHJ(黄金分割法求一维函数的极值)7.高级优化算法1）粒子群优化算法(求解无约束优化问题)1>PSO(基本粒子群算法)2>YSPSO(待压缩因子的粒子群算法)3>LinWPSO(线性递减权重粒子群优化算法)4>SAPSO(自适应权重粒子群优化算法)5>RandWSPO(随机权重粒子群优化算法)6>LnCPSO(同步变化的学习因子)7>AsyLnCPSO(异步变化的学习因子)(算法还有bug)8>SecPSO(用二阶粒子群优化算法求解无约束优化问题)9>SecVibratPSO(用二阶振荡粒子群优化算法求解五约束优化问题)10>CLSPSO(用混沌群粒子优化算法求解无约束优化问题)11>SelPSO(基于选择的粒子群优化算法)12>BreedPSO(基于交叉遗传的粒子群优化算法)13>SimuAPSO(基于模拟退火的粒子群优化算法)2)遗传算法1>myGA(基本遗传算法解决一维约束规划问题)2>SBOGA(顺序选择遗传算法求解一维无约束优化问题)3>NormFitGA(动态线性标定适应值的遗传算法求解一维无约束优化问题)4>GMGA(大变异遗传算法求解一维无约束优化问题)5>AdapGA(自适应遗传算法求解一维无约束优化问题)6>DblGEGA(双切点遗传算法求解一维无约束优化问题)7>MMAdapGA(多变异位自适应遗传算法求解一维无约束优化问题)自己编写的马尔科夫链程序A 代表一组数据序列一维数组本程序的操作对象也是如此t=length(A); % 计算序列“A”的总状态数B=unique(A); % 序列“A”的独立状态数顺序，“E”E=sort(B,'ascend');a=0;b=0;c=0;d=0;for j=1:1:ttLocalization=find(A==E(j)); % 序列“A”中找到其独立状态“E”的位置for i=1:1:length(Localization)if Localization(i)+1>tbreak; % 范围限定elseif A(Localization(i)+1)== E(1)a=a+1;elseif A(Localization(i)+1)== E(2)b=b+1;elseif A(Localization(i)+1)== E(3)c=c+1;% 依此类推，取决于独立状态“E”的个数elsed=d+1;endendT(j,1:tt)=[a,b,c,d]; % “T”为占位矩阵endTT=T;for u=2:1:ttTT(u,:)= T(u,:)- T(u-1,:);endTT; % 至此，得到转移频数矩阵Y=sum(TT,2);for uu=1:1:ttTR(uu,:)= TT(uu,:)./Y(uu,1);endTR % 最终得到马尔科夫转移频率/概率矩阵% 观测序列马尔科夫性质的检验：N=numel(TT);uuu=1;Col=sum(TT,2); % 对列求和Row=sum(TT,1); % 对行求和Total=sum(Row); % 频数总和for i=1:1:ttfor j=1:1:ttxx(uuu,1)=sum((TT(i,j)-(Row(i)*Col(j))./Total).^2./( (Row(i)*Col(j)). /Total));uuu=uuu+1; % 计算统计量x2endendxx=sum(xx)。
最小二乘法、最佳均方逼近、随机拟合及其MATLAB程序
2曲线拟合的线性最小二乘法及其MATLAB?序例2给出一组数据点(X i, y i)列入表2中，试用线性最小二乘法求拟合曲线, 估计其误差，作出拟合曲线•解 (1)在MATLAB工作窗口输入程序>> x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.126.50 68.04];plot(x,y , 'r*'),lege nd( '实验数据(xi,yi)' )xlabel( 'x' ), ylabel( 'y'),title( '数据点(xi,yi) 的散点图’)运行后屏幕显示数据的散点图(略)(3)编写下列MATLAB程序计算f(x)在(X j,yj处的函数值，即输入程序>> syms al a2 a3 a4x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];fi=a1.*x.A3+ a2.*x.A2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi =[ -125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4,a4, 1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4]编写构造误差平方和的MATLAB程序>> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04];fi=[-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4, a4,1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4];fy=fi-y； fy2=fy.A2; J=sum(fy.A2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)A2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)A2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20F2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25F2+(a4+91/10F2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)A2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)A2+(19683/1000 *a1+729/100*a2+27/10*a3+a4-13/2F2+(5832/125*a1+324/25*a2 +18/5*a3+a4-1701/25F2为求31,32,33,34使J达到最小，只需利用极值的必要条件-丄 0 a k (k 1,2,3,4)，得到关于31,32,33,34的线性方程组，这可以由下面的MATLAB程序完成，即输入程序>> syms a1 a2 a3 a4 J=(-125/8*a1+25/4*32-5/2*a3+34+1929/10)A2+(-4913/1000*a1+289/100*32-17/10*33+34...+171/2)A2+(-1331/1000*a1+1 21/100*a2-11/10*a3+34+723/20)A2+(-64/125*31+16/25*32-4/5*a3+34+663/25)A2+(34+91/10)A2+(1/1000*31+1/100*32+1/10*a3+a4+843/100)A2+(27/8*31+9/4*32+3/2*a3+34+328/25)A2+(19683/ 1000*a1+729/100*32+27/10*a3+34-13/2)A2+(5832/125*31+324/2 5*a2+18/5*a3+a4-1701/25)A2;Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3);Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),Ja41=simple(Ja4),运行后屏幕显示J分别对31, 32 ,33 ,34的偏导数如下Ja1仁56918107/10000*31+32097579/25000*32+1377283/2500*33+23667/250*34-8442429/625J321 =32097579/25000*31+1377283/2500*32+23667/250*33 +67*34+767319/625 J331 = 1377283/2500*31+23667/250*32+67*33+18/5*34-232638/125J341 = 23667/250*31+67*32+18/5*33+18*34+14859/25解线性方程组J311 =0，J321 =0，J331 =0，J341 =0，输入下列程序>>A=[56918107/10000, 32097579/25000, 1377283/2500,23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18];B=[8442429/625, -767319/625, 232638/125, -14859/25];C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*xA3-7988544102557579/562949953421312*xA2+1804307491277693/281474976710656*x -4648521160813215/562949953421312 故所求的拟合曲线为f (x) 5.0911 x314.1905 x2 6.4102 x 8.2574 .(4)编写下面的MATLAB 程序估计其误差，并作出拟合曲线和数据的图形.输入程序>> xi=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04]; n=length(xi);f=5.0911.*xi.A3-14.1905.*xi.A2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6;F=5.0911.*x.A3-14.1905.*x.A2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.A2; Ew=max(fy),E1=sum(fy)/n, E2=sqrt((sum(fy2))/n) plot(xi,y,'r*'), hold on, plot(x,F, 'b-'),hold off legend('数据点(xi,yi)','拟合曲线y=f(x)'),xlabel('x'), ylabel('y'),title(例2的数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据(X i，yj与拟合函数f的最大误差平均误差E i和均方根误差E2及其数据点(X j,yj和拟合曲线y=f(x)的图形(略).Ew = E1 = E2 =3.105 4 0.903 4 1.240 96函数逼近及其MATLAB?序最佳均方逼近的MATLAB^程序function [yy1,a,WE]=zjjfbj(f,X,Y,xx) m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1);if n~=length(Y) error( 'X和丫的维数应该相同') end for j=1:m for k=1:mb(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i))*feval(f(k,:),X(i));endendc(j)=0;for i=1:nc(j)=c(j)+feval(f(j,:),X(i))*Y(i);endenda=b\c;WE=0;for i=1:nff=0;for j=1:m ff=ff+a(j)*feval(f(j,:),X(i)); endWE=WE+(Y(i)-ff)*(Y(i)-ff);endif nargin==3return ;endyy=[];for i=1:ml=[];for j=1:length(xx) l=[l,feval(f(i,:),xx(j))]; end yy=[yy l'];endyy=yy*a; yy1=yy'; a=a';WE;例6. 1对数据X和Y,用函数y 1,y x, y x2进行逼近，用所得到的逼近函数计算在x 6.5处的函数值，并估计误差.其中X=(1 3 4 5 6 7 8 9); Y=(-11 -13 -11 -7 -1 7 17 29). 解在MATLA工作窗口输入程序>> X=[ 1 3 4 5 6 7 8 9]; Y=[-11 -13 -11 -7 -17 17 29];f=['fun0';'fun1';'fun2'];[yy,a,WE]=zjjfbj(f,X,Y,6.5) 运行后屏幕显示如下yy =2.75000000000003a =-7.00000000000010 -4.999999999999951.00000000000000WE =7.172323350269439e-027例 6.2 对数据X 和丫，用函数 y 1, y x, y x2，y cosx，y e x，y sinx进行逼近，其中X=(0 0.50 1.00 1.50 2.00 2.50 3.00 )，丫=(0 0.47940.8415 0.9815 0.9126 0.5985 0.1645 ) .解在MATLA工作窗口输入程序>> X=[ 0 0.50 1.00 1.50 2.00 2.50 3.00];丫=[0 0.4794 0.8415 0.9815 0.9126 0.1645];f=['fun0';'fun1';'fun2';'fun3';'fun4';'fun5'];xx=0:0.2:3;[yy,a,WE]=zjjfbj(f,X,Y,plot(X,Y,'ro',xx,yy,'b-')运行后屏幕显示如下(图略)yy = Columns 1 through 7-0.0005 0.2037 0.3939 0.5656 0.8348 0.9236Columns 8 through 140.9771 0.9926 0.9691 0.9069 0.6766 0.5191Columns 15 through 160.3444 0.1642 0.5985xx), 0.7141 0.8080a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653 WE = 1.5769e-004 即，最佳逼近函数为y=0.3828+0.4070*x-0.3901*xA2+0.0765*exp(x) +0.5653*sin(x) .8随机数据点上的二元拟合及其MATLA 程序例 8 设节点 ( X,Y,Z ) 中的 X 和 Y 分别是在区间 [ 3,3] 和 [ 2.5,3.5]上的 5022个 随 机 数 ， Z 是 函 数 Z=7-3 x 3e -x2 -y2 在 (X,Y ) 的 值 ， 拟 合 点 ( X I ,Y I ) 中的 X I =-3:0.2:3, Y I =-2.5:0.2:3.5. 分别用二元拟合方法中最近邻内插法、三 角基线性内插法、三角基三次内插法和 MATLAB4 网格化坐标方法计算在 ( X I ,Y I ) 处的值，作出它们的图形，并与被拟和曲面进行比较 .解 (1 )最近邻内插法 .输入程序>> x=rand(50,1); y=rand(50,1); %生成50个一元均匀分布随机数x 和y , x , y .X=-3+(3-(-3))*x; %利用x 生成的随机变量.title( ' 用最近邻内插法拟合函数 z =7-3 xA3 exp(-xA2 - yA2) 的曲面和节点的图形 ' )%legend( ' 拟合曲面 ',' 节点 (xi,yi,zi)' )hold on%在当前图形上添加新图形面及其插值乙(略).(2)三角基线性内插法 输入程序>> x=rand(50,1); y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x; %利用x 生成 上的随机变量.Y=-2.5+(3.5-(-2.5))*y;%利用 y 生成 上的随机变量 .Z=7-3* X.A3 .* exp(-X.A2 - Y.A2);%在每个随机点( X,Y )处计算Z 的值.-0.4598*cos(x)Y=-2.5+(3.5-(-2.5))*y;Z=7-3* X.A3 .* exp(-X.A2 - Y.A2);处计算Z 的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); ZI=griddata(X,Y,Z,XI,YI, (XI,YI )处的插值 ZI.mesh(XI,YI, ZI)xlabel( 'x' ), ylabel( 'y'%利用 y 生成的随机变量 .%在每个随机点( X,Y )%将坐标( XI,YI )网格化 .'nearest' ) %计算在每个插值点 %作二元拟合图形 . ), zlabel( 'z' ),plot3(X,Y,Z, 'bo' ) (X,Y,Z). hold of 运行后屏幕显示用最近邻内插法拟合函数%用兰色小圆圈画出每个节点 %结束在当前图形上添加新图形 .22Z=7-3 x 3e-x2 -y2在两组不同节点处的曲X1=-3:0.2:3;title('用三角基线性内插法拟合函数z =7-3 x A3 exp(-x A2 -y A2)的曲面和节点的图形’) %legend( ' 拟合曲面',' hold on plot3(X,Y,Z, 'bo' ) (X,Y,Z). hold of22运行后屏幕显示用三角基线性内插法拟合函数Z=7-3 x 3e-x -y在两组不同节点处的曲面和节点的图形及其插值 乙(略).(3)三角基三次内插法 . 输入程序>> x=rand(50,1); y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .title( ' 用三角基三次内插法拟合函数 z =7-3 xA3 exp(-xA2 -22运行后屏幕显示用三角基三次内插法拟合函数 Z=7-3 x 3e -x -y 在两组不同节点处的曲面和节点的图形及其插值 Z I (略).( 4 ) MATLAB 4网格化坐标方法 . 输入程序>> x=rand(50,1); y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成 上的随机变量.Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); ZI=griddata(X,Y,Z,XI,YI,XI,YI )处的插值 ZI.mesh(XI,YI, ZI)xlabel( 'x' ), ylabel( 'y%将坐标( XI,YI )网格化 .'linear' ) %计算在每个插值点 %作二元拟合图形 . ), zlabel( 'z' ),节点 (xi,yi,zi)' )%在当前图形上添加新图形 .%用兰色小圆圈画出每个节点 %结束在当前图形上添加新图形 . X=-3+(3-(-3))*x; %利用x 生成上的随机变量.Y=-2.5+(3.5-(-2.5))*y;Z=7-3* X.A3 .* exp(-X.A2 - Y.A2);处计算Z 的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); ZI=griddata(X,Y,Z,XI,YI,( XI,YI )处的插值 ZI.mesh(XI,YI, ZI) xlabel( 'x' ), ylabel( 'y'%利用y 生成上的随机变量•%在每个随机点( X,Y )%将坐标( XI,YI )网格化 .'cubic' ) %计算在每个插值点 %作二元拟合图形 . ), zlabel( 'z' ),yA2) 的曲面和节点的图形 ' )%legend( ' 拟合曲面','hold on 节点 (xi,yi,zi)' )%在当前图形上添加新图形plot3(X,Y,Z,'bo' ) (X,Y,Z).hold of%用兰色小圆圈画出每个节点%结束在当前图形上添加新图形 .22运行后屏幕显示用MATLAB 网格化坐标方法拟合函数Z=7-3 x 3e -x-y 在两组不同 节点处的曲面和节点的图形及其插值 ZI (略).22(5) 作被拟合曲面Z=7-3x 3e -x-y 和节点的图形. 输入程序>> x=ra nd(50,1); y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x; %利用x 生成随机变量.Y=-2.5+(3.5-(-2.5))*y;%利用y 生成随机变量.Z=7-3* X.A3 .* exp(-X.A2 - Y.A2);%在每个随机点( X,Y )处计算Z 的值.X1=-3.:0.1:3.;Y1=-2.5:0.1:3.5;[XI,YI] = meshgrid(X1,Y1);%将坐标( XI,YI )网格化 .ZI=7-3* XI.A3 .* exp(-XI.A2 - YI.A2); mesh(XI,YI, ZI) %作二元拟合图形 .xlabel( 'x'), ylabel('y' ), zlabel( 'z' ),title( ' 被拟合函数 z =7-3 xA3 exp(-xA2 - yA2) 的曲面和节点的图形 ' )%legend('被拟合函数曲面','节点(xi,yi,zi)' )hold on%在当前图形上添加新图形 .plot3(X,Y,Z, 'bo' )%用兰色小圆圈画出每个节点 (X,Y,Z).hold of%结束在当前图形上添加新图形 .22运行后屏幕显示被拟合函数 Z=7-3 x 3e -x-y 的曲面和节点的图形及其函数值 ZI(略) .Y=-2.5+(3.5-(-2.5))*y;Z=7-3* X.A3 .* exp(-X.A2 - Y.A2);处计算Z 的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5; [XI,YI] = meshgrid(X1,Y1); ZI=griddata(X,Y,Z,XI,YI, (XI,YI )处的插值 ZI.mesh(XI,YI, ZI) xlabel( 'x' ), ylabel( 'y' %利用y 生成上的随机变量.%在每个随机点( X,Y )'v4'%将坐标( XI,YI )网格化 . ) %计算在每个插值点%作二元拟合图形 . ), zlabel( 'z' ),' 用 MATLAB 4 网格化坐标方法 拟合函数 z =7-3 xA3 的曲面和节点的图形 ' )%legend( ' 拟合曲面 ',' 节点 (xi,yi,zi)' ) hold on %在当前图形上添加新图形 . plot3(X,Y,Z, 'bo' ) %用兰色小圆圈画出每个节点(X,Y,Z).hold oftitle(exp(-x A2 - y A2)'bo' ) %结束在当前图形上添加新图形 .。
最小二乘法matlab程序
最小二乘法matlab程序最小二乘法是一种统计模型，它可以被用来拟合一元函数数据，或者拟合非线性曲线。
它的基本思想是找到一组参数，使得拟合的曲线与实际数据的差距最小。
本文将介绍如何使用Matlab实现一个最小二乘法的程序，并与现有的一些现成的最小二乘法的matlab程序进行比较，找出其优缺点。
首先，要使用最小二乘法拟合曲线，需要准备一组输入数据，一般可以将其表示为两个向量，分别是自变量x和因变量y。
这些数据可以是由测量和实验得到的，也可以是由人工输入的，但无论如何都要确保它们的准确性。
接下来，就可以使用Matlab输入数据进行处理，用最小二乘法计算出最拟合的曲线及其参数。
具体步骤主要分为三步：第一步是计算输入数据的均值和方差，包括自变量x和因变量y的均值和方差；第二步是计算自变量x和因变量y的关系，即最小二乘拟合曲线的系数；第三步是验证拟合的曲线的准确性，如果不满意，可以重新调整参数，以获得较好的拟合效果。
此外，Matlab除了提供自带的最小二乘法函数外，还支持第三方开发者开发现成的matlab程序，用于解决最小二乘法的问题。
这些程序中有一些是开源的，另一些则是出售的。
其中开源的有LEAST，CURVEFIT，CURVEFITTOOL等，而出售的有MATLAB Curve Fitting Toolbox，Optimization Toolbox和Statistics Toolbox等。
它们的突出特点是速度快，代码简洁，容易上手，适用于多种拟合类型。
然而，各种matlab程序也有自身的缺点，最明显的就是当输入数据非常庞大时，它们的计算能力就无法跟上，速度就会变慢。
此外，使用出售的matlab程序可能相对昂贵，而且有时需要安装某些复杂的库文件，这也是一种麻烦。
因此，使用最小二乘法拟合曲线时，可以参考现有的matlab程序，也可以自己编写matlab代码，同时要考虑到程序的可靠性、效率和可行性。
本文介绍的matlab程序的最大优势是它不需要依赖第三方的软件，而且能够满足大多数用户的需求，使得最小二乘法可以在短时间内被成功运用。
最小二乘法圆拟合及matlab程序
Q(a,b, c)
a
2( X i2 Yi2 aX i bYi c) X i 0
①
Q(a,b, c)
b
2( X i2 Yi2 aX i bYi c)Yi 0 ②
Q(a,b, c)
c
2( X i2 Yi2 aX i bYi c) 0 ③
最小二乘法圆拟合
1
最小二乘法拟合圆曲线： R2 (x A)2 ( y B)2
R2 x2 2Ax A2 y2 2By B2
令a=-2A,b=-2B, c A2 B2 R2
则：圆的另一形式为：
x2 y2 ax by c 0
2
A a 只需求出参数a,b,c即可以求的圆半径参数： 2
A a 2
B b 2
R 1 a2 b2 4c 2
9
t=0:0.01:pi; a=20;%设定圆心X轴数值 b=30;%设定圆心Y轴数值 r=5;%设定圆半径数值 x=a+r*cos(t)+randn(1,315); y=b+r*sin(t)+randn(1,315); plot(x,y); hold on; x=x(:); y=y(:); m=[x y ones(size(x))]\[-(x.^2+y.^2)]; xc = -.5*m(1)%拟合圆心X轴数值 yc = -.5*m(2)%拟合圆心Y轴数值 R = sqrt((m(1)^2+m(2)^2)/4-m(3))%拟合半径数值 plot(xc,yc,'r-x',(xc+R*cos(t)),(yc+R*sin(t)),'r-'); axis equal;
最小二乘法圆拟合及matlab程序
X i2 Yi2 +aX i bYi c
3
令Q(a,b,c)为
的平方和：
i
Q(a, b, c) i2 [( Xi2 Yi2 aXi bYi c)]2
下面求参数a,b,c使得Q(a,b,c)的值最小即可
4
F(a,b,c)对a,b,c求偏导，令偏导等于0，得到极值点，比较所有极值点的函 数值即可得到最小值。
② × N- ③ × Yi
且令 C (N Xi2 Xi Xi )
D (N XiYi Xi Yi )
E N
X
3 i
N
X iYi2
( X i2 Yi2 )
Xi
G (N Yi2 Yi Yi )
H N Yi3 N Xi2Yi ( X i2 Yi2 ) Yi
最小二乘法拟合圆曲线： R2 (x A)2 ( y B)2
R2 x2 2Ax A2 y2 2By B2
令a=-2A,b=-2B, c A2 B2 R2
则：圆的另一形式为：
x2 y2 ax by c 0
1
只需求出参数a,b,c即可以求的圆半径参数：
a A
2
B a 2
Q(a,b, c)
a
2( X i2 Yi2 aX i bYi c) X i 0
①
Q(a,b, c)
b
2( X i2 Yi2 aX i bYi c)Yi 0 ②
Q(a,b, c)
c
2( X i2 Yi2 aX i bYi c) 0 ③
5
由 ① × N- ③ × Xi
9
6
解得： Ca+Db+E=0
Da+Gb+H=0
a
最小二乘法曲线拟合的Matlab程序
最⼩⼆乘法曲线拟合的Matlab程序⽅便⼤家使⽤的最⼩⼆乘法曲线拟合的Matlab程序⾮常⽅便⽤户使⽤,直接按提⽰操作即可;这⾥我演⽰⼀个例⼦:(红⾊部分为⽤户输⼊部分,其余为程序运⾏的结果,结果图为Untitled.fig,Untitled2.fig) 请以向量的形式输⼊x,y.x=[1,2,3,4]y=[3,4,5,6]通过下⾯的交互式图形，你可以事先估计⼀下你要拟合的多项式的阶数，⽅便下⾯的计算.polytool()是交互式函数,在图形上⽅[Degree]框中输⼊阶数,右击左下⾓的[Export]输出图形回车打开polytool交互式界⾯回车继续进⾏拟合输⼊多项式拟合的阶数m = 4Warning: Polynomial is not unique; degree >= number of data points. > In polyfit at 72In zxecf at 64输出多项式的各项系数a = 0.0200000000000001a = -0.2000000000000008a = 0.7000000000000022a = 0.0000000000000000a = 2.4799999999999973输出多项式的有关信息 SR: [4x5 double]df: 0normr: 2.3915e-015Warning: Zero degrees of freedom implies infinite error bounds.> In polyval at 104In polyconf at 92In zxecf at 69观测数据拟合数据x y yh1.0000 3.0000 3.00002.0000 4.0000 4.00003 5 54.0000 6.0000 6.0000剩余平⽅和 Q = 0.000000相关指数 RR = 1.000000请输⼊你所需要拟合的数据点，若没有请按回车键结束程序.输⼊插值点x0 = 3输出插值点拟合函数值 y0 = 5.0000>>结果：untitled.figuntitled2.fig⼀些matlab优化算法代码的分享代码的⽬录如下：欢迎讨论1.约束优化问题：minRosen(Rosen梯度法求解约束多维函数的极值)(算法还有bug) minPF(外点罚函数法解线性等式约束) minGeneralPF(外点罚函数法解⼀般等式约束)minNF(内点罚函数法)minMixFun(混合罚函数法)minJSMixFun(混合罚函数加速法)minFactor(乘⼦法)minconPS(坐标轮换法)(算法还有bug)minconSimpSearch(复合形法)2.⾮线性最⼩⼆乘优化问题minMGN(修正G-N法)3.线性规划：CmpSimpleMthd(完整单纯形法)4.整数规划(含0-1规划)DividePlane(割平⾯法)ZeroOneprog(枚举法)5.⼆次规划QuadLagR(拉格朗⽇法)ActivedeSet(起作⽤集法)6.辅助函数(在⼀些函数中会调⽤)minNT(⽜顿法求多元函数的极值)minMNT(修正的⽜顿法求多元函数极值)minHJ(黄⾦分割法求⼀维函数的极值)7.⾼级优化算法1）粒⼦群优化算法(求解⽆约束优化问题)1>PSO(基本粒⼦群算法)2>YSPSO(待压缩因⼦的粒⼦群算法)3>LinWPSO(线性递减权重粒⼦群优化算法)4>SAPSO(⾃适应权重粒⼦群优化算法)5>RandWSPO(随机权重粒⼦群优化算法)6>LnCPSO(同步变化的学习因⼦)7>AsyLnCPSO(异步变化的学习因⼦)(算法还有bug)8>SecPSO(⽤⼆阶粒⼦群优化算法求解⽆约束优化问题)9>SecVibratPSO(⽤⼆阶振荡粒⼦群优化算法求解五约束优化问题)10>CLSPSO(⽤混沌群粒⼦优化算法求解⽆约束优化问题)11>SelPSO(基于选择的粒⼦群优化算法)12>BreedPSO(基于交叉遗传的粒⼦群优化算法)13>SimuAPSO(基于模拟退⽕的粒⼦群优化算法)2)遗传算法1>myGA(基本遗传算法解决⼀维约束规划问题)2>SBOGA(顺序选择遗传算法求解⼀维⽆约束优化问题)3>NormFitGA(动态线性标定适应值的遗传算法求解⼀维⽆约束优化问题)4>GMGA(⼤变异遗传算法求解⼀维⽆约束优化问题)5>AdapGA(⾃适应遗传算法求解⼀维⽆约束优化问题)6>DblGEGA(双切点遗传算法求解⼀维⽆约束优化问题)7>MMAdapGA(多变异位⾃适应遗传算法求解⼀维⽆约束优化问题)⾃⼰编写的马尔科夫链程序A 代表⼀组数据序列⼀维数组本程序的操作对象也是如此t=length(A); % 计算序列“A”的总状态数B=unique(A); % 序列“A”的独⽴状态数顺序，“E”E=sort(B,'ascend');a=0;b=0;c=0;d=0;Localization=find(A==E(j)); % 序列“A”中找到其独⽴状态“E”的位置for i=1:1:length(Localization)if Localization(i)+1>tbreak; % 范围限定elseif A(Localization(i)+1)== E(1)a=a+1;elseif A(Localization(i)+1)== E(2)b=b+1;elseif A(Localization(i)+1)== E(3)c=c+1;% 依此类推，取决于独⽴状态“E”的个数elsed=d+1;endendT(j,1:tt)=[a,b,c,d]; % “T”为占位矩阵endTT=T;for u=2:1:ttTT(u,:)= T(u,:)- T(u-1,:);endTT; % ⾄此，得到转移频数矩阵Y=sum(TT,2);for uu=1:1:ttTR(uu,:)= TT(uu,:)./Y(uu,1);endTR % 最终得到马尔科夫转移频率/概率矩阵% 观测序列马尔科夫性质的检验：N=numel(TT);uuu=1;Col=sum(TT,2); % 对列求和Row=sum(TT,1); % 对⾏求和Total=sum(Row); % 频数总和for i=1:1:ttxx(uuu,1)=sum((TT(i,j)-(Row(i)*Col(j))./Total).^2./( (Row(i)*Col(j)). /Total)); uuu=uuu+1; % 计算统计量x2endendxx=sum(xx)。
最小二乘法曲线拟合的Matlab程序
最小二乘法曲线拟合的Matlab程序最小二乘法是一种常用的数学优化技术，它通过最小化误差的平方和来找到最佳函数匹配。
在曲线拟合中，最小二乘法被广泛使用来找到最佳拟合曲线。
下面的Matlab程序演示了如何使用最小二乘法进行曲线拟合。
% 输入数据x = [1, 2, 3, 4, 5];y = [2.2, 2.8, 3.6, 4.5, 5.1];% 构建矩阵A = [x(:), ones(size(x))]; % 使用x向量和单位矩阵构建矩阵A% 使用最小二乘法求解theta = (A' * A) \ (A' * y); % 利用最小二乘法的公式求解% 显示拟合曲线plot(x, theta(1) * x + theta(2), '-', 'LineWidth', 2); % 画出拟合曲线hold on; % 保持当前图像plot(x, y, 'o', 'MarkerSize', 10, 'MarkerFaceColor','b'); % 在图像上画出原始数据点xlabel('x'); % 设置x轴标签ylabel('y'); % 设置y轴标签legend('拟合曲线', '原始数据点'); % 设置图例这个程序首先定义了一组输入数据x和y。
然后，它构建了一个矩阵A，这个矩阵由输入数据x和单位矩阵构成。
然后，程序使用最小二乘法的公式来求解最佳拟合曲线的参数。
最后，程序画出拟合曲线和原始数据点。
这个程序使用的是线性最小二乘法，适用于一次曲线拟合。
如果你的数据更适合非线性模型，例如二次曲线或指数曲线，那么你需要使用非线性最小二乘法。
Matlab提供了lsqcurvefit函数，可以用于非线性曲线拟合。
例如：% 非线性模型 y = a * x^2 + b * x + cfun = @(theta, x) theta(1) * x.^2 + theta(2) * x +theta(3);guess = [1, 1, 1]; % 初始猜测值% 使用lsqcurvefit函数求解theta = lsqcurvefit(fun, guess, x, y);% 显示拟合曲线plot(x, fun(theta, x), '-', 'LineWidth', 2); % 画出拟合曲线hold on; % 保持当前图像plot(x, y, 'o', 'MarkerSize', 10, 'MarkerFaceColor','b'); % 在图像上画出原始数据点xlabel('x'); % 设置x轴标签ylabel('y'); % 设置y轴标签legend('拟合曲线', '原始数据点'); % 设置图例这个程序定义了一个非线性函数fun，然后使用lsqcurvefit函数来求解最佳拟合曲线的参数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Q ( a ,,) b c 2 2 2 ( X Y a X b Y ) Y 0 i i i i c i ② b
Q ( a ,,) b c 2 2 2 ( X Ya X b Y )0 ③ i i i i c c
6
由 ① × N- ③ × Xi
最小二乘法圆拟合
1
( x A ) ( y B ) 最小二乘法拟合圆曲线： R
2 2
2 2 2 2 2 R x 2 A x A y 2 B y B
2
令a=-2A,b=-2B,
2 2 2 c A B R
则：圆的另一形式为：
x y a x b y c0
2 i 2 2 i i
2
下面求参数a,b,c使得Q(a,b,c)的值最小即可
5
F(a,b,c)对a,b,c求偏导，令偏导等于0，得到极值点，比较所有极值点的函 数值即可得到最小值。
Q ( a ,, b c ) 2 2 2 ( X Y a X Y ) X 0 i i i b i c i a ①
7
3 2 2 2 E N X N X Y ( X Y ) X i i i i i i
解得： Ca+Db+E=0 Da+Gb+H=0
HD EG a 2 CG D
HC ED b 2 D GC
2 2 ( X Y a Xb Y i i) i i c N
2 2
2
只需求出参数a,b,c即可以求的圆半径参数：
a A 2 a B 2
1 2 2 R a b 4 c 2
3
( X , Y ) , i ( 1 , 2 , 3 N ) 中 到 圆 心 的 距 离 为 d : ii i
2 2 2 d ( X A ) ( YB ) i i
样本集
( X i , Yi )
点
dR ( X AY )( B )R
i i i i
到圆边缘的距离的平方与和半径平方的差为： 2 2 2 2 2
2 i 2 i
X Y + a X b Yc i i
4
令 Q ( a , b , c ) 为 的 平 方 和 ： i
Q ( a , b , c ) [ ( X Y a X b Y c ) ] i i
② × N- ③ × Y i
且令
2 C ( N X XX ) i i i
D ( N X Y XY ) ii i i
2 G ( N Y YY ) i i i
3 2 2 2 H N Y N X Y ( X Y ) Y i ii i i i
8
得A、B、R的估计拟合值：
a A 2
b BБайду номын сангаас 2
1 2 2 R a b 4 c 2
9
t=0:0.01:pi; a=20;%设定圆心X轴数值 b=30;%设定圆心Y轴数值 r=5;%设定圆半径数值 x=a+r*cos(t)+randn(1,315); y=b+r*sin(t)+randn(1,315); plot(x,y); hold on; x=x(:); y=y(:); m=[x y ones(size(x))]\[-(x.^2+y.^2)]; xc = -.5*m(1)%拟合圆心X轴数值 yc = -.5*m(2)%拟合圆心Y轴数值 R = sqrt((m(1)^2+m(2)^2)/4-m(3))%拟合半径数值 plot(xc,yc,'r-x',(xc+R*cos(t)),(yc+R*sin(t)),'r-'); axis equal;