function [rbar,Rhat,qbar,Qhat] = sect_7_prob_3_noise_adaptive_filtering(z, phi, H, Q, R, x0, P0)
%   
% Written by:
% -- 
% John L. Weatherwax                2006-08-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

N = length(z)-1; % we need *one* extra measurment in z to compute q_k since it needs \hat{x}_{k+1} 

% Perform Kalman filtering on this system using the given values of Q and R
%
[xhat_plus, P_plus, VV, loglik, innovations, xpred, P_minus, S, Sinv] = kalman_filter(z, phi, H, Q, R, x0, P0);

% calculate the mean of the r_k also known as the innovations: 
%
rbar = zeros(1,N);
for ii=1:N
  start_index = 1; %max(round(0.9*ii),1);
  rbar(ii) = mean( innovations(start_index:ii) );
end

% calculate the estimate of the measurement noise covariance R: 
%
Rhat = zeros(1,N);
ii = 1;
pt1 =   ( innovations(ii) - rbar(ii) )^2;
pt2 = - H * P_minus(ii) * H';
Rhat(1) = pt1 + pt2; 
for ii=2:N
  start_index = 1; %max(round(0.9*ii),1);
  pt1 =   sum( ( innovations(start_index:ii) - rbar(ii) ).^2 )/(ii-1);
  pt2 = - sum( H * P_minus(start_index:ii) * H' )/ii;
  Rhat(ii) = pt1 + pt2; 
end

% calculate the forcing residuals q_k: 
%
q_k = zeros(1,N);
for ii=1:N,
  q_k(ii) = xhat_plus(ii+1) - phi * xhat_plus(ii);
end

% calculate the mean of the forcing residuals: 
%
qbar = zeros(1,N);
for ii=1:N,
  start_index = 1; %max(round(0.9*ii),1);
  qbar(ii) = mean( q_k(start_index:ii) );
end

% calculate the estimate of the process noise covariance Q: 
%
Qhat = zeros(1,N);
ii = 1;
pt1 =   ( q_k(ii) - qbar(ii) )^2;
pt2 = - phi * P_plus(ii) * phi';
Qhat(1) = pt1 + pt2; 
for ii=2:N
  start_index = 1; %max(round(0.9*ii),1);
  pt1 =   sum( ( q_k(start_index:ii) - qbar(ii) ).^2 )/(ii-1);
  pt2 = - sum( phi * P_plus(start_index:ii) * phi' )/ii;
  Qhat(ii) = pt1 + pt2; 
end