Contents

clear all

Generate data from a "true" system with parameters w_0

select some values for the parameters w_0

M = 4;
w_0 = [1.5 -0.8 2.3 0.5]';
% generate the input u(t) as an iid Gaussian sequence
Nd = 200;
Nw = 100;
u = randn(Nd,1);
e = randn(Nd,1);
y  = zeros(Nd,1);
% the output of the sistem w_0 is y(t)= u(t)^T*w_0 + e(t)
for i = (M+1):Nd
    y(i) = w_0' * u(i:-1:(i-M+1)) + e(i);
end
% remove the warming up period 1:Nw
u(1:Nw) = [];
y(1:Nw) = [];
N  = Nd-Nw;

Estimate the parameters w based on the data u and y

Select the data matrix A by using the Covariance method

Ac=[];
for j = M:-1:1
    Ac = [Ac u(j+(0:(N-M)))];
end
yc = y(M:N);
% Estimate the LS parameters
hat_w = (Ac'*Ac)\(Ac'*yc);
[w_0 hat_w]

% use the compact formulations to get the LS solution with i2 = N
i1 = M; i2 = N;
Phi = zeros(M,M); psi = zeros(M,1);
for i = i1:i2
    vect_u = u(i:-1:(i-M+1));
    Phi = Phi + vect_u* vect_u';
    psi = psi + vect_u*y(i);
end
hat_w2 =Phi\psi;
Pn = inv(Phi);
[w_0 hat_w hat_w2]

% use the compact formulations to get the LS solution with i2 = N-1
i1 = M; i2 = N-1;
Phi = zeros(M,M); psi = zeros(M,1);
for i = i1:i2
    vect_u = u(i:-1:(i-M+1));
    Phi = Phi + vect_u* vect_u';
    psi = psi + vect_u*y(i);
end
hat_w3 =Phi\psi;
Pn1 = inv(Phi);
param_true_estN_estN1 = [w_0 hat_w2 hat_w3]

% check the recursive equation for Pn
i= N;
vect_u = u(i:-1:(i-M+1));
T = Pn1*vect_u;
Pnrec = Pn1- (Pn1*vect_u*vect_u'*Pn1')/(1+vect_u'*Pn1*vect_u);
diffP = Pnrec-Pn
% check the recursive equation for w
alpha_n = y(N)-hat_w3'*vect_u;
w_rec = hat_w3 + Pnrec*vect_u*alpha_n;
diff_wn = w_rec-hat_w2
ans =

    1.5000    1.4270
   -0.8000   -0.9892
    2.3000    2.3528
    0.5000    0.5748


ans =

    1.5000    1.4270    1.4270
   -0.8000   -0.9892   -0.9892
    2.3000    2.3528    2.3528
    0.5000    0.5748    0.5748


param_true_estN_estN1 =

    1.5000    1.4270    1.4342
   -0.8000   -0.9892   -0.9931
    2.3000    2.3528    2.3600
    0.5000    0.5748    0.5679


diffP =

   1.0e-17 *

   -0.3469    0.0108    0.0651    0.0054
    0.0108    0.6939    0.0217    0.0867
    0.0651    0.0217    0.3469    0.0054
    0.0054    0.0867    0.0054    0.3469


diff_wn =

   1.0e-15 *

    0.4441
   -0.4441
         0
    0.1110

Use the RLS algorithm for recursively computing w for i = 1...N

Initialize

w = zeros(M,1);
delta = 100; lambda = 1;
P = delta*eye(M);
W = [];
for i = (M+1):N
    vect_u = u(i:-1:(i-M+1));
    piv = vect_u'*P;
    gamma = lambda + piv*vect_u;
    k = piv'/gamma;
    alpha = y(i) - vect_u'*w;
    w = w + k*alpha;
    Pprime = k*piv;
    P = (P-Pprime)/lambda;
    % store the parameters
    W(:,i) = w;
end
param_true_estN_estN1 = [w_0 hat_w2 hat_w3 w]
diffw = hat_w2-w
param_true_estN_estN1 =

    1.5000    1.4270    1.4342    1.4137
   -0.8000   -0.9892   -0.9931   -0.9856
    2.3000    2.3528    2.3600    2.3184
    0.5000    0.5748    0.5679    0.5499


diffw =

    0.0133
   -0.0037
    0.0344
    0.0249

Plot the norm of w-w_0

param_err = [];
for i = 1:N
    w = W(:,i);
    param_err(i) = sqrt( (w-w_0)'*(w-w_0));
end
plot(param_err)
xlabel('Iteration i')
ylabel('|| w-w_0||')