%遞推隨機牛頓參數(shù)估計（RSNA）
clear;?
clc;
a=[1?-1.5?0.7]‘;?%對象參數(shù)
b=[1?0.5]‘;?
d=3;?
na=length（a）-1;?%na、nb為A、B階次
nb=length（b）-1;?
L=400;?%仿真長度
uk=zeros（d+nb1）;?%輸入初值：uk（i）表示u（k-i）
yk=zeros（na1）;?%輸出初值
xik=zeros（na1）;?%白噪聲初值
etak=zeros（d+nb1）;?%白噪聲初值
u=randn（L1）;?%輸入采用白噪聲序列
xi=sqrt（0.1）*randn（L1）;?%白噪聲序列
eta=sqrt（0.25）*randn（L1）;?%白噪聲序列
theta=[a（2:na+1）;b];?%對象參數(shù)真值
thetae_1=zeros（na+nb+11）;?%參數(shù)估計初值
Rk_1=eye（na+nb+1）;?
%-------隨機牛頓算法-------------
for?k=1:L
????phi=[-yk;uk（d:d+nb）];
????e（k）=a‘*[xi（k）;xik]-b‘*etak（d:d+nb）;
????y（k）=phi‘*theta+e（k）;?%采集輸出數(shù)據(jù)
???
????R=Rk_1+（phi*phi‘-Rk_1）/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;
???????????
???%更新數(shù)據(jù)
????thetae_1=thetae（: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
disp（‘估計參數(shù)值‘）
thetae（:400）
plot（[1:L]thetae）;
xlabel（‘k‘）;?ylabel（‘參數(shù)估計a、b‘）;
legend（‘a(chǎn)_1‘‘a(chǎn)_2‘‘b_0‘‘b_1‘）;?axis（[0?L?-2?2]）;