function varargout = computepositionposterior( EstModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, timePointsVec, nTimePoints, logInitialStateDist, logTransitionMat, nStates, startTime, endTime, logSmoothedStatePost )
    % This function uses an HMM approximation to the contiuous time model
    % to compute smoothing posteriors over position on a discrete time grid
    % (given spike train observatons).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % PositionData - this is required for some of the fields it contains.
    % No position observations are required.
    %
    % logLikelihood - spike trains only likelihood on each discrete time
    % point. Optional as we can compute it here.
    %
    % logInitialStateDist - to be supplied separately from ModelParams, and
    % appropriate to the interval of data (e.g. the invariant distribution
    % of the process).
    %
    % logTransitionMat - to be supplied separately from ModelParams.
    % Required only for the exact approach to computing the smoothed
    % posteriors by imposing a time grid on the data and considering the
    % model as an HMM. Must be computed (from the generator matrix) for the
    % time discretisation, which must also be specified (time interval must
    % be known).
    %
    % logPosPost - 3-D matrix. Rows are indexed by time point; columns by
    % position (discrete support points); pages by spatial dimension.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 04.02.2015 - created, coded.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 8,
        error('Insufficient arguments!')
    end
    if nargin < 11,
        nStates = length(EstModelParams.logInitialStateDist);
    end
    if nargin < 12,
        startTime = DataSpec.startTime;
    end
    if nargin < 13,
        endTime = DataSpec.endTime;
    end
    if isempty(timePointsVec),
        % Setup time discretisation.
        if isempty(nTimePoints),
            error('At least one of timePointsVec or nTimePoints must be nonempty!')
        end
        timeStep = (endTime - startTime) / nTimePoints;
        timePointsVec = (startTime:timeStep:endTime)';
        % We can't have a discretisation point actually AT the end time,
        % because we exclude data at that point by convention.
        timePointsVec = timePointsVec(1:(end - 1));
    else
        if size(timePointsVec, 2) ~= 1,
            error('timePointsVec must be a column vector!');
        end
        nTimePoints = length(timePointsVec);
        timeStep = (endTime - startTime) / nTimePoints;
    end

    % Must now also compute the transition matrix, knowing the time
    % step.
    if (nargin < 10) || isempty(logTransitionMat),
        logTransitionMat = log(matrixexponentialmulti(timeStep .* EstModelParams.generatorMat));
    end

    if ModelSpec.positionModelIndicator > 0,
        error('positionModelIndicator must be 0 (turn off position model in decoding data ModelSpec)!')
    end

    % We will not use augmented states.
    ModelSpec.useAugStateModel = false;
    
    if nargin < 14,
        [~, logSmoothedStatePost] = decodestate( EstModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, 1, startTime, endTime, [], nTimePoints, logInitialStateDist, logTransitionMat );
    end

    % Get invariant position transition probabilities.
    factorMatInds = bsxfun(@plus, PositionData.validDiscretePositionsVector, (permute(EstModelParams.asympMeanPosInd, [2, 1]) - 1) .* PositionData.nValidDiscretePositions);
    transformationFactors = reshape(PositionData.spaceTransformationFactor(factorMatInds(:)), [PositionData.nValidDiscretePositions, 1, nStates]);
    clear factorMatInds;
    asympMeanGridCentreCoords = reshape(PositionData.positionIndToGridCentreCoordsMap(EstModelParams.asympMeanPosInd(:), :)', [1, 2, nStates]);
    transformedSupportPoints = bsxfun(@times, bsxfun(@minus, PositionData.positionIndToGridCentreCoordsMap, asympMeanGridCentreCoords), transformationFactors);
    % Correlated model.
    if ModelSpec.positionJitterModelIndicator == 3,
        logPositionProbs = mymvnpdfmultimulti(transformedSupportPoints, zeros([PositionData.nValidDiscretePositions, 2]), EstModelParams.posCov, true);
    else
    % Uncorrelated model.
        covMat = zeros([2, 2, nStates]);
        covMat(1, 1, :) = permute(EstModelParams.posCov(1, :), [1, 3, 2]);
        covMat(2, 2, :) = permute(EstModelParams.posCov(2, :), [1, 3, 2]);
        logPositionProbs = mymvnpdfmultimulti(transformedSupportPoints, zeros([PositionData.nValidDiscretePositions, 2]), covMat, true);
    end
    % Normalise.
    logPositionProbs = normaliselogdistributionsmatrix(logPositionProbs, 1);

    logPosPost = sumlogprobmat(bsxfun(@plus, permute(logPositionProbs, [3, 1, 2]), permute(logSmoothedStatePost, [1, 3, 2])), 3);
    [~, mapPosTrajectory] = max(logPosPost, [], 2);
    mapPosTrajectory = [timePointsVec, mapPosTrajectory];
    [ logYPosPost, logXPosPost ] = marginalisepositionlogposteriors( logPosPost, PositionData );

    if nargout == 3,
        varargout = {logPosPost, logYPosPost, logXPosPost};
    elseif nargout == 4,
        varargout = {logPosPost, logYPosPost, logXPosPost, [timePointsVec, mapPosTrajectory]};
    elseif nargout == 5,
        varargout = {logPosPost, logYPosPost, logXPosPost, [timePointsVec, mapPosTrajectory], logTransitionMat};
    end
end
