1、系统辨识与建模,第三讲 参数估计的批量法,最小二乘算法 参数估计的一般性质 最小二乘、加权最小二乘估计性质 噪声方差估计 广义最小二乘 偏倚校正算法 辅助变量法 多步最小二乘 相关最小二乘,最小二乘算法,考虑差分方程：A( )y(k)=B( )u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m,），将其写成线性回归模型：y(k)=-a1y(k-1)-any(k-n)+b1u(k-)+bmu(k-m-+1)+ w(k)=a1, , an,b1 , bmT(k)=-y(k-1),-y(k-n),u(k-),u(k-m-+1)Ty(k)=T(k)+w(k) （以下令=1）,汇总的
2、观察误差,y(k)=-a1y(k-1)-any(k-n)+b1u(k-1)+bmu(k-m)+ w(k) y(k+1)=-a1y(k)-any(k-n+1)+b1u(k)+bmu(k-m+1)+w(k+1) y(k+2)=-a1y(k+1)-any(k-n+2)+b1u(k+1)+bmu(k-m+2)+w(k+2)T(k)=- y(k-1)-y(k-n) u(k-1)u(k-m) T(k+1)=- y(k)-y(k-n+1) u(k)u(k-m+1) T(k+2)=- y(k+1)-y(k-n+2) u(k+1)u(k-m+2) T(k+N)=- y(k+N-1)-y(k-n+N) u(k+N
3、-1)u(k-m+N),最小二乘算法,若我们的观测数据可写出N个这样的等式，YN=+WN ，T=(k)，(N+k-1)T YN=T+T WN=(T)-1T YN - (T)-1T WNLS=(T)-1TYN 条件1：E(T WN)=0 （当w(k)为白色时,条件1满足） 条件2: T可逆,当w(k)为白色时,条件1满足,T WN=- y(k-1) -y(k) -y(k+1) -y(k+N-2) w(k)- y(k-2) -y(k-1) -y(k) -y(k+N-1) w(k+1) .- y(k-n) -y(k-n+1) -y(k-n+2)-y(k-n+N-1) .u(k-1) u(k) u(k
4、+1) u(k+N-2) . .u(k-m) u(k-m+1) u(k-m+2) u(k-m+N-1) w(k+N-1) y(k)与w(k)、w(k-1)、w(k-2)、相关。当w(k+1)与w(k)、w(k-1)、w(k-2)、不相关时， y(k)与w(k+1)也不相关。,最小二乘算法,以上结果等同于求使：J=(y(k)-T(k)2 最小的，因此称为最小二乘算法。=0T YN=T-正则方程 T-正则矩阵(对称、正定、可逆),加权最小二乘,若对各次观测数据加不同的权，即求使J=k(y(k)-T(k)2最小的,则得到参数的加权最小二乘估计：LS=(T)-1TYN 的对角元由k构成,参数估计的一般
5、性质,1无偏性 如果E( )=0,( = - ) 或E( )= ，则称估计为无偏的。2有效性 如果无偏估计 满足Cov( )=M-1，则称估计为有效的。其中： 称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。 在一般情况下，有Crammer-Rao不等式： Cov( )=E( - )( - )TM-1,有效性例,对于z=H+w，若噪声w为零均值、协方差阵w=E(WWT)的正态分布噪声，即：WN(0,w)，则输出信号ZN(E(H),w)，即，因此： 于是，M=E(HTw-1H)，其Crammer-Rao不等式为：Cov( )E(HTw-1H)-1。 有效估计也称为最小方差估
6、计、马尔可夫估计。,参数估计的一般性质,3一致性 若估计 为渐近无偏的( E( )=0），且Var( )=0，则称 为一致估计。Var( )=E( - )T( - ),最小二乘、加权最小二乘估计性质,最小二乘、加权最小二乘估计性质,影响最小二乘计算结果的因素:,u(k) T是否可逆 噪声w(k) 大,则 的方差大 白色零均值, 是无偏估计 数据总量N N越大, 的方差越小,最小二乘算法的MATLAB程序,读入数据,读入结构,构造矩阵和Y,计算T和TY,计算Ls,function zta,m,tao=ls(tt) %最小二乘法for MISO, tt的格式为：第一列是系统输出数据，其它列是对应的
7、输入数据 ll=size(tt); %得到数据维数 r=ll(2)-1; m(1)=input(输入 A(z)的阶次); %指定模型结构 tao(1)=0; for i=1:r i %给出输入编号 m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数 tao(i+1)=input(输入B(z)的时延); end n=m(1)+max(tao); %算出一个方程最多使用的数据 lll=ll(1)-n; %算出可列出的方程数 in1=r+1; %构造观测数据矩阵ff kn=0;,for k=1:in1 %每一行中的变量循环 for i=1:lll %列循环 for j=
8、1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延 if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);end if k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号 end;end; kn=kn+m(k); %算出行变量的启始位置 end; for i=1:lll yy(i)=tt(i+n,1); %构造输出向量 end; fa=ff*ff; %最小二乘计算 fay=ff*yy; zz=inv(fa); zta=zz*fay %显示参数和结构 m，tao,Ls.m,Ls1.m,clear all; %最小二
9、乘法for MISO load y3; tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 clear uyr; plot(tt(:,1)ls(tt);,仿真例,1. 无噪声模型：数据文件 y3.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k) 辨识结果（给定结构：m =2 1，tao = 0 2） zta =1.50000.70003.2000,2. 白噪声模型：数据文件 y5.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+w(k
10、) 辨识结果（给定结构：m =2 1，tao = 0 2） zta =1.50270.70323.1935,3. 有色噪声模型：数据文件 y6.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+ 辨识结果（给定结构：m =2 1，tao =0 2） zta =0.97570.17823.3955,4. 白噪声模型b：数据文件 y5b.maty(k)= + w(k)辨识结果（给定结构：m =2 1，tao =0 2） zta =1.16270.37503.3907,噪声方差估计,e=YN- =+WN- =+WN-(T)-1TYN =WN-(T)-1T WN=(I-(T)-
11、1T) WN=B WN， 因为BT=B ,B2=B 所以，E(eTe)=E(WN TB WN)= E(WN T WN- WN T(T)-1T WN)= N- (Tr(T)-1T)= (N-dim)=(y(k)- )2/(N-dim) 在有色噪声环境下，最小二乘估计是有偏的。下面的一些算法对最小二乘进行改进。,追迹:对角线元素之和,广义最小二乘,考虑差分方程：A( )y(k)=B( )u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m）。 如果噪声模型C( )已知，显然用C( )对输入/输出数据进行滤波，则可得到满足最小二乘估计无偏条件的模型：A( ) (k)= B( ) (
12、k)+ w(k)，其中： (k) =C( )y(k)， (k) =C( )u(k)。 在C( )未知时，我们可考虑采用迭代估计的方法去求得。,广义最小二乘的计算步骤,1 令C0( )=1，i=0下标表示迭代次数；0=1000000 2 计算 (k) = Ci( )y(k)， (k) = Ci( )u(k)；i=i+1; 3用最小二乘估计A( ) (k)=B( ) (k)+w(k)中的参数； 4 用估计模型 、 以及各时刻的观测数据，计算出残差 : (k)= ( )y(k)- ( )u(k) 5 计算i= 2 (k)及=i-1-i，如果小于一定数，则结束辨识。否则转下一步。 6 对于噪声模型C(
13、 ) (k)= w(k)，用最小二乘估计出参数，得到更新的Ci( )后，返回2。 以上算法的每一次循环中都要进行滤波和两次求逆。下面的算法将在计算工作量上有所改进。,function zta,m,tao=ls0(tt,m,tao,ll) %最小二乘法for MISO,tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 %m为各多项式中参数个数，应与tt的列数一致；tao为时延；ll=size(tt); n=max(m)+max(tao); %算出一个方程最多使用的数据 lll=ll(1)-n; %算出可列出的方程数 in1=ll(2); %构造观测数据矩阵ff kn=0; for k=
14、1:in1 %每一行中的变量循环 for i=1:lll %列循环 for j=1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延,if k1, ff(i,j+kn)=tt(i+n-jtao+1,k);end if k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号 end;end; kn=kn+m(k); %算出行变量的启始位置 end; for i=1:lll yy(i)=tt(i+n,1); %构造输出向量 end; fa=ff*ff; %最小二乘计算 fay=ff*yy; zz=inv(fa); zta=zz*fay
15、%显示参数和结构 M tao,Ls0.m,clear all; %广义最小二乘法for MISO load y6; tt(:,1)=uyr(1,:); %读入数据赋给tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 clear uy; plot(tt(:,1) ll=size(tt); %得到数据维数 in1=ll(2); r=in1-1; m(1)=input(输入 A(z)的阶次); %指定模型结构 tao(1)=0;,for i=1:r i %给出输入编号 m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数
16、个数 tao(i+1)=input(输入B(z)的时延); End taomax=max(tao+m); mc=input(输入噪声模型C(z)的阶次); c=1;d=1; %初始化滤波器 lb=1;xci=100000;xci0=1000000; while abs(xci0-xci)0.001 %滤波循环 xci0=xci; tt1=filter(c,d,tt); %输入输出滤波. zta,m,tao=ls0(tt1,m,tao,ll); %主最小二乘,Gls.m,for k=1:taomax %计算输出估计y(k)=tt(k,1); %设定初始输出end for k=1+taomax:l
17、l(1) mm=0; for i=1:in1if i=1for j=1:m(i)fb(mm+j)=-tt(k-j,1);end;elsefor j=1:m(i)if k-tao(i)-j+10 fb(mm+j)=tt(k-tao(i)-j+1,i);elsefb(mm+j)=0;endend;end;mm=mm+m(i); end; y(k)=fb*zta; %用LS估计参数求输出估计 end,e=tt(:,1)-y; %计算方程误差 ee=0; %计算损失函数值 for k=1:ll(1)ee=ee+e(k)*e(k); end xci=ee/ll(1) taoc=0; %估计噪声模型 lc
18、=size(e); cc,mc,taoc=ls0(e,mc,taoc,lc); c(2:mc+1)=cc %显示参数和结构 lb=lb+1; end %结束滤波循环 plot(e)，ylabel(error)，xlabel(time),广义最小二乘例,1. 有色噪声模型：数据文件 y6.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+ 辨识结果（给定结构：m =2 1，tao =0 2） zta =1.49500.69433.2075 c = 1.0000 -1.7118 0.8069 迭代次数： 6,2. 白噪声模型b：数据文件 y5b.maty(k)= + w(k
19、)辨识结果（给定结构：m =2 1，tao =0 2） zta =1.50160.69923.1881 c = 1.0000 -1.5067 1.4236 -1.0215 0.5519 -0.1689 迭代次数： 7,偏倚校正算法,仍考虑差分方程：A( )y(k)=B( )u(k)+ w(k)，其中w(k)为白噪声。令(k)= w(k)，则 A( )y(k)= B( )u(k)+(k)。分别写成回归模型： Y=+,=C+W,组合起来有Y=, +W，其最小二乘解为： = ，利用分块矩 阵求逆公式有：=(T)-1TY-(T)-1T=D-1 TMY M=I-(T)-1T D= TM 须要注意， 的求
20、取仍然是一个迭代过程。,=Ay-Bu,偏倚校正算法的计算步骤,1 令C( )=1，用最小二乘法求 =LS=(T)-1TY，并保留、=(T)-1T以及M=I-。 2 计算 =Y- ，并依据 =C +W构造 ，计算 D= TM 。 3 计算 =D-1 TMY，并计算 = 及 =LS- 。 4 若参数已收敛，则结束辨识，否则转2。 以上算法的一次循环中没有滤波，且只有一次求逆。如果将第3步中 的计算改为： =( T )-1 T ，则还可省去D的计算。（这一改进由夏天长首先给出。）本法可能会出现收敛慢的情况，可用对 求均值来解决,Gls2.m,clear all; %偏倚校正法for MISO loa
21、d y6; tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 clear uyr;plot(tt(:,1) ll=size(tt); %得到数据维数 r=ll(2)-1; m(1)=input(输入 A(z)的阶次); %指定模型结构 tao(1)=0; for i=1:r i %给出输入编号 m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数 tao(i+1)=input(输入B(z)的时延); end,mc=input(输入噪声模型C(z)的阶次)
22、; n=m(1)+max(tao); %算出一个方程最多使用的数据 lll=ll(1)-n; %算出可列出的方程数 lb=1;xci=100000;xci0=1000000;c=1; in1=r+1; %构造观测数据矩阵ff kn=0; for k=1:in1 %每一行中的变量循环 for i=1:lll %列循环 for j=1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延 if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);end if k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号 end;end;
23、 kn=kn+m(k); %算出行变量的启始位置 end;for i=1:lll yy(i)=tt(i+n,1); %构造输出向量 end;,fa=ff*ff; %最小二乘计算 fay=ff*yy; zz=inv(fa); zaa=zz*ff; zta=zz*fay %显示参数和结构 M,tao ztab=zta-zta; ztaa=zta; while abs(xci0-xci)0.05 xci0=xci; e(1:mc)=0; %计算输出误差 e(mc+1:lll+mc)=yy-ff*zta; %计算损失函数值 ee=0; for k=1:lll+mcee=ee+e(k)*e(k); en
24、d xci=ee/(lll+mc) %估计噪声模型 for i=1:mc for j=1:lll,fc(j,i)=-e(j+mc-i); end,end cy=e(mc+1:lll+mc); fd=fc*fc; %最小二乘计算 fdcy=fc*cy; zzc=inv(fd); c(2:mc+1)=zzc*fdcy %显示参数和结构 ztab=ztab+zaa*fc*c(2:mc+1); zta=ztaa-ztab/lb %zta=ztaa-zaa*fc*c(2:mc+1) lb=lb+1; end %结束滤波循环,偏倚校正算法例,1. 有色噪声模型：数据文件 y6.mat (1+1.5q-1+
25、0.7q-2)y(k)=q-23.2u(k)+ 辨识结果（给定结构：m =2 1，tao =0 2） zta =1.48790.68313.1976 c =1.0000 -1.7030 0.7983 迭代次数： 20,辅助变量法,分析最小二乘法中，在Y=+W的各项乘上T,然后利用T W的期望值为零得到参数的无偏估计。受此启发，若在Y=+ 的各项乘上T，使其满足以下两个条件：1.T 的期望值为零；2. T可逆，则也可得到参数的无偏估计。下面讨论辅助变量的选取： 设模型为A( )y(k)= B( )u(k)+ (k)，若u(k)与 (k)不相关： a 选取辅助模型D( )z(k)= F( )u(k
26、)，用z(k)、u(k)构造； b 若系统的纯时延已知，则可用u(k-)、u(k)构造； c 用u(k)、u(k)构造 d (k)=D( )w(k),若D( )的阶次n已知,则可用y(k-n)、u(k)构造； e 先求出最小二乘解LS=(T)-1TY ，然后依据 ( )z(k)= ( )u(k)计算出输出估计z(k)，再用z(k)、u(k)构造；,clear all; %辅助变量最小二乘法for MISO load y6; tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 clea
27、r uyr; ll=size(tt); %得到数据维数 r=ll(2)-1; m(1)=input(输入 A(z)的阶次); %指定模型结构 tao(1)=0; for i=1:r i %给出输入编号 m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数 tao(i+1)=input(输入B(z)的时延); end,n=m(1)+max(tao); %算出一个方程最多使用的数据 lll=ll(1)-n; %算出可列出的方程数 in1=r+1; %构造观测数据矩阵ff kn=0; for k=1:in1 %每一行中的变量循环 for i=1:lll %列循环 for
28、j=1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延 if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);end if k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号 end;end; kn=kn+m(k); %算出行变量的启始位置 end; for i=1:lll yy(i)=tt(i+n,1); %构造输出向量 end;,Iv_ls,fa=ff*ff; %最小二乘计算 fay=ff*yy; zz=inv(fa); zta=zz*fay %显示参数和结构%建立辅助变量 %a=1 zta(1:m(1);
29、%b=0 zta(m(1)+1:m(1)+m(2); a=1 1.7 0.72; b=0 0 1; tf=filter(b,a,tt(:,2); kn=0; for k=1:in1 %每一行中的变量循环 for i=1:lll %列循环 for j=1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延 if k1 fh(i,j+kn)=tt(i+n-jtao+1,k);end %if k=1, fh(i,j)=-tt(i+n-jtao+1,2);end %以u(k)为辅助变量,%if k=1, fh(i,j)=-tt(i+n-jtao-tao(2)+1,2)
30、;end %以u(k-tao)为辅助变量 if k=1, fh(i,j)=-tf(i+n-j);end %以辅助模型输出为辅助变量 end;end; kn=kn+m(k); %算出行变量的启始位置 end; for i=1:lll yy(i)=tt(i+n,1); %构造输出向量 end; fa=fh*ff; %辅助变量最小二乘计算 fay=fh*yy; zz=inv(fa); zta=zz*fay %显示参数和结构 M，tao,辅助变量法例,有色噪声模型：数据文件 y6.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+ 辨识结果（给定结构：m =2 1，tao =0
31、 2） 1. 使用辅助模型：a=1 1.7 0.72;b=0 0 1; zta =1.49560.69623.2234,2. 以LS为辅助模型： b=0 0 3.3955; a=1 0.9757 0.1782; zta =1.4908 0.6895 3.2235 3.以u(k)为辅助变量Xzta = 2.7082 -1.8511 3.9692 4.以u(k-tao)为辅助变量（方程病态，溢出）X 5.以y(k-tao)为辅助变量Xzta =4.3079 2.9131 3.8725,2. 白噪声模型b：数据文件 y5b.maty(k)= + w(k)辨识结果（给定结构：m =2 1，tao =0
32、 2） 1. 使用辅助模型：a=1 1.7 0.72;b=0 0 1; zta =1.4740 0.6845 3.36922. 以LS为辅助模型： b=0 0 3.3907; a=1 1.1627 0.3750; zta =1.4456 0.6389 3.3777,多步最小二乘,考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。模型可改写为C()A()y(k)=C()B()u(k)+w(k）或 D()y(k)=F()u(k)+w(k)，其中: D()=C()A()，F()= C()B()-*。 此模型可用最小二乘解出D()、F()。这是第一步。 第二步可有两种不
33、同方法： a 解同次幂方程组 由*式，可得关系：A( )F( )=B( )D( )，两边分别展开，并按 的同次幂相等规则，可列出na+nb+nc个方程：F=H， 其中=a1,.ana,b1,.bnbT，F=f1,.fnb+nc,0,.0,0 0 0. 0 1 0 0. 0 0-f1 0 0. 0 -d1 1 0 0 0 H= -f2 f1 0 . .0 -d2 -d1 1. .0 0. -d1 1. -f2 -f1 .-d2 d1 na+nb+nc 行. . -f3 -f2. . .-fnb+nc. .0 -fnb+nc . -dna+nc.0 0 -fnb+nc. 0 -dna+nc.0 0
34、 0. -fnb+nc 0 .-dna+ncna列 nb列,用最小二乘可求得的无偏估计，即A()、B()。此法也可用来求C()。 b 传递函数等价降阶 由*式，可有： F()/D()=B()/A()，此说明两个传递函数是等价的。对F()/D()施加激励信号u(k)，可得输出z(k)= F()/D()u(k)。用最小二乘法处理 z(k)，u(k)，选择合适的阶次，可得A()、B()的无偏估计。,clear all; %最小二乘法for MISO load y6; tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出
35、数据，其它列是对应的输入数据 clear uyr; ll=size(tt); %得到数据维数 r=ll(2)-1; zta,m,tao=ls10(tt); %辨识高阶模型,%- a=1 zta(1:m(1); %输出仿真 b=zta(m(1)+1:m(1)+m(2); tf=filter(b,a,tt(:,2); tt(:,1)=tf; clear zta tao ff fa fay zz; %- s=第二步，计算低阶模型 ls(tt);,Ms_ls,多步最小二乘例-传递函数等价降阶,有色噪声模型：数据文件 y6.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+ 辨识
36、结果（给定结构：m =2 1，tao =0 2） zta =1.49380.69383.2194,相关最小二乘,设模型为A()y(k)=B()u(k)+E(k)，若u(k)与E(k)不相关，则用u(k-j)乘模型中的各项并求期望，得：A()Ruy(j)= B()Ru(j)，用最小二乘法可得A()、B()的无偏估计。 注意:使用相关最小二乘时,扰动信号序列的周期不要太长,以保证由相关函数组成的模型在求解时不发生病态。另外，计算相关函数时不要以信号序列周期的整数倍来计算。,clear all; %相关最小二乘法for MISO load y6; tt(:,1)=uyr(1,:); %读入数据，并赋
37、给变量tt。 tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据 clear uyr; ll=size(tt); %得到数据维数 r=ll(2)-1; m(1)=input(输入 A(z)的阶次); tao(1)=0; %指定模型结构 for i=1:r i %给出输入编号 m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数,tao(i+1)=input(输入B(z)的时延); end n=m(1)+max(tao); %算出一个方程最多使用的数据 mz=27-1; lll=mz-n; %算出可列出的方程数 in1=r
38、+1; %计算相关函数 for k=1:in1 for j=1:mz rt(j,k)=0; for i=j:j+ll(1)-mz rt(j,k)=rt(j,k)+tt(i,k)*tt(i-j+1,2); %所有变量对U1求相关，长度为27-1，一个M序列周期 end rt(j,k)=rt(j,k)/(ll(1 )-mz); end end,Cov_ls,kn=0; %构造观测数据矩阵ff %最小二乘求解 for k=1:in1 %每一行中的变量循环 for i=1:lll %列循环 for j=1:m(k) %每行变量中的观测数据循环 jtao=j+tao(k); %构造时考虑时延 if k1
39、 ff(i,j+kn)=rt(i+n-jtao+1,k);end if k=1, ff(i,j)=-rt(i+n-jtao,k);end %输出变量变号 end;end; kn=kn+m(k); %算出行变量的启始位置 end; for i=1:lll yy(i)=rt(i+n,1); %构造输出向量 end; fa=ff*ff; %最小二乘计算 fay=ff*yy; zz=inv(fa); zta=zz*fay %显示参数和结构 M，tao,相关最小二乘例,有色噪声模型：数据文件 y6.mat (1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+ 辨识结果（给定结构：m =2
40、1，tao =0 2） zta =1.47070.67113.2319,第四讲 辨识原理,随机逼近法 模型参考自适应辨识方法 极大似然法 预报误差估计法 Bayse估计 1 极大验后参数估计法 2 条件期望参数估计法,随机逼近法,设广义误差e(k)是参数估计值的函数，参数辨识问题可通过极小化e(k)的方差来实现。即求参数使下列准则函数最小：J()=1/2Ee2(k)。J()的负梯度为： =E-e(k) 。如果可求解 =0，则可求得参数的估计。但当e(k)的分布未知时，实际上是不可求解的。 在计算数学中，求二次函数的极小值常采用迭代法。首先给出参数的一个估计值，以二次函数在该参数估计值处的负梯度
41、为修正方向，选取适当的步长后，修正参数估计值，直到收敛。,随机逼近法,仿此，我们有：(k+1)= (k)+(k) 。 如果在求 时不求期望，则得到一个随机的迭代算法，称之为随机逼近法。考虑线性回归模型：y(k)=T(k)+e(k)，其中e(k)是零均值随机噪声。 J()=1/2E y(k)-T(k)2,=E(k) y(k)-T(k)。,例,A 如果假定e(k)是各态遍历的，则梯度为零可用E(k) y(k)-T(k)= (k) y(k)-T(k)=0来代替，由此得到了最小二乘法。 B 应用随机逼近法，可得： (k+1)=(k)+(k)(k)y(k)-T(k)(k)。 (k)的选取应保证迭代收敛，
42、可选取满足如下条件的(k)： (k)0且 ； ，例如(k)=1/k；,例,C 若以准则函数的二阶导数（即海赛矩阵）之逆来参与选择修正方向，则称为牛顿法： (k+1)=(k)+(k)R(k)-1(k)y(k)-T(k)(k) 其中： R(k)= =E(k)T(k) =R(k-1)+(k)(k)T(k)-R(k-1) 牛顿法的优点是收敛速度快。,模型参考自适应辨识方法,考虑线性回归模型：y(k)=T(k)+e(k)，其中e(k)是零均值随机噪声。以输出估计误差为反馈信号，以PI调节器的方式来修正参数： =I+P， I(k)=I(k-1)+P(k)(k)，P(k)=Q(k)(k) 其中P为对称正定矩
43、阵，Q满足P/2+Q0 (k)为广义误差，最简单的取法为 (k)= y(k)-T(k)(k-1)= 0(k),极大似然法,输出z是一个随机变量，它的概率密度p(z|)取决于参数。当获得观测序列ZL=z(1),z(2),.,z(L)T时，由该观测序列组成的联合概率密度p(ZL|)应当取得最大值。（当一个随机事件发生了，我们有理由相信，外部条件一定处于使这随机事件发生的概率最大时的状态。）那么，的极大似然估计就是使p(ZL|)|ML=max的参数估计值。 由于p(ZL|)中ZL已知，因而它只是参数的函数，故称它为的似然函数。有时也记作L(ZL|)。,极大似然法,极大似然原理可用下列等价的表示方式：
44、，求解极大似然估计的下一步是要给出p(ZL|)的具体描述。,独立观测情况,设z(1),z(2),.,z(L)是一组在独立观测条件下获得的随机序列，即各观测值是互相独立的，则p(ZL|)可简化为： p(ZL|)= p(z(1)|)p(z(2)|).p(z(L)|)= p(z(k)|)，其对数似然函数为：L(ZL|)= ln p(z(k)|)。,独立观测情况,设z(k)N(m,2)，即： p(z(k)|)= ，负对数似然函数为： - L(ZL|)= +Lln+(L/2)ln2 当已知时，准则函数就是：J()= 当未知时，可先由min(J()求 及Jmin，再由min(-L(ZL|)求 。，,非独立
45、观测,若z(1),z(2),.,z(L)是非独立观测条件下获得的随机序列，即观测值z(k)是在已有观测z(1),z(2),.,z(k)的基础上得到的，则p(ZL|)应按条件概率的乘法规则写成： p(ZL|)= p(z(L)|ZL-1,)p(ZL-1 |)，以此类推，有：p(ZL|)= p(z(k)|Zk-1,)。,非独立观测,设z(k)N(mk,k)，并取mk为z(k)的条件均值： (k)= Ez(k)| Zk-1，即： p(z(k)|Zk-1,)= ，负对数为： -L(ZL|)= +(L/2)ln2。 当k已知时，就是：J()= 当k未知时，可先由min(J()求 及Jmin，再由min(- L(ZL|)求 。,
