function [ x_post, P_prior, NP_prior ] = stimConst( lambda, tau, LB, A, C, Q, R, y, x0_, P0 )

T = length(y);
[outs,states] = size(C);

[ K,S ] = solveGains( lambda, LB, A, C, Q, R );

x_post = zeros(states,T);                           % stima
P_prior = zeros(states,states,T);                   % varianza a priori
NP_prior = zeros(1,T);                              % traccia varianza a priori 

for t=1:LB+1                                        % da 0 a N nel paper

    P = P0;                                         % inizializzazione ad ogni ciclo
    x = x0_ + gamma_tk(t,1,tau(1))*...              % da xo a priori a x0 a posteriori
        K(:,:,1)*( y(:,1) - C*x0_ );                %
    
    for k=2:t
        g = gamma_tk(t,k-1,tau(k-1));
        K_ = K(:,:,t-k+1);
        P = A*(eye(4) - g*K_*C)*P*(eye(4) - g*K_*C)'*A'+...
            Q+g*A*K_*R*K_'*A';
        x = A*x+g*K_*(y(:,k-1)-C*A*x);
    end
        
    P_prior(:,:,t)=P;                               % assegnazione a fine ciclo
    NP_prior(t)=trace(P);                           %
    x_post(:,t)=x;                                  %
end

for t=LB+2:T
    P1 = P_prior(:,:,t-LB-1);                       % inizializzazione ad ogni ciclo
    x1 = x_post(:,t-LB-1);                          %
    
    for k=0:LB-1
        g = gamma_tk(t,t-LB+k,tau(t-LB+k));
        K_ = K(:,:,LB-k);
        P1 = A*(eye(4) - g*K_*C)*P1*...
            (eye(4) - g*K_*C)'*A'+Q+g*A*K_*R*K_'*A';
        x1 = A*x1+g*K_*(y(:,t-LB+k)-C*A*x1);
    end
    
    P_prior(:,:,t)=P1;                              % assegnazione a fine ciclo
    NP_prior(t)=trace(P1);                          %
    x_post(:,t)=x1;                                 %
end
