文档库 最新最全的文档下载
当前位置:文档库 › matlab最小二乘法

matlab最小二乘法

matlab最小二乘法
matlab最小二乘法

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,继续循环。

图1 程序流程图

(4) Matlab仿真程序、输出参数估计值、

参数估计变化轨迹图像、结果分析

仿真程序如下:

X=[1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10]; Y=[9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.8 3.8 9.0];

%实验输入数据、实验输出数据

syms a b c % 定义待辨识参数

theta=[a;b;c]; %theta包含待辨识参数a,b,c

theta1=zeros(3,1); %对象参数初始化

P=10^6*eye(3) %构造初始P阵

for k=1:10 %仿真步长范围1到10

phi=[X(k)*X(k);X(k);1];

%y=aX*X+bX+c=phi'*theta

%theta=[a;b;c];phi=

[X(k)*X(k);X(k);1]

K=P*phi/(1+phi'*P*phi); %递推最小二乘法K阵的递推公式

theta=theta1+K*(Y(k)-phi'*theta1); %theta的递推公式

P=(eye(3)-K*phi')*P; %递推最小二乘法P阵的递推公式

theta1=theta; %theta的最终估计向量

theta2(:,k)=theta; %theta估计向量矩阵化,目的是为了

%下面的plot仿真图像输出

end

theta1 %输出参数估计值

plot([1:10],theta2) %输出参数逐步递推估计的轨迹图像

xlabel('k'); %设置横坐标为步长k

ylabel('参数估计a,b,c'); %纵坐标为估计参数a,b,c legend('a','b','c'); %标示相应曲线对应的参数

axis([1 10 -10 20]); %设置坐标轴范围

P =

1000000 0 0

0 1000000 0

0 0 1000000

输出参数估计值、参数估计变化轨迹图像:

theta1 =

0.4575

-5.0734

13.3711

图 2 参数估计逐步变化轨迹图像

结果分析:

通过matlab仿真可知,由递推最小二乘法辨识到的参数为:

a=0.4575;b=-5.0734;c=13.3711

所以Y=0.4575-5.0734X+13.3711 。通过输入输出的数据验算,模型参数a,b,c的递推最小二乘估计在误差允许的范围内基本符合要求。从上面参数估计逐步变化轨迹图像看出,参数a,b,c的递推估计曲线变化比较明显,曲线也不够平滑,这可能是由于题目只给出了10组实验输入输出数据,数据长度不算大,导致了参数估计的误差增大和估计精度不够高。

曲线拟合的最小二乘法matlab举例

曲线拟合的最小二乘法 学院:光电信息学院 姓名:赵海峰 学号: 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')

最小二乘法曲线拟合 原理及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. 设拟合多项式为: 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)

php九九乘法表

"; for($j = 1;$j <= 9;$j++){ echo ""; for($z = 0;$z <9 - $j;$z++){ echo " "; } for($i = $j;$i >= 1;$i--){ echo "$i*$j=".$i*$j.""; } echo ""; } echo "";*/ //第三种 /*echo "

"; for($j = 9;$j >= 1;$j--){ echo ""; for($z = 1;$z <= 9 - $j;$z++ ){ echo ""; //echo "$i*$j=".$i*$j."";

} for($i = 1;$i <= $j;$i++){ echo "

"; } echo ""; } echo "
$i*$j=".$i*$j."
";*/ /*第二种 for($j = 9;$j >= 1;$j--){ for($i = 1; $i <= $j; $i++){ //echo "$i x $j = ".$i * $j; //echo " "; echo "$i x $j =".$i*$j."  "; } echo "
"; } */ /*第一种 for($j = 1;$j <= 9;$j++){

Matlab最小二乘法曲线拟合的应用实例

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实现最小二乘递推算法辨识系统参数

用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循环语句

matlab 基本语句 1.循环语句for for i=s1:s3:s2 循环语句组 end 解释:首先给i赋值s1;然后,判断i是否介于s1与s2之间;如果是,则执行循环语句组,i=i+s3(否则,退出循环.);执行完毕后,继续下一次循环。 例:求1到100的和,可以编程如下: sum=0 for i=1:1:100 sum=sum+i end 这个程序也可以用while语句编程。 注:for循环可以通过break语句结束整个for循环. 2.循环语句while 例:sum=0;i=1; while(i<=100) sum=sum+i;i=i+1; end 3.if语句 if(条件) 语句 end if(条件) 语句 else 语句 end if(条件) 语句 elseif 语句 end 4.关系表达式:

=,>,<,>=,<=,==(精确等于) 5.逻辑表达式:|(或),&(且) 6.[n,m]=size(A)(A为矩阵) 这样可以得到矩阵A的行和列数 n=length(A),可以得到向量A的分量个数;如果是矩阵,则得到矩阵A的行与列数这两个数字中的最大值。 7.!后面接Dos命令可以调用运行一个dos程序。 8.常见函数: poly():为求矩阵的特征多项式的函数,得到的为特征多项式的各个系数。如 a=[1,0,0;0,2,0;0,0,3],则poly(a)=1 -6 11 -6。相当于poly(a)=1入^3+(-6)入^2+11入+(-6)。 compan():可以求矩阵的伴随矩阵. sin()等三角函数。 MATLAB在数学建模中的应用(3) 一、程序设计概述 MATLAB所提供的程序设计语言是一种被称为第四代编程语言的高级程序设计语言,其程序简洁,可读性很强,容易调试。同时,MATLAB的编程效率比C/C++语言要高得多。 MATLAB编程环境有很多。常用的有: 1.命令窗口 2.word窗口 3.M-文件编辑器,这是最好的编程环境。 M-文件的扩展名为“.m”。M-文件的格式分为两种: ①λ M-脚本文件,也可称为“命令文件”。 ② M-函数文件。这是matlab程序设计的主流。λ 保存后的文件可以随时调用。 二、MATLAB程序结构 按照现代程序设计的观点,任何算法功能都可以通过三种基本程序结构来实现,这三种结构是:顺序结构、选择结构和循环结构。其中顺序结构是最基本的结构,它依照语句的自然顺序逐条地执行程序的各条语句。如果要根据输入数据的实际情况进行逻辑判断,对不同的结果进行不同的处理,可以使用选择结构。如果需要反复执行某些程序段落,可以使用循环结构。 1 顺序结构 顺序结构是由两个程序模块串接构成。一个程序模块是完成一项独立功能的逻辑单元,它可以是一段程序、一个函数,或者是一条语句。 看图可知,在顺序结构中,这两个程序模块是顺序执行的,即先执行<程序

最小二乘法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)

最小二乘法的多项式拟合matlab实现

最小二乘法的多项式拟 合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

九九乘法表的C语言代码

九九乘法表的C语言代码,黄路平编写与2012.3.6 代码一:#include int main() { int i=1,j; for (i=1,j=1;j<=9;j++) { if( j==1) printf("%d*%d=%d\n",i,j,i*j); else {for (i=1;i<=j;i++) printf("%d*%d=%d\t",i,j,i*j); printf("\n"); } } } 代码二:switch语句 #include int main() { int i=1,j; for (i=1,j=1;j<=9;j++) { switch(j) { case 1:printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 2: for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 3:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 4:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 5:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break;

case 6:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 7:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 8:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 9:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; } } }

