function [ logLikelihood ] = computepositionlikelihooddiscretetime( PositionData, isCorrelatedModel, isForPartialMl, asympMeanPosInd, driftRate, posCov, stationaryPosCov, jumpTimes, trajectoryMatBinInds, endTime, maxNTransitions, nStates, nParticles, lateStartTime )
    % Compute likelihood over multiple successive intervals of time for
    % subsequences of discrete position observations from the
    % Ornstein-Uhlenbeck model.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % posCov - the total covariance, i.e. in the 'continuously observed'
    % model, the product of the observation covariance and the jitter
    % covariance. We assume spatial dimensions are independent, so this is
    % a 2*nStates*nParticles matrix.
    %
    % Note: spatial dimensions are assumed to be independent.
    %
    % We assume the first position observation coincides with the
    % initialisation time, t_0 (startTime) (unless using the 'late start
    % time').
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 24.10.2014 - created.
    % 26.10.2014 - finished coding functionality.
    % 01.11.2014 - removed observations at t_0 (startTime) from
    % consideration, as per our convention.
    % 04.11.2014 - corrected the way we implement this: we were still using
    % the observation at t_0.
    % 11.11.2014 - noticed an error in the way we compute the squared
    % distances from the mean: corrected. Refactored computation of
    % consecutiveDistances; is now a field of PositionData. Combined
    % spatial dimensions of sampling statistics to reduce calls of
    % accumarray.
    % 11.12.2014 - replaced coordsTrajectory with posTimeVec.
    % 22.01.2015 - added the option of supplying a start time, needed for
    % computing likelihoods on data intervals that do not start at or
    % before the first position observation.
    % 19.02.2015 - renamed 'fixedStartTime' 'lateStartTime'. Changed the
    % way we compute the likelihood of the initial position observation. We
    % were giving it a probability of 1 in all cases; now we use the
    % stationary distribution for the position process with instantaneous
    % error variance.
    % 27.02.2015 - corrected a bug: we were not squaring the 'squared
    % distances' (bracket in wrong place).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 10,
        error('Insufficient arguments!')
    end
    if nargin < 11,
        maxNTransitions = size(trajectoryMatBinInds, 1);
    elseif size(trajectoryMatBinInds, 1) ~= maxNTransitions,
        error('Size of trajectoryMatBinInds does not match maxNTransitions!')
    end
    if nargin < 12,
        nStates = size(asympMeanPosInd, 1);
    elseif size(asympMeanPosInd, 1) ~= nStates,
        error('Size of asympMeanPosInd does not match nStates!')
    end
    if nargin < 13,
        nParticles = size(asympMeanPosInd, 2);
    elseif size(asympMeanPosInd, 2) ~= nParticles,
        error('Size of asympMeanPosInd does not match nParticles!')
    end
    if (nargin == 14) && (lateStartTime > PositionData.posTimeVec(1)),
        haveLateStartTime = true;
    else
        haveLateStartTime = false;
    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 isCorrelatedModel,
        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
    else
        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
    end
    if ~isfield(PositionData, 'posTimeVec'),
        error('posTimeVec must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'consecutiveDistances'),
        error('consecutiveDistances must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'linearisedPositionTrajectory'),
        error('linearisedPositionTrajectory must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'gridCentreCoordsTrajectory'),
        error('gridCentreCoordsTrajectory 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

    % Most of this code is directed at computing the sum of squared
    % distances of position observations from asymptotic means within the
    % intervals between jumps for each particle.
    
    % Strip data.
    if haveLateStartTime,
        validObsTimeBinInds = PositionData.posTimeVec >= lateStartTime;
        % If no position observations coincide with the fixed startTime,
        % which is after the first observation, we must include the
        % "initial" position, the one immediately before the start time.
        if ~any(PositionData.posTimeVec == lateStartTime),
            validObsTimeBinInds = validObsTimeBinInds | [validObsTimeBinInds(2:end); false];
        end
        PositionData.posTimeVec = PositionData.posTimeVec(validObsTimeBinInds);
        PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds);
        PositionData.gridCentreCoordsTrajectory = PositionData.gridCentreCoordsTrajectory(validObsTimeBinInds, :);
        PositionData.consecutiveDistances = PositionData.consecutiveDistances(validObsTimeBinInds(2:end), :);
    end
    validObsTimeBinInds = PositionData.posTimeVec < endTime;
    PositionData.posTimeVec = PositionData.posTimeVec(validObsTimeBinInds);
    PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds);
    PositionData.gridCentreCoordsTrajectory = PositionData.gridCentreCoordsTrajectory(validObsTimeBinInds, :);

    nObs = length(PositionData.posTimeVec);

    % Get position transformation factors...
    posMeanFactorMatInds = bsxfun(@plus, PositionData.validDiscretePositionsVector, (permute(asympMeanPosInd, [3, 1, 2]) - 1) .* PositionData.nValidDiscretePositions);
    transformationFactors = reshape(PositionData.spaceTransformationFactor(posMeanFactorMatInds(:)), [PositionData.nValidDiscretePositions, 1, nStates, nParticles]);
    % ...and real space coords of asymptotic means.
    asympMeanGridCentreCoords = reshape(PositionData.positionIndToGridCentreCoordsMap(asympMeanPosInd(:), :)', [1, 2, nStates, nParticles]);
    
    transformedSupportPoints = bsxfun(@times, bsxfun(@minus, PositionData.positionIndToGridCentreCoordsMap, asympMeanGridCentreCoords), transformationFactors);

    % Correlated model.
    if isCorrelatedModel,
        logPositionProbs = mymvnpdfmultimulti(transformedSupportPoints, zeros([PositionData.nValidDiscretePositions, 2]), posCov, true);
    % Uncorrelated model.
    else
        covMat = zeros([2, 2, nStates, nParticles]);
        covMat(1, 1, :, :) = permute(posCov(1, :, :), [1, 4, 2, 3]);
        covMat(2, 2, :, :) = permute(posCov(2, :, :), [1, 4, 2, 3]);
        logPositionProbs = mymvnpdfmultimulti(transformedSupportPoints, zeros([PositionData.nValidDiscretePositions, 2]), covMat, true);
    end
    % Normalise.
    logPositionProbs = normaliselogdistributionsmatrix(logPositionProbs, 1);
    
    logLikelihood = logPositionProbs(PositionData.linearisedPositionTrajectory, :, :);
    if (~isequal(size(logLikelihood), [nObs, nStates, nParticles])) && (~isequal(size(logLikelihood), [nObs, nStates])),
        size(logLikelihood)
        [nObs, nStates, nParticles]
        error('logLikelihood not this size?')
    end
end
