function [ asympMeanPosInd, nObsEachState ] = sampleasymptoticmeanposdiscretetime( PriorParams, ModelSpec, PositionData, updateStepTime, driftRate, posCov, nStates, nParticles, stateEachObsTime, validObsTimeBinInds )
    % Sample the asymptotic mean from the full conditional given sampled
    % MJP trajectory.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Note: we assume that the first position observation occurs at
    % initialisation time, t_0 (startTime).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 27.10.2014 - created, coded functionality.
    % 03.11.2014 - discovered an error in calculations. Corrected approach
    % to computing sampling distribution, now based on first computing the
    % posterior mean position (in full 2-D to ensure posterior mean is
    % 'valid').
    % 08.11.2014 - refactored compuation of sampling statistics out to
    % computesamplingstatistics.m.
    % 17.11.2014 - refactored out validObsTimeBinInds.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 10,
        error('Insufficient arguments!')
    end
    if updateStepTime < 0,
        error('Invalid updateStepTime!')
    end
    if ~isfield(ModelSpec, 'positionJitterModelIndicator'),
        error('positionJitterModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(PositionData, 'linearisedPositionTrajectory'),
        error('linearisedPositionTrajectory must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'spaceTransformationFactor'),
        error('spaceTransformationFactor must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'positionIndToGridCentreCoordsMap'),
        error('positionIndToGridCentreCoordsMap must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'nValidDiscretePositions'),
        error('nValidDiscretePositions must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'validDiscretePositionsVector'),
        error('validDiscretePositionsVector must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'sandMats'),
        error('sandMats must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'resShifterMats'),
        error('resShifterMats must be a field of PositionData!')
    end
    if size(driftRate, 1) ~= nStates,
        error('Size of driftRate does not match nStates!')
    end
    if size(driftRate, 2) ~= nParticles,
        error('Size of driftRate does not match nParticles!')
    end
    if (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
        if size(posCov, 2) ~= nStates,
            error('Size of posCov does not match nStates!')
        end
        if size(posCov, 3) ~= nParticles,
            error('Size of posCov does not match nParticles!')
        end
    elseif (ModelSpec.positionJitterModelIndicator == 3),
        if size(posCov, 3) ~= nStates,
            error('Size of posCov does not match nStates!')
        end
        if size(posCov, 4) ~= nParticles,
            error('Size of posCov does not match nParticles!')
        end

    end

    % Strip data.
    PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds);
    nObs = size(PositionData.linearisedPositionTrajectory, 1);

    if ~isequal(size(stateEachObsTime), [nObs, nParticles]),
        error('stateEachObsTime must be a matrix with one row for each position observation (in this data interval) and one column for each particle!')
    end
    if sum(validObsTimeBinInds) ~= nObs,
        error('validObsTimeBinInds must be a binary column vector with nObs ones!')
    end
   
    % Setup prior probabilities.
    logPosteriorProbs = repmat(PriorParams.asympMeanPriorDist, [1, nStates, nParticles]);
    if nObs < 1,
        % If too few observations in this interval, sample from uniform prior.
        asympMeanPosInd = samplefromcategoricaldistributions(permute(logPosteriorProbs, [2, 3, 1]), 3);
        return;
    end
    logPosteriorProbs = log(logPosteriorProbs);

    % State at each position observation, excluding the initial
    % observation, as per our convention.
%     currentStateEachObsTime = stateEachObsTime(2:end, :);

    % First we need to compute relevant statistics.

    % Count of visits to each position whilst in each state, including
    % initial position but excluding last position (we already have the
    % corresponding statistics excluding initial position but including
    % last position), posCountEachState.
%     previousPosCountEachState = accumarray([repmat(PositionData.linearisedPositionTrajectory(1:(end - 1)), [nParticles, 1]), currentStateEachObsTime(:), currentTimeParticleInds(:)], 1, [PositionData.nValidDiscretePositions, nStates, nParticles]);
    particleInds = repmat((1:nParticles), [nObs, 1]);
    posCountEachState = accumarray([repmat(PositionData.linearisedPositionTrajectory, [nParticles, 1]), stateEachObsTime(:), particleInds(:)], 1, [PositionData.nValidDiscretePositions, nStates, nParticles]);
    nObsEachState = sum(posCountEachState, 1);
    clear currentStateEachObsTime;

    % Identify state/particle conbinations for which we have observations
    % and hence have posterior not equal to the prior.
    sampleFromPostBinInds = permute(nObsEachState > 0, [2, 3, 1]);
    sampleFromPostBinIndsBig = repmat(nObsEachState > 0, [PositionData.nValidDiscretePositions, 1, 1]);
    nToSampleFromPost = sum(sum(sampleFromPostBinInds));
    
    if nToSampleFromPost > 0,
        % Compute the posterior mean discrete position.
        posCountEachState = reshape(posCountEachState(sampleFromPostBinIndsBig), [PositionData.nValidDiscretePositions, nToSampleFromPost]);
%         previousPosCountEachState = reshape(previousPosCountEachState(sampleFromPostBinIndsBig), [PositionData.nValidDiscretePositions, nToSampleFromPost]);
        previousPosCountEachState = zeros(size(posCountEachState));
        if nStates == 1,
            postMean = computeposteriormeanasympmean( [], 1, PositionData, driftRate(sampleFromPostBinInds)', posCountEachState, previousPosCountEachState, nToSampleFromPost );
        else
            postMean = computeposteriormeanasympmean( [], 1, PositionData, driftRate(sampleFromPostBinInds), posCountEachState, previousPosCountEachState, nToSampleFromPost );
        end
        clear posCountEachState previousPosCountEachState;

        % Evaluate posterior pdf at transformed positions corresponding to
        % the valid discrete maze positions (sample space for asymptotic
        % mean). Transformation is relative to the posterior mean.

        factorMatInds = bsxfun(@plus, PositionData.validDiscretePositionsVector, (postMean' - 1) .* PositionData.nValidDiscretePositions);
        transformationFactors = reshape(PositionData.spaceTransformationFactor(factorMatInds(:)), [PositionData.nValidDiscretePositions, 1, nToSampleFromPost]);
        clear factorMatInds;

        postMeanGridCoords = reshape(PositionData.positionIndToGridCentreCoordsMap(postMean, :)', [1, 2, nToSampleFromPost]);
        transformedSupportPoints = bsxfun(@times, bsxfun(@minus, PositionData.positionIndToGridCentreCoordsMap, postMeanGridCoords), transformationFactors);

        % Posterior covariance (stored and used as precision).
        if (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            postPrecision = bsxfun(@rdivide, nObsEachState, posCov);
            % Now state by particle by 2.
            postPrecision = permute(postPrecision, [2, 3, 1]);
            yPostPrecision = postPrecision(:, :, 1);
            xPostPrecision = postPrecision(:, :, 2);
            clear postPrecision;

            % Now Position by particle by 2.
            transformedSupportPoints = permute(transformedSupportPoints, [1, 3, 2]);
            % These are position * state * particle.
            yTransformedSupportPoints = transformedSupportPoints(:, :, 1);
            xTransformedSupportPoints = transformedSupportPoints(:, :, 2);
            clear postMean postMeanGridCoords transformationFactors transformedSupportPoints;

            % Multiply (independent) spatial dimensions to get unnormalised
            % sampling probabilities on valid discrete positions.
            % Probs will go in first dimension.
            % We will do this in log domain to avoid underflow. Note we don't
            % need the 2 * pi term since we are normalising anyway.
            % Initialise with the prior probabilities.
            if nStates == 1,
                yLogPosteriorProbs = bsxfun(@times, yTransformedSupportPoints .^ 2, (-1/2) .* yPostPrecision(sampleFromPostBinInds));
                xLogPosteriorProbs = bsxfun(@times, xTransformedSupportPoints .^ 2, (-1/2) .* xPostPrecision(sampleFromPostBinInds));
            else
                yLogPosteriorProbs = bsxfun(@times, yTransformedSupportPoints .^ 2, (-1/2) .* yPostPrecision(sampleFromPostBinInds)');
                xLogPosteriorProbs = bsxfun(@times, xTransformedSupportPoints .^ 2, (-1/2) .* xPostPrecision(sampleFromPostBinInds)');
            end
            clear yPostPrecision xPostPrecision yTransformedSupportPoints xTransformedSupportPoints;

            logPosteriorProbs(sampleFromPostBinIndsBig) = yLogPosteriorProbs(:) + xLogPosteriorProbs(:);
        elseif ModelSpec.positionJitterModelIndicator == 3,
            % In this case posCov is a matrix.
            sampleFromPostCovMatBinInds = repmat(permute(sampleFromPostBinInds, [3, 4, 1, 2]), [2, 2, 1, 1]);
%             postPosCov = zeros([2, 2, nToSampleFromPost]);
            postPosCov = bsxfun(@rdivide, reshape(posCov(sampleFromPostCovMatBinInds), [2, 2, nToSampleFromPost]), reshape(nObsEachState(permute(sampleFromPostBinInds, [3, 1, 2])), [1, 1, nToSampleFromPost]));
            logPosteriorProbs(sampleFromPostBinIndsBig) = mymvnpdfmultimulti(transformedSupportPoints, zeros([PositionData.nValidDiscretePositions, 2]), postPosCov, true);
        end
    end
    % Normalise.
    logPosteriorProbs = bsxfun(@minus, logPosteriorProbs, sumlogprobmat(logPosteriorProbs, 1));
    % Sample.
    asympMeanPosInd = permute(samplefromcategoricaldistributions(exp(logPosteriorProbs), 1), [2, 3, 1]);

end