最小二乘法Matlab自编函数实现及示例.docx

、最小二乘拟合原理 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 一次拟合仿真图

最小二乘法MATLAB程序及结果

最小二乘递推算法的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;

九九乘法表源代码(vb)

Private Sub Command1_Click() For i = 1 To 4 For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command2_Click() For i = 1 To 4 Print Tab(i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command3_Click() For i = 1 To 4 Print Tab(5 - i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command4_Click() For i = 1 To 9 For j = 1 To 9 Print i; "*"; j; "="; i * j; Next j Print Next i End Sub Private Sub Command5_Click() Dim se As String Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j

Print Tab((j - 1) * 9); se; Next j Picture1.Print Next i End Sub Private Sub Command6_Click() End End Sub Private Sub Command7_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = i To 9 se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command8_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To i se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command9_Click() Cls End Sub Private Sub Picture1_Click() Dim se As String Picture1.Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j Picture1.Print Tab((j - 1) * 9); se; Next j

VB程序设计-九九乘法表

VB程序设计-九九乘法表 (作者:草原飞狼 2014年5月31日) 声明:仅供学习与交流使用,高手请飘过,谢谢!所有代码都是个人亲自编写并调试成功。 布局

运行界面(1) 运行界面(2) 源代码如下: Private Sub Command1_Click() Rem 乘法表算法一 Dim i, j As Integer For i = 1 To 9 For j = 1 To i expss = i & "*" & j & "=" & i * j Form1.Picture1.Print Tab((j - 1) * 8 + 1); '控制每个输出表达式的输出位置 Form1.Picture1.Print expss; Next j Form1.Picture1.Print Next i End Sub Private Sub Command2_Click() Rem 乘法表算法二 Dim i, j As Integer For i = 9 To 1 Step -1 For j = 1 To i expss = j & "*" & i & "=" & i * j

Form1.Picture1.Print Tab((j - 1) * 10 + 1); '控制每个输出表达式的输出位置 Form1.Picture1.Print expss; Next j Form1.Picture1.Print Next i End Sub Private Sub Command3_Click() Rem 清空 Form1.Picture1.Cls End Sub Private Sub Command4_Click() Rem 退出 Dim int_msg As Integer int_msg = MsgBox("单击“是”退出程序,单击“否”返回程序!", vbYesNo + vbQuestion + vbDefaultButton1, "退出提示") If int_msg = vbYes Then Unload Me End If End Sub Private Sub Form_Load() Rem 初始化 With Picture1 .FontName = "宋体" .FontSize = 12 .FontBold = True .ForeColor = &H80FF& End With End Sub

曲线拟合_线性最小二乘法及其MATLAB程序

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最小二乘法

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程序

方便大家使用的最小二乘法曲线拟合的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程序及结果

最小二乘一次完成算法的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

相关文档
相关文档 最新文档