function [phi, sigma_v, sigma_w, psiVec, phiNewVec,sigmavNewVec,...
    sigmawNewVec] = onlineEM(phi, psiVec, xBuffer, deltaYBuffer,...
    linkLocs, linkDist, counter, phiVec, sigmavVec, sigmawVec, sigma)
% PURPOSE : Online EM to calculate the sufficient statistics and update
%           the unknown model parameters.

% INPUTS  : - phi = the attenuation paramter in pixel-free model.
%           - psiVec = unknown parameter set.
%           - xBuffer = particle locations collected in previous particle
%             filtering runs for this online EM run.
%           - deltaYBuffer = RSS collected in previous time steps for this 
%                            online EM run.
%           - linkLocs = locations of node pair for each link.
%           - linkDist = Distance of a link.
%           - counter = count the number of particle filtering runs.
%           - phiVec = vector to store \phi in each online EM runs.
%           - sigmavVec = vector to store transition noise in each 
%                         online EM runs.
%           - sigmawVec = vector to store measurement noise in each 
%                         online EM runs.
%           - sigma = variable controls rate of decay in pixel-free model

% OUTPUTS : - phi = updated attenuation paramter.
%           - sigma_v = updated transition noise.
%           - sigma_w = updated measurement noise.
%           - psiVec = updated vector to store unknown parameter set.
%           - phiNewVec = updated vector to store \phi in each online EM
%                         runs.
%           - sigmavNewVec = updated vector to store transition noise in 
%                            each online EM runs.
%           - sigmawNewVec = updated vector to store measurement noise in  
%                            each online EM runs.

% AUTHORS  : Xi Chen and Yunpeng Li
%            @ECE Department, McGill University
%            http://www.tsp.ece.mcgill.ca/Networks/index.html
% DATE     : Jun 29, 2011


%1st dimension denotes different tests, 2nd dimension
%is different targets/particles, in single target's case, 
%3rd dimension is the cartesion coordinates)
L = size(xBuffer,3);
N = size(xBuffer,2);
gamma = sqrt(L/counter);
psiNewVec = zeros(N,4);

%loop on different instance, 1:N
for i = 1:N
    %compute sufficient statsitcs
    xBufferOneInstance(:,:) = xBuffer(1:2,i,:);

    xDiff = xBufferOneInstance(:,2:end) - xBufferOneInstance(:,1:end-1);

    YParticlesWithoutPhi = simData(...
        xBufferOneInstance,linkLocs,linkDist,sigma);
    %Calculate sufficient statistics
    psiNewVec(i,1) = (norm(xDiff,2))^2;
    psiNewVec(i,2) = (norm(deltaYBuffer' - phi *...
        YParticlesWithoutPhi, 'fro'))^2;
    psiNewVec(i,3) = sum(sum((YParticlesWithoutPhi.*(deltaYBuffer'))));
    psiNewVec(i,4) = (norm(YParticlesWithoutPhi,'fro'))^2; 
end

if nnz(psiNewVec(:,4)) ~= 0
    psiVec = (1-gamma) * psiVec + gamma * 1/N * sum(psiNewVec);
    phi = psiVec(3) / psiVec(4);
else
    psiVec(:,1:2) = (1-gamma) * psiVec(:,1:2) +...
        gamma * 1/N * sum(psiNewVec(:,1:2));
end

sigma_v = sqrt(psiVec(1) / (2*(L-1)));
if sigma_v == 0
    sigma_v = eps;
end
numLink = size(deltaYBuffer, 2);
sigma_w = sqrt(psiVec(2) / (numLink * L));

phiNewVec = zeros(1,size(phiVec,2)+1);
phiNewVec(1:end-1) = phiVec;
phiNewVec(end) = phi;

sigmavNewVec = zeros(1,size(sigmavVec,2)+1);
sigmavNewVec(1:end-1) = sigmavVec;
sigmavNewVec(end) = sigma_v;

sigmawNewVec = zeros(1,size(sigmawVec,2)+1);
sigmawNewVec(1:end-1) = sigmawVec;
sigmawNewVec(end) = sigma_w;
iteration = 1/gamma^2
true;