function jitterCov = samplecovariancematrix( PriorParams, nStates, nParticles, nObsEachState, distancesFromMean, stateEachObsTime )
    % (For discrete time position model only at the moment): sample
    % correlated anisotropic covariance matrix.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % distancesFromMean - nObs by 2 by nParticles.
    %
    % nObsEachState - 1 by nStates by nParticles.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 11.04.2015 - created.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    nObs = size(distancesFromMean, 1);
    
    if nObs < 1,
        jitterCov = myiwishrndmultimulti(repmat(PriorParams.jitterCovPriorScaleMat, [1, 1, nStates, nParticles]), repmat(PriorParams.jitterCovPriorDof, [1, 1, nStates, nParticles]));
        return;
    end
    
%     sampleSumSquaredDeviations = zeros([2, 2, nStates, nParticles]);

    % This will store the squared deviations and cross-deviations.
    squaredDeviations = zeros([3 * nObs, nParticles]);
    % Squared term (each dimension).
    squaredDeviations(1:(2 * nObs), :) = reshape(distancesFromMean .^ 2, [2 * nObs, nParticles]);
    % Cross-product term.
    squaredDeviations((2 * nObs + 1):(3 * nObs), :) = permute(distancesFromMean(:, 1, :) .* distancesFromMean(:, 2, :), [1, 3, 2]);

    dimensionInds = [ones([nObs, nParticles]); 2 .* ones([nObs, nParticles]); ones([nObs, nParticles])];
    whichProductInd = [ones([nObs, nParticles]); 2 .* ones([2 * nObs, nParticles])];
    paticleInds = repmat(1:nParticles, [3 * nObs, 1]);
    stateInds = repmat(stateEachObsTime, [3, 1]);
    ssIndexingMat = [dimensionInds(:), whichProductInd(:), stateInds(:), paticleInds(:)];
    
    sampleSumSquaredDeviations = accumarray(ssIndexingMat, squaredDeviations(:), [2, 2, nStates, nParticles]);
    sampleSumSquaredDeviations(2, 1, :, :) = sampleSumSquaredDeviations(1, 2, :, :);
    
    covPosteriorScaleMat = bsxfun(@plus, PriorParams.jitterCovPriorScaleMat, sampleSumSquaredDeviations);
    covPosteriorDofs = permute(nObsEachState + PriorParams.jitterCovPriorDof, [1, 4, 2, 3]);

    jitterCov = myiwishrndmultimulti(covPosteriorScaleMat, covPosteriorDofs);
end