function [R, mu] = recursive_covariance(x, R, mu, k)
% function recursive_covariance
% 
% A script that employs an adaptive update of the mean and covariance
% matrix,.
%
% %First set up the image
% figh = figure(1); clf; 
% I = zeros(256, 256);
% B1 = roipoly(I);
% clf(figh);
% 
% % get the xy values of the image patch
% [yind, xind] = find(B1);
% 
% % this is the actual covariance matrix.
% real_R = cov([xind, yind]);
% real_mu = mean([xind, yind])';
% 
% sse_r = zeros(length(xind), 1);
% sse_mu = zeros(length(xind), 1);
% % Initialize the covariance and mean
% R = zeros(2,2);
% mu = zeros(2,1);
% 
% % this k denotes the sample number, so we could imagine having another for loop where
% % k would not start at 1.
% k =1;
% for i=1:length(xind)
%     x = [xind(i); yind(i)];
%    [R, mu] = recursive_covariance(x, R, mu, k);
%     k = k+1;
%     sse_r(i) = sum(sum((R - real_R).^2));
%     sse_mu(i) = sum((mu - real_mu).^2);
% end
% figure(figh); 
% subplot(2,1,1);
% plot(sse_r);
% grid on; 
% axis tight;
% title('Covariance error');
% 
% subplot(2,1,2);
% plot(sse_mu);
% grid on; 
% axis tight;
% title('Mean error');
% axis tight;
%
% Dr. A. I. Hanna (2007)


mu = mu_step(mu, x, k);
R = R_step(R, x, mu, k);
return;
%%%
%
%%%
function R = R_step(R, x, mu, k)
x = x - mu;
R = ((k-1)/k)*R + (x*x')/k;
%%%
%
%%%
function mu = mu_step(mu, x, k)
mu = ((k-1)/k)*mu + x/k;
