clear
clc
a=[1 -1.5 0.7]'; %对象参数
b=[1 0.5]';
d=3;
na=length(a)-1; %a b阶次
nb=length(b)-1;
L=500; %仿真长度
uk=zeros(d+nb,1);yk=zeros(na,1);%输入输出初值
xik=zeros(na,1); %噪声初值
etak=zeros(d+nb,1);
u=randn(L,1); %输入白噪声序列
xi=sqrt(0.1)*randn(L,1); %方差为0.1的噪声
eta=sqrt(0.25)*randn(L,1);
theta=[a(2:na+1);b]; %真值theta=[a1 ...an b0...bn]
thetae_1=zeros(na+nb+1,1);
thetae_11=zeros(na+nb+1,1);
Rk_1=eye(na+nb+1);
%o=randn(1,L);
%--------随机牛顿算法---------------------------------------------
for k=1:L
phi=[-yk;uk(d:d+nb)];
%e(k)=o(k);
e(k)=a'*[xi(k);xik]-b'*etak(d:d+nb);
y(k)=phi'*theta+e(k); %采集输出数据
%y1(k)=phi'*theta;
R=Rk_1+(phi*phi'-Rk_1)/k; %式25 Hessian矩阵在k时刻的估计值
dR=det(R);
if abs(dR)<10^(-6) %避免矩阵R非奇异
R=eye(na+nb+1);
end
IR=inv(R);
thetae(:,k)=thetae_1+IR*phi*(y(k)-phi'*thetae_1)/k; %更新第k列数据 式21
% thetae1(:,k)=thetae_11+IR*phi*(y1(k)-phi'*thetae_11)/k;
thetae_1=thetae(:,k);
% thetae_11=thetae1(:,k);
Rk_1=R;
for i=d+nb:-1:2
uk(i)=uk(i-1);
etak(i)=etak(i-1);
end
uk(1)=u(k);
etak(1)=eta(k);
for i=na:-1:2
yk(i)=yk(i-1);
xik(i)=xik(i-1);
end
yk(1)=y(k);
xik(1)=xi(k);
end
thetae(:,L)
plot(1:L,thetae)
%hold on
%plot(1:L,thetae1)
title('a b参数估计值')
legend('a_1','a_2','b_0','b_1')
axis([0 L -2 2])