实验二
#include "stdio.h"
float gs(float a[20][20],float b[20],int n )
{
int i,j,k,l;
float s;
k=1;
while(k!=n+1)
{
if(a[k][k]!=0)
{
for(i=k+1;i<=n+1;i++)
{
a[i][k]=a[i][k]/a[k][k];
b[i]=b[i]-a[i][k]*b[k];
for(j=k+1;j<=n+1;j++)
a[i][j]=a[i][j]-a[i][k]*a[k][j];
}
}
k=k+1;
}
for(k=n+1;k>=1;k--)
{
s=0;
for(l=k+1;l<=n+1;l++)
s=s+a[k][l]*b[l];
b[k]=(b[k]-s)/a[k][k];
}
return 0;
}
int main()
{
float a[20][20]={0.0};//定义a矩阵
float c[20][20];//定义c矩阵
float ct[20][20];//定义ct矩阵
float x[20];//定义数组用于存放x的数据
float y[20];//定义数组用于存放y的数据
float b[20]={0.0};//定义b矩阵
int i,j,k,m,n;
printf("输入所求函数的最高次数n:\n");//输入n(求线性的函数输入1。。)scanf("%d",&n);
printf("输入测试数据的组数m:\n");//输入测试数据的组数scanf("%d",&m);
printf("输入x的测试数据%d个:\n",m);//输入x的测试数据m个for(i=1;i<=m;i++)
scanf("%f",&x[i]);
printf("输入y的测试数据%d个:\n",m);//输入y的测试数据m个for(i=1;i<=m;i++)
scanf("%f",&y[i]);
for(i=1;i<=m;i++)//c矩阵第一列赋值为1
c[i][1]=1.0;
//求C[][]
for(j=2;j<=n+1;j++)
for(i=1;i<=m;i++)
c[i][j]=x[i]*c[i][j-1];
//输出C[][]
printf("C矩阵如下:\n");
for(i=1;i<=m;i++)
for(j=1;j<=n+1;j++)
{
printf("%f ",c[i][j]);
if(j==n+1)
printf("\n");
}
//求c的转置矩阵CT[][]
for(i=1;i<=m;i++)
for(j=1;j<=n+1;j++)
ct[j][i]=c[i][j];
//输出CT[][]
printf("CT矩阵如下:\n");
for(i=1;i<=n+1;i++)
for(j=1;j<=m;j++)
{
printf("%f ",ct[i][j]);
if(j==m)
printf("\n");
}
//求a[][]
for(i=1;i<=n+1;i++)
for(j=1;j<=n+1;j++)
for(k=1;k<=m;k++)
a[i][j]+=ct[i][k]*c[k][j];
//输出a[][]
printf("a矩阵如下:\n");
for(i=1;i<=n+1;i++)
for(j=1;j<=n+1;j++)
{
printf("%f ",a[i][j]);
if(j==n+1)
printf("\n");
}
//求b[]
for(i=1;i<=n+1;i++)
for(k=1;k<=m;k++)
b[i]+=ct[i][k]*y[k];
//输出b[]
printf("b矩阵如下:\n");
for(i=1;i<=n+1;i++)
printf("%f ",b[i]);
gs(a,b,n);//调用高斯函数求方程组的解
printf("\n\n");
printf("输出求得的函数的系数为:\n");
for(i=1;i<=n+1;i++)//输出求得的函数的系数printf("a%d=%f ",i-1,b[i]);
printf("\n\n");
return 0;
}
曲线拟合的最小二乘法 学院:光电信息学院 姓名:赵海峰 学号: 200820501001 一、曲线拟合的最小二乘法原理: 由已知的离散数据点选择与实验点误差最小的曲线 S( x) a 0 0 ( x) a 1 1(x) ... a n n ( x) 称为曲线拟合的最小二乘法。 若记 m ( j , k ) i (x i ) j (x i ) k (x i ), 0 m (f , k ) i0 (x i )f (x i ) k (x i ) d k n 上式可改写为 ( k , jo j )a j d k ; (k 0,1,..., n) 这个方程成为法方程,可写成距阵 形式 Ga d 其中 a (a 0,a 1,...,a n )T ,d (d 0,d 1,...,d n )T , 、 数值实例: 下面给定的是乌鲁木齐最近 1个月早晨 7:00左右(新疆时间 )的天气预报所得 到的温度数据表,按照数据找出任意次曲线拟合方程和它的图像。 它的平方误差为: || 2 | 2 ] x ( f
(2008 年 10 月 26~11 月 26) F 面应用Matlab 编程对上述数据进行最小二乘拟合 三、Matlab 程序代码: x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1]; %三次多项式拟合% %九次多项式拟合% %十五次多项式拟合% %三次多项式误差平方和 % %九次次多项式误差平方和 % %十五次多项式误差平方和 % %用*画出x,y 图像% %用红色线画出x,b1图像% %用绿色线画出x,b2图像% %用蓝色o 线画出x,b3图像% 四、数值结果: 不同次数多项式拟和误差平方和为: r1 = 67.6659 r2 = 20.1060 r3 = 3.7952 r1、r2、r3分别表示三次、九次、十五次多项式误差平方和 拟和曲线如下图: a 仁polyfit(x,y,3) a2= polyfit(x,y,9) a3= polyfit(x,y,15) b1= polyval(a1,x) b2= polyval(a2,x) b3= polyval(a3,x) r1= sum((y-b1).A 2) r2= sum((y-b2).A2) r3= sum((y-b3).A2) plot(x,y,'*') hold on plot(x,b1, 'r') hold on plot(x,b2, 'g') hold on plot(x,b3, 'b:o')
曲线拟合(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. 设拟合多项式为: 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(残差)用于生成预测值的误差估计。 [p,s,mu]=polyfit(x,y,n)在拟合过程中,首先对x进行数据标准化处理,以在拟合中消除量纲等影响,mu包含标准化处理过程中使用的x的均值和标准差。 polyval( )为多项式曲线求值函数,调用格式: y=polyval(p,x) [y,DELTA]=polyval(p,x,s) y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。 [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。 如下给定数据的拟合曲线: x=[0.5,1.0,1.5,2.0,2.5,3.0], y=[1.75,2.45,3.81,4.80,7.00,8.60]。 解:MATLAB程序如下: x=[0.5,1.0,1.5,2.0,2.5,3.0]; y=[1.75,2.45,3.81,4.80,7.00,8.60]; p=polyfit(x,y,2) x1=0.5:0.05:3.0; y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b') 运行结果如图1 计算结果为: p =0.5614 0.8287 1.1560 即所得多项式为y=0.5614x^2+0.08287x+1.15560 图1 最小二乘法曲线拟合示例 对比检验拟合的有效性: 例:在[0,π]区间上对正弦函数进行拟合,然后在[0,2π]区间画出图形,比较拟合区间和非拟合区间的图形,考察拟合的有效性。 在MATLAB中输入如下代码: clear x=0:0.1:pi; y=sin(x); [p,mu]=polyfit(x,y,9)
MATLAB机械工程 最小二乘法曲线拟合的应用实例 班级: 姓名: 学号: 指导教师:
一,实验目的 通过Matlab上机编程,掌握利用Matlab软件进行数据拟合分析及数据可视化方法 二,实验内容 1.有一组风机叶片的耐磨实验数据,如下表所示,其中X为使用时间,单位为小时h,Y为磨失质量,单位为克g。要求: 对该数据进行合理的最小二乘法数据拟合得下列数据。 x=[10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 2 0000 21000 22000 23000]; y=[24.0 26.5 29.8 32.4 34.7 37.7 41.1 42.8 44.6 47.3 65.8 87.5 137.8 174. 2] 三,程序如下 X=10000:1000:23000; Y=[24.0,26.5,29.8,32.4,34.7,37.7,41.1,42.8,44.6,47.3,65.8,87.5,137.8,17 4.2] dy=1.5; %拟合数据y的步长for n=1:6 [a,S]=polyfit(x,y,n); A{n}=a;
da=dy*sqrt(diag(inv(S.R′*S.R))); Da{n}=da′; freedom(n)=S.df; [ye,delta]=polyval(a,x,S); YE{n}=ye; D{n}=delta; chi2(n)=sum((y-ye).^2)/dy/dy; end Q=1-chi2cdf(chi2,freedom); %判断拟合良好度 clf,shg subplot(1,2,1),plot(1:6,abs(chi2-freedom),‘b’) xlabel(‘阶次’),title(‘chi2与自由度’) subplot(1,2,2),plot(1:6,Q,‘r’,1:6,ones(1,6)*0.5) xlabel(‘阶次’),title(‘Q与0.5线’) nod=input(‘根据图形选择适当的阶次(请输入数值)’); elf,shg, plot(x,y,‘kx’);xlabel(‘x’),ylabel(‘y’); axis([8000,23000,20.0,174.2]);hold on errorbar(x,YE{nod},D{nod},‘r’);hold off title(‘较适当阶次的拟合’) text(10000,150.0,[‘chi2=’num2str(chi2(nod))‘~’int2str(freedom(nod))])
用matlab实现最小二乘递推算法辨识系统参 数 自动化系统仿真实验室指导教师: 学生姓名班级计082-2 班学号撰写时间: 全文结束》》-3-1 成绩评定: 一.设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作; 3、了解并掌握Matlab中一些函数的作用与使用;二.设计要求最小二乘递推算法辨识系统参数,利用matlab编程实现,设初始参数为零。z(k)-1、5*z(k-1)+0、7*z(k-2)=1*u(k-1)+0、5*u(k-2)+v(k); 选择如下形式的辨识模型:z(k)+a1*z(k- 1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k);三.实验程序 m=3;N=100;uk=rand(1,N);for i=1:Nuk(i)=uk(i)*(-1)^(i-1);endyk=zeros(1,N); for k=3:N yk(k)=1、5*yk(k-1)-0、 7*yk(k-2)+uk(k-1)+0、5*uk(k-2); end%j=100;kn=0;%y=yk(m:j);%psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j- 2)];%pn=inv(psi*psi);%theta=(inv(psi*psi)*psi*y);theta=[0 ;0;0;0];pn=10^6*eye(4);for t=3:Nps=([yk(t-1);yk(t-
2);uk(t-1);uk(t-2)]);pn=pn- pn*ps*ps*pn*(inv(1+ps*pn*ps));theta=theta+pn*ps*(yk(t)-ps*theta);thet=theta;a1=thet(1);a2=thet(2);b1=thet(3);b2= thet(4); a1t(t)=a1;a2t(t)=a2;b1t(t)=b1;b2t(t)=b2;endt=1:N;plot(t,a 1t(t),t,a2t(t),t,b1t(t),t,b2t(t));text(20,1、 47,a1);text(20,-0、67,a2);text(20,0、97,b1);text(20,0、47,b2);四.设计实验结果及分析实验结果图:仿真结果表明,大约递推到第步时,参数辨识的结果基本到稳态状态,即a1=1、5999,b1=1,c1=0、5,d1=-0、7。五、设计感受这周的课程设计告一段落了,时间短暂,意义重大。通过这次次练习的机会,重新把matlab课本看了一遍,另外学习了系统辨识的有关内容,收获颇丰。对matlab的使用更加纯熟,也锻炼了自己在课本中搜索信息和知识的能力。在设计过程中虽然遇到了一些问题,但经过一次又一次的思考,一遍又一遍的检查终于找出了原因所在,也暴露出了前期我在这方面的知识欠缺和经验不足。同时我也进一步认识了matlab软件强大的功能。在以后的学习和工作中必定有很大的用处。
南京信息工程大学实验(实习)报告实验课程数学建模实验名称_ 最小二乘法__ 实验日期 _ 指导老师 专业统计学年级 小组成员 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 实验目的:学会MATLAB软件中曲线拟合方法。 实验内容及要求: 问题1:多项式回归 某种合金中的主要成分为金属A与金属B,经过实验与分析发现,这两种金属成分之和x与膨胀系数y之间有一定的关系。由下面的数据建立描述这种关系的数学表示。 金属成分和x=[37.0 37.5 38.0 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 42.5 43.0]; 膨胀系数y=[3.40 3.00 3.00 2.27 2.10 1.83 1.53 1.70 1.80 1.90 2.35 2.54 2.90]; 注:使用命令:a=polyfit(x,y,n) %求出n阶拟合多项式y=f(x)的系数; y1=polyval(a,x1) %求出f(x)在x1点的函数值,其中x1=37.0:0.5:43.0; plot(x,y,'*r',x1,y1,'-b') %比较原数据和拟合曲线效果; 问题2:非线性回归 设观测到的数据如下: x=20:10:210; y=[0.57 0.72 0.81 0.87 0.91 0.94 0.95 0.97 0.98 0.99 1.00 0.99 0.99 1.00 1.00 0.99 1.00 1.00 0.99 1.00]; 取回归函数为y=b(1)*(1-exp(-b(2)*x)),试估计参数b(1)、b(2)。 注:使用命令: [b,r,j]=nlinfit(x,y,fun,b0); %非线性回归,其中b0为参数初始值,可取 b0=[2,0.1],fun=inline('b(1)*(1-exp(-b(2)*x))','b','x')为内联函数; nlintool(x,y,fun,b0) %绘制非线性回归图。 程序及运行结果: 1、 >> x1=37:0.5:43; y1=[3.40 3.00 3.00 2.27 2.10 1.83 1.53 1.70 1.80 1.90 2.35 2.54 2.90]; >> plot(x1,y1)
最小二乘法的多项式拟 合m a t l a b实现 Document serial number【NL89WT-NY98YT-NC8CB-NNUUT-NUT108】
用最小二乘法进行多项式拟合(matlab 实现) 西安交通大学 徐彬华 算法分析: 对给定数据 (i=0 ,1,2,3,..,m),一共m+1个数据点,取多项式P(x),使 函数P(x)称为拟合函数或最小二乘解,令似的 使得 其中,a0,a1,a2,…,an 为待求未知数,n 为多项式的最高次幂,由此,该问题化为求 的极值问题。由多元函数求极值的必要条件: j=0,1,…,n 得到: j=0,1,…,n 这是一个关于a0,a1,a2,…,an 的线性方程组,用矩阵表示如下:
因此,只要给出数据点 及其个数m ,再给出所要拟合的参数n ,则即可求出未知数矩阵(a0,a1,a2,…,an ) 试验题1 编制以函数 为基的多项式最小二乘拟合程序,并用于对下列数据作三次多项式最小二乘拟合(取权函数wi ≡1) x i y i 总共有7个数据点,令m=6 第一步:画出已知数据的的散点图,确定拟合参数n; x=::;y=[,,,,,,]; plot(x,y,'*') xlabel 'x 轴' ylabel 'y 轴' title '散点图' hold on {} n k k x 0=
因此将拟合参数n设为3. 第二步:计算矩阵 A= 注意到该矩阵为(n+1)*(n+1)矩阵, 多项式的幂跟行、列坐标(i,j)的关系为i+j-2,由此可建立循环来求矩阵的各个元素,程序如下: m=6;n=3; A=zeros(n+1); for j=1:n+1 for i=1:n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)^(j+i-2) end end
、最小二乘拟合原理 x= xl x2 ... xn y= yl y2 ... yn 求m 次拟合 ?力* y 卅…I ZA ; A T A = ZX 茁 X x i - X x i +1 ,- ? ? ? [函Oi …备F =⑷矿丄? A T y 所以m 次拟合曲线为y = a 0 +勿?怎+吐■审+???? +如■牙皿 二、 Matlab 实现程序 function p=funLSM (x, y, m) %x z y 为序列长度相等的数据向量,m 为拟合多项式次数 format short; A=zeros(m+l,m+l); for i=0:m for j=0:m A(i + 1, j + 1)=sum(x.A (i+j)); end b(i+1)=sum(x.A i.*y); end a=A\b 1; p=fliplr (a'); 三、 作业 题1:给出如下数据,使用最小二乘法球一次和二次拟合多项式(取小数点后3位) X 1.36 1.49 1.73 1.81 1.95 2.16 2.28 2.48 Y 14.094 15.069 16.844 17.378 18.435 19.949 20.963 22.495 解:
? x=[1.36 1.49 1.73 1. 81 1. 95 2. 16 2. 28 2. 48]: ? y=[14.094 15.069 16.844 17. 378 18.435 19.949 20.963 22.495]; >> p=funLSM(x, y? 1) P = 7.4639 3.9161 >> p=funLSM(x, y? 2) P = 0.3004 6.3145 4.9763 一次拟合曲线为: y = 7.464x+ 3.91S 二次拟合曲线为: y = +6.315^4-4.976 一次拟合仿真图
第6章曲线拟合的最小二乘法 6.1 拟合曲线 通过观察或测量得到一组离散数据序列,当所得数据比较准确时,可构造插值函数逼近客观存在的函数,构造的原则是要求插值函数通过这些数据点,即。此时,序列与 是相等的。 如果数据序列,含有不可避免的误差(或称“噪音”),如图6.1 所示;如果数据序列无法同时满足某特定函数,如图6.2所示,那么,只能要求所做逼近函数最优地靠近样点,即向量与的误差或距离最小。按与之间误差最小原则作为“最优”标准构造的逼近函数,称为拟合函数。 图6.1 含有“噪声”的数据 图6.2 一条直线公路与多个景点 插值和拟合是构造逼近函数的两种方法。插值的目标是要插值函数尽量靠近离散点;拟合的目标是要离散点尽量靠近拟合函数。 向量与之间的误差或距离有各种不同的定义方法。例如: 用各点误差绝对值的和表示: 用各点误差按模的最大值表示: 用各点误差的平方和表示: 或(6.1)
其中称为均方误差,由于计算均方误差的最小值的方法容易实现而被广泛采用。按 均方误差达到极小构造拟合曲线的方法称为最小二乘法。本章主要讲述用最小二乘法构造拟合曲线的方法。 在运筹学、统计学、逼近论和控制论中,最小二乘法都是很重要的求解方法。例如,它是统计学中估计回归参数的最基本方法。 关于最小二乘法的发明权,在数学史的研究中尚未定论。有材料表明高斯和勒让德分别独立地提出这种方法。勒让德是在1805年第一次公开发表关于最小二乘法的论文,这时高斯指出,他早在1795年之前就使用了这种方法。但数学史研究者只找到了高斯约在1803年之前使用了这种方法的证据。 在实际问题中,怎样由测量的数据设计和确定“最贴近”的拟合曲线?关键在选择适当的拟合曲线类型,有时根据专业知识和工作经验即可确定拟合曲线类型;在对拟合曲线一无所知的情况下,不妨先绘制数据的粗略图形,或许从中观测出拟合曲线的类型;更一般地,对数据进行多种曲线类型的拟合,并计算均方误差,用数学实验的方法找出在最小二乘法意义下的误差最小的拟合函数。 例如,某风景区要在已有的景点之间修一条规格较高的主干路,景点与主干路之间由各具特色的支路联接。设景点的坐标为点列;设主干路为一条直线 ,即拟合函数是一条直线。通过计算均方误差最小值而确定直线方程(见图6.2)。 6.2线性拟合和二次拟合函数 线性拟合 给定一组数据,做拟合直线,均方误差为 (6.2) 是二元函数,的极小值要满足 整理得到拟合曲线满足的方程:
最小二乘递推算法的MATLAB仿真 针对辨识模型,有z(k)-+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k)模型结构,对其进行最小二乘递推算法的MATLAB仿真,对比真值与估计值。更改a1、a2、b1、b2参数,观察结果。 仿真对象:z(k)-1.5*z(k-1)+0.7*z(k-2)=u(k-1)+0.5*u(k-2)+v(k) 程序如下: L=15; y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的初始值 for i=1:L; %移位循环 x1=xor(y3,y4); x2=y1; x3=y2; x4=y3; y(i)=y4; %取出作为输出信号,即M序列 if y(i)>0.5,u(i)=-0.03; %输入信号 else u(i)=0.03; end y1=x1;y2=x2;y3=x3;y4=x4; end figure(1); stem(u),grid on z(2)=0;z(1)=0; for k=3:15; z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %输出采样信号 end c0=[0.001 0.001 0.001 0.001]'; %直接给出被识别参数的初始值 p0=10^6*eye(4,4); %直接给出初始状态P0 E=0.000000005; c=[c0,zeros(4,14)]; e=zeros(4,15); for k=3:15; %开始求k h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %开始求k的值 d1=z(k)-h1'*c0;c1=c0+k1*d1; e1=c1-c0; e2=e1./c0; %求参数的相对变化 e(:,k)=e2; c0=c1; c(:,k)=c1; p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出P(k)的值 p0=p1;
普通最小二乘法(OLS ) 普通最小二乘法(Ordinary Least Square ,简称OLS ),是应用最多的参数估计方 法,也是从最小二乘原理出发的其他估计方法的基础,是必须熟练掌握的一种方法。 在已经获得样本观测值 i i x y ,(i=1,2,…,n )的情况下 (见图2.2.1中的散点),假如模型(2.2.1)的参数估计量 已经求得到,为^0β和^ 1β,并且是最合理的参数估计量,那 么直线方程(见图2.2.1中的直线) i i x y ^ 1^0^ββ+= i=1,2,…,n (2.2.2) 应该能够最好地拟合样本数据。其中 ^ i y 为被解释变量的估计值,它是由参数估计量和解释 变量的观测值计算得到的。那么,被解释变量的估计值与观测值应该在总体上最为接近,判断的标准是二者之差的平方和最小。 ),()(102 2101ββββQ u x y Q i i n i i ==--=∑∑= ()() ),(min ????1 02 1 102 12?,?1 1 ββββββββQ x y y y u Q n i i n i i i =--=-==∑∑∑== (2.2.3) 为什么用平方和?因为二者之差可正可负,简单求和可能将很大的误差抵消掉,只有平方和才能反映二者在总体上的接近程度。这就是最小二乘原则。那么,就可以从最小二乘原则和样本观测值出发,求得参数估计量。 由于 2 1 ^ 1^01 2 ^ ))(()(∑∑+--=n i i n i i x y y y Q ββ= 是 ^ 0β、^ 1β的二次函数并且非负,所以其极小值总是存在的。根据罗彼塔法则,当Q 对^ 0β、 ^ 1β的一阶偏导数为0时,Q 达到最小。即
1 曲线拟合的线性最小二乘法及其MATLAB 程序 例7.2.1 给出一组数据点),(i i y x 列入表7–2中,试用线性最小二乘法求拟合曲线,并用(7.2),(7.3)和(7.4)式估计其误差,作出拟合曲线. 表7–2 例7.2.1的一组数据),(y x 解 (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.12 6.50 68.04]; plot(x,y,'r*'), legend('实验数据(xi,yi)') xlabel('x'), ylabel('y'), title('例7.2.1的数据点(xi,yi)的散点图') 运行后屏幕显示数据的散点图(略). (3)编写下列MA TLAB 程序计算)(x f 在),(i i y x 处的函数值,即输入程序 >> syms a1 a2 a3 a4 x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4 运行后屏幕显示关于a 1,a 2, a 3和a 4的线性方程组 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.^2; J=sum(fy.^2) 运行后屏幕显示误差平方和如下 J= (-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+2 89/100*a2-17/10*a3+a4+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a 2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/ 2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2 为求4321,,,a a a a 使J 达到最小,只需利用极值的必要条件0=??k a J )4,3,2,1(=k ,
生活中的计算方法应用实例——— 最小二乘法,用MATLAB实现1. 数值实例 下面给定的是某市最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度 天数 1 2 3 4 5 6 7 8 9 10 温度9 10 11 12 13 14 13 12 11 9 天数11 12 13 14 15 16 17 18 19 20 温度10 11 12 13 14 12 11 10 9 8 天数21 22 23 24 25 26 27 28 29 30 温度7 8 9 11 9 7 6 5 3 1 下面用MATLAB编程对上述数据进行最小二乘拟合,按照数据找出任意次曲线拟合方程和它的图像。 2、程序代码 x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7, 6,5,3,1]; a1=polyfit(x,y,3) %三次多项式拟合% a2= polyfit(x,y,9) %九次多项式拟合% a3= polyfit(x,y,15) %十五次多项式拟合% b1= polyval(a1,x) b2= polyval(a2,x) b3= polyval(a3,x) r1= sum((y-b1).^2) %三次多项式误差平方和% r2= sum((y-b2).^2) %九次次多项式误差平方和% r3= sum((y-b3).^2) %十五次多项式误差平方和% plot(x,y,'*') %用*画出x,y图像% hold on plot(x,b1, 'r') %用红色线画出x,b1图像% hold on plot(x,b2, 'g') %用绿色线画出x,b2图像% hold on plot(x,b3, 'b:o') %用蓝色o线画出x,b3图像% 3、数值结果 不同次数多项式拟合误差平方和为: r1=67.6659
4. 设某物理量Y与X 满足关系式Y=aX2+bX+c,实验获得一批数据如 下表,试辨识模型参数a,b和c 。(50分) X 1.01 2.03 3.02 4.015 6.027.038.049.0310 Y9.6 4.1 1.30.40.050.10.7 1.8 3.89.0单,最后给出结果及分析。 (1) 问题描述: 由题意知,这是一个已知模型为Y=aX2+bX+c,给出了10组实验输入输出数据,要求对模型参数a,b,c进行辨识。这里对该模型参数辨识采用递推最小二乘法。 (2) 参数估计原理 对该模型参数辨识采用递推最小二乘法,即RLS(recurisive least square),它是一种能够对模型参数进行在线实时估计的辨识方法。 其基本思想可以概括为:新的估计值=旧的估计值+修正项 下面将批处理最小二乘法改写为递推形式即递推最小二乘参数估计的计算方法。 批处理最小二乘估计为,设k时刻的批处理最小二乘估计为: 令 K时刻的最小二乘估计可以表示为 = =;式中,因为要推导出P(k)和K(k)的递推方程,因此这里介绍一下矩阵求逆引理:设A、(A+BC)和(I+)均为非奇异方阵,则通过运用矩阵求逆引理,把复杂的矩阵求逆转化为标量求倒数,大大减小了计算量。与间的递推关系。最终得到递推最小二乘参数递推估计公式如下: (3)程序流程图 (如右图1所示) 递推最小二乘法(RLS)步骤如下: 已知:、和d。 Step 1 :设置初值和P(0),输入初始数据; Step2 :采样当前输出y(k)、和输入u(k) Step3 :利用上面式计算、和; Step4 :kk+1,返回step2,继续循环。
方便大家使用的最小二乘法曲线拟合的Matlab程序 非常方便用户使用,直接按提示操作即可;这里我演示一个例子:(红色部分为用户输入部分,其余为程序运行的结果,结果图为Untitled.fig,Untitled2.fig) 请以向量的形式输入x,y. x=[1,2,3,4] y=[3,4,5,6] 通过下面的交互式图形,你可以事先估计一下你要拟合的多项式的阶数,方便下面的计算. polytool()是交互式函数,在图形上方[Degree]框中输入阶数,右击左下角的[Export]输出图形 回车打开polytool交互式界面 回车继续进行拟合 输入多项式拟合的阶数m = 4 Warning: Polynomial is not unique; degree >= number of data points. > In polyfit at 72 In zxecf at 64 输出多项式的各项系数 a = 0.0200000000000001 a = -0.2000000000000008 a = 0.7000000000000022 a = 0.0000000000000000 a = 2.4799999999999973 输出多项式的有关信息 S R: [4x5 double] df: 0 normr: 2.3915e-015 Warning: Zero degrees of freedom implies infinite error bounds. > In polyval at 104 In polyconf at 92 In zxecf at 69 观测数据拟合数据 x y yh 1.0000 3.0000 3.0000 2.0000 4.0000 4.0000 3 5 5 4.0000 6.0000 6.0000 剩余平方和 Q = 0.000000 标准误差 Sigma = 0.000000 相关指数 RR = 1.000000 请输入你所需要拟合的数据点,若没有请按回车键结束程序. 输入插值点x0 = 3 输出插值点拟合函数值 y0 = 5.0000
最小二乘一次完成算法的MATLAB仿真 针对辨识模型,有z(k)-+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k)模型结构,对其进行最小二乘一次完成算法的MATLAB仿真,对比真值与估计值。更改a1、a2、b1、b2参数,观察结果。 仿真对象:z(k)-1.5*z(k-1)+0.7z(k-2)=u(k-1)+0.5*u(k-2)+v(k) 程序如下: u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1];%输入信号为一个周期的M序列 z=zeros(1,16); for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%以理想输出值作为观测值 end subplot(3,1,1) stem(u) subplot(3,1,2) i=1:1:16; plot(i,z) subplot(3,1,3) stem(z),grid on u,z %显示输入信号与输出观测信号 L=14; HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8); -z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵HL赋值 ZL=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)] %给样本矩阵ZL 赋值 c1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c=c2*c3 a1=c(1),a2=c(2),b1=c(3),b2=c(4) 程序运行结果如下: u = -1 1 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 1 z = Columns 1 through 9 0 0 0.5000 0.2500 0.5250 2.1125 4.3012 6.4731 6.1988 Columns 10 through 16 3.2670 -0.9386 -3.1949 - 4.6352 -6.2165 - 5.5800 -2.5185 HL = 0 0 1.0000 -1.0000 -0.5000 0 -1.0000 1.0000 -0.2500 -0.5000 1.0000 -1.0000
用最小二乘法进行多项式拟合(matlab实现) 西安交通大学 徐彬华 算法分析: 对给定数据(i=0 ,1,2,3,..,m),一共m+1个数据点,取多项式P(x),使 函数P(x)称为拟合函数或最小二乘解,令似的使得 其中,a0,a1,a2,…,an为待求未知数,n为多项式的最高次幂,由此,该问题化为求 的极值问题。由多元函数求极值的必要条件: j=0,1,…,n 得到: j=0,1,…,n
因此,只要给出数据点 及其个数m ,再给出所要拟合的参数n ,则即可求出未知数矩阵(a0,a1,a2,…,an ) 试验题1 编制以函数 为基的多项式最小二乘拟合程序,并用于对 下列数据作三次多项式最小二乘拟合(取权函数wi ≡1) x i -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 y i -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552 总共有7个数据点,令m=6 第一步:画出已知数据的的散点图,确定拟合参数n; x=-1.0:0.5:2.0;y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552]; plot(x,y,'*') xlabel 'x 轴' ylabel 'y 轴' title '散点图' hold on {}n k k x 0=
因此将拟合参数n设为3. 第二步:计算矩阵 A= 注意到该矩阵为(n+1)*(n+1)矩阵, 多项式的幂跟行、列坐标(i,j)的关系为i+j-2,由此可建立循环来求矩阵的各个元素,程序如下: m=6;n=3; A=zeros(n+1); for j=1:n+1
4.最小二乘法线性拟合 我们知道,用作图法求出直线的斜率a 和截据b ,可以确定这条直线所对应的经验公式,但用作图法拟合直线时,由于作图连线有较大的随意性,尤其在测量数据比较分散时,对同一组测量数据,不同的人去处理,所得结果有差异,因此是一种粗略的数据处理方法,求出的a 和b 误差较大。用最小二乘法拟合直线处理数据时,任何人去处理同一组数据,只要处理过程没有错误,得到的斜率a 和截据b 是唯一的。 最小二乘法就是将一组符合Y=a+bX 关系的测量数据,用计算的方法求出最佳的a 和b 。显然,关键是如何求出最佳的a 和b 。 (1) 求回归直线 设直线方程的表达式为: bx a y += (2-6-1) 要根据测量数据求出最佳的a 和b 。对满足线性关系的一组等精度测量数据(x i ,y i ),假定自变量x i 的误差可以忽略,则在同一x i 下,测量点y i 和直线上的点a+bx i 的偏差d i 如下: 111bx a y d --= 222bx a y d --= n n n bx a y d --= 显然最好测量点都在直线上(即d 1=d 2=……=d n =0),求出的a 和b 是最理想的,但测量点不可能都在直线上,这样只有考虑d 1、d 2、……、d n 为最小,也就是考虑d 1+d 2+……+d n 为最小,但因d 1、d 2、……、d n 有正有负,加起来可能相互抵消,因此不可取;而|d 1|+ |d 2|+……+ |d n |又不好解方程,因而不可行。现在采取一种等效方法:当d 12+d 22+……+d n 2 对a 和b 为最小时,d 1、d 2、……、d n 也为最小。取(d 12+d 22+……+d n 2 )为最小值,求a 和b 的方法叫最小二乘法。 令 ∑== n i i d D 1 2=21 1 2][i i n i n i i b a y d D --== ∑∑== (2-6-2) D 对a 和b 分别求一阶偏导数为: ][211∑∑==---=??n i i n i i x b na y a D ][21 2 11∑∑∑===---=??n i i n i i n i i i x b x a y x b D
最小二乘法拟合探究 吴春晖 (中国海洋大学海洋环境学院山东青岛 266100) 摘要: 本文的拟合对象为含多个变量的待定系数的多项式。通过最小二乘法对多项式作出拟合,以向量矩阵的形式来解出待定的系数。在matlab中,通过算法,写出具体的解法。之后,先对最小二乘法的准确性作出检验,分析该方法在应对复杂情况的误差。在检验该方法的可行性之后,对给定的变量值进行拟合与解题。同时,本文将对基于Laguerre多项式的最小二乘法进行分析检验, 关键词:最小二乘法拟合多变量Laguerre多项式 引言: 在之前的计算方法中,在给出已知节点后,如果需要根据给出的节点来确定未知节点的值,我们需要运用插值。在对插值的精准性进行分析后,我们发现不同插值方式的误差都极大,而且插值所得出的函数的特征由插值方式所决定,并不能反映具体的节点原来可能的规律与分布。所以,拟合的方法相比插值而言,并不要求函数值在原节点处的值相等,却能在一定程度上反映原函数的规律。在该文中,我们主要运用最小二乘法进行拟合。
目录 第一章matlab最小二乘法拟合程序 (3) 1.1最小二乘法拟合的数学法 (3) 1.2 编写最小二乘法的matlab拟合程序 (3) 1.2.1程序算法 (3) 1.2.2 最小二乘法拟合的程序 (4) 1.3程序的分析说明 (4) 第二章最小二乘拟合法的检验及应用 (5) 2.1 最小二乘法拟合的检验 (5) 2.2最小二乘法拟合的实际应用 (7) 第三章Laguerre多项式的最小二乘拟合 (8) 3.1 算法与程序 (8) 3.2检验与分析 (9) 第四章最小二乘法拟合的分析总结 (11)