function [ varargout ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectoryMat, jumpTimes, trajectoryMatBinInds, interjumpIntervals, nTrueStates, updateStepTime, maxNTransitions, nParticles )
    % Given a matrix of true state samples over time, and the observation
    % data, compute the various sample summary quantities necessary for
    % sampling new parameter values, including: the state occupancy
    % durations for each sampled MJP trajectory.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % transitionCountsMat - has one fewer columns than rows: we ignore
    % self-transitions.
    %
    % firstArrivals - has one fewer columns than rows: we ignore
    % self-transitions.
    %
    % stateTrajectoryMat - a T by H matrix where T is the maximum number of
    % jumps for any particle and H is the number of particles. In each
    % column is a sampled trajectory of states, augmented if using
    % augmented state model.
    %
    % AlgorithmParams - struct that must contain the following fields:
    % - timeDiscretisation
    % - discreteTimeBinWidth - only if discrete time algorithm used.
    %
    % PositionData - struct with the following fields:
    % - linearisedPositionTrajectory - the shortest time interval binned
    % position observations, as linear positions.
    % - nValidPositions
    % - startTime - time (seconds) of earliest observation used for
    % analysis.
    % - posObsInterval - the shortest length of time, in secods, between
    % any two position observations we have. Use this to index posTimeVec
    % (and the equivalent matrix of linearised coordinates computed
    % elsewhere) using:
    % index = floor((t - startTime) / posObsInterval) + 1;
    % - for any time t, in seconds, greater than the earliest time
    % possible, startTime.
    % These time bins are the smallest interval between any two successive
    % position observations. This allows us a high precision of observed
    % position at any (arbitrary and continuos) time.
    % - positionObsTimesVec - column vector of the times (in seconds)
    % corresponding to each position observation in
    % linearisedPositionTrajectory.
    %
    % jumpTimes - cell array with one cell for each column of
    % trueStateTrajectoryMat, with each cell being a column vector
    % detailing the actual times (in seconds) for which each
    % stream/particle made the corresponding sample in
    % trueStateTrajectoryMat. For discrete time sampling these cells will
    % all be the same, and the times will be evenly spaced.
    %
    % covMatExtendedStateLinearInds - column vector of linear indices
    % for 2 by 2 covariance matrices, stacked and augmented by state and
    % particle indices. Used for linearly indexing covariance matrices for
    % particular state/particle combinations.
    %
    % extendedStateLinearInds - column vector of the linear indices of the
    % 'extended' states. That is, the linear indices of a 2-D matrix where
    % true states index rows and streams index columns for the state/stream
    % combinations we have samples for.
    %
    % nTrueStates - total number of 'true' states.
    %
    % nValidPositions - total number of valid positions in the discretised
    % environment.
    %
    % nParticles - number of particles for which we have trajectory
    % samples.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 16.10.2014 - copied from old version, analyzestatesamplescell.m.
    % 20.10.2014 - modified parameter names and structures to conform to
    % current version. Added transitionCountsMat and firstArrivals to
    % output for when using augmented states model. Have not checked the
    % position model sampling
    % statistics - for now, we will not use the position model.
    % 21.20.2014 - added transitionCountsMat to output for non-aug states
    % model - needed here also. Also added spikeCountsMat computation to
    % function.
    % 23.10.2014 - changed spikesByParticleVec so that rather than growing
    % this vector we index it. Also introduced possiblity of empty spike
    % trains.
    % 08.11.2014 - added position model parameters' sampling statistics.
    % 11.12.2014 - replaced coordsTrajectory with posTimeVec.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 9,
        error('Insufficient arguments!')
    end
    if nargin < 10,
        maxNTransitions = size(stateTrajectoryMat, 1);
    else
        if size(stateTrajectoryMat, 1) ~= maxNTransitions,
            error('Number of rows implied by stateTrajectoryMat does not match maxNTransitions!')
        end
    end
    if nargin < 11,
        nParticles = size(stateTrajectoryMat, 2);
    else
        if size(stateTrajectoryMat, 2) ~= nParticles,
            error('Number of pages implied by stateTrajectoryMat does not match nParticles!')
        end
    end
    if ~isequal(size(trajectoryMatBinInds), size(stateTrajectoryMat)),
        error('trajectoryMatBinInds and stateTrajectoryMat must be the same size/dimensions!')
    end
    if ~isequal(size(jumpTimes), size(interjumpIntervals)),
        error('jumpTimes and interjumpIntervals must be the same size/dimensions!')
    end
    if size(jumpTimes, 1) ~= nParticles,
        error('jumpTimes must be a cell array with the same number of cells, in a column, as stateTrajectoryMat has columns (streams/particles)!')
    end
    if updateStepTime < 0,
        error('updateStepTime must be positive!')
    end
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'spikeTrainModelIndicator'),
        error('spikeTrainModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'positionModelIndicator'),
        error('positionModelIndicator must be a field of ModelSpec!')
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        if ~isfield(SpikeTrainData, 'nSpikeTrains'),
            error('nSpikeTrains must be a field of SpikeTrainData!')
        end
        if ~isfield(SpikeTrainData, 'spikeTimesArray'),
            error('spikeTimesArray must be a field of SpikeTrainData!')
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        if ~isfield(PositionData, 'nValidDiscretePositions'),
            error('nValidDiscretePositions must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'posTimeVec'),
            error('posTimeVec must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'linearisedPositionTrajectory'),
            error('linearisedPositionTrajectory must be a field of PositionData!')
        end
    end
    if ModelSpec.useAugStateModel,
        requiredNargout = 3;
    else
        requiredNargout = 2;
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        requiredNargout = requiredNargout + 1;
    end
    if ModelSpec.positionModelIndicator > 0,
        requiredNargout = requiredNargout + 6;
    end
    if nargout ~= requiredNargout,
        error('Incorrect outputs!')
    end
    
    % This only happens when we need these structures to use parameter
    % sampling functions for sampling from the prior.
    if maxNTransitions < 1,
        stateOccupancyDurations = zeros([nTrueStates, nParticles]);
        transitionCountsMat = zeros([nTrueStates, nTrueStates - 1, nParticles]);
        if ModelSpec.useAugStateModel,
            firstArrivals = false([nTrueStates, nTrueStates - 1, nParticles]);
            varargout = {stateOccupancyDurations, transitionCountsMat, firstArrivals};
        else
            varargout = {stateOccupancyDurations, transitionCountsMat};
        end
        return;
    end

    % Convert sampled process to 'true' state trajectories - if the one we
    % have is augmented.
    if ModelSpec.useAugStateModel,
        trueStateTrajectoryMat = nan(size(stateTrajectoryMat));
        trueStateTrajectoryMat(trajectoryMatBinInds) = ModelSpec.augStateComponents(stateTrajectoryMat(trajectoryMatBinInds), 1);
    else
        trueStateTrajectoryMat = stateTrajectoryMat;
    end

    % True state occupancy durations - the total time spent in each true
    % state.
    interjumpIntervalsVec = cell2mat(interjumpIntervals);
    particleIndsMat = repmat(1:nParticles, [maxNTransitions, 1]);
    particleIndsMat = particleIndsMat(trajectoryMatBinInds);
    if maxNTransitions > 1,
        indsMat = [trueStateTrajectoryMat(trajectoryMatBinInds), particleIndsMat(:)];
    else
        indsMat = [trueStateTrajectoryMat(trajectoryMatBinInds)', particleIndsMat(:)];
    end
    stateOccupancyDurations = accumarray(indsMat, interjumpIntervalsVec, [nTrueStates, nParticles]);

    % True state transition counts matrix. State transition counts exclude
    % same-state transitions (no transitions to self).
    if maxNTransitions > 1,
        if ModelSpec.modelIndicator == 1,
            transitionIndexingMat = [reshape(trueStateTrajectoryMat(1:(maxNTransitions - 1), :), [(maxNTransitions - 1) * nParticles, 1]), reshape(trueStateTrajectoryMat(2:maxNTransitions, :), [(maxNTransitions - 1) * nParticles, 1]), reshape(repmat(1:nParticles, [maxNTransitions - 1, 1]), [nParticles * (maxNTransitions - 1), 1])];
            transitionCountsMat = accumarray(transitionIndexingMat, 1, [nTrueStates, nTrueStates, nParticles]);
        elseif ModelSpec.modelIndicator == 2,
            % Useful indexing structures.
            diagonalBinInds = repmat(eye(nTrueStates) & true([nTrueStates, nTrueStates]), [1, 1, nParticles]);
            offDiagonalBinInds = true([nTrueStates, nTrueStates, nParticles]) & ~diagonalBinInds;

            trajectoryMatBinIndsAndNext = trajectoryMatBinInds(1:(maxNTransitions - 1), :) & trajectoryMatBinInds(2:maxNTransitions, :);
            pageInds = repmat(1:nParticles, [maxNTransitions - 1, 1]);
            pageInds = pageInds(trajectoryMatBinIndsAndNext);
            transitionIndexingMat = [trueStateTrajectoryMat([trajectoryMatBinIndsAndNext; false([1, nParticles])]), trueStateTrajectoryMat([false([1, nParticles]); trajectoryMatBinIndsAndNext]), pageInds(:)];
            % Should have zeros down main diagonal (no self-transitions).
            transitionCountsMatWithZeros = accumarray(transitionIndexingMat, 1, [nTrueStates, nTrueStates, nParticles]);
            % Remove diagonal.
            % This permute / unpermute operation required to get transition
            % rates into correct positions.
            transitionCountsMatWithZeros = permute(transitionCountsMatWithZeros, [2, 1, 3]);
            transitionCountsMat = reshape(transitionCountsMatWithZeros(offDiagonalBinInds), [nTrueStates - 1, nTrueStates, nParticles]);
            transitionCountsMat = permute(transitionCountsMat, [2, 1, 3]);
        end

        if ModelSpec.useAugStateModel,
            % First arrivals matrix: the matrix of indicators of from which
            % (true) source state did each trajectory first arrive at each
            % (true) destination state.
            firstArrivals = formfirstarrivalsmatrix( trueStateTrajectoryMat, nTrueStates, maxNTransitions, nParticles );
            if ModelSpec.modelIndicator == 1,
                firstArrivalsWithZeros = false([nTrueStates, nTrueStates, nParticles]);
                firstArrivals = permute(firstArrivals, [2, 1, 3]);
                binInds = repmat(eye(nTrueStates) == 0, [1, 1, nParticles]);
                firstArrivalsWithZeros(binInds) = firstArrivals(:);
                firstArrivals = firstArrivalsWithZeros;
                firstArrivals = permute(firstArrivals, [2, 1, 3]);
            end
        end
    else
        if ModelSpec.modelIndicator == 1,
            transitionCountsMat = zeros([nTrueStates, nTrueStates, nParticles]);
            if ModelSpec.useAugStateModel,
                firstArrivals = false([nTrueStates, nTrueStates, nParticles]);
            end
        elseif ModelSpec.modelIndicator == 2,
            transitionCountsMat = zeros([nTrueStates, nTrueStates - 1, nParticles]);
            if ModelSpec.useAugStateModel,
                firstArrivals = false([nTrueStates, nTrueStates - 1, nParticles]);
            end
        end
    end
    
    
    % Spike train model parameters.
    if ModelSpec.spikeTrainModelIndicator > 0,
        if updateStepTime > 0,
            % First get counts of spikes in each time interval for each
            % particle, and for each spike train.
            % Loop over spike trains.
            nTotalJumps = sum(sum(trajectoryMatBinInds));
            spikesByParticleVec = zeros([nTotalJumps * SpikeTrainData.nSpikeTrains, 1]);

            % We add the SMC step time - the upper time limit for this
            % sampling interval - to the end of the vectors of jump times,
            % which we use as interval bounds between which we count
            % spikes.
            jumpTimes = cellfun(@(timesVec) [timesVec; updateStepTime], jumpTimes, 'UniformOutput', false);
            for iTrain = 1:SpikeTrainData.nSpikeTrains,
                if isempty(SpikeTrainData.spikeTimesArray{iTrain}),
                    continue;
                end
                % Returns a cell array, one cell for each particle, of
                % vectors of the counts in each of the time intervals
                % specified for each particle.
                countsThisTrainBetweenSampleTimes = cellfun(@(timeVec) histc(SpikeTrainData.spikeTimesArray{iTrain}, timeVec, 1), jumpTimes, 'UniformOutput', false);
                % Strip off the last interval counts: these occur at (and
                % beyond) the end time.
                countsThisTrainBetweenSampleTimes = cellfun(@(iCountsVec) iCountsVec(1:(end - 1)), countsThisTrainBetweenSampleTimes, 'UniformOutput', false);
                indsVec = ((iTrain - 1) * nTotalJumps + 1):(iTrain * nTotalJumps);
                spikesByParticleVec(indsVec) = cell2mat(countsThisTrainBetweenSampleTimes);
            end

            % Now we must total up the number of spikes observed whilst in
            % each state.
            % We do this via accumarray(). First get the indices.
            binIndsMat = repmat(trajectoryMatBinInds, [1, 1, SpikeTrainData.nSpikeTrains]);
            spikeTrainIndsMat = repmat(permute(1:SpikeTrainData.nSpikeTrains, [1, 3, 2]), [maxNTransitions, nParticles, 1]);
            spikeTrainIndsMat = spikeTrainIndsMat(binIndsMat);
            particleIndsMat = repmat(1:nParticles, [maxNTransitions, 1, SpikeTrainData.nSpikeTrains]);
            particleIndsMat = particleIndsMat(binIndsMat);
            if maxNTransitions == 1,
                indsMat = [repmat(trueStateTrajectoryMat(trajectoryMatBinInds)', [SpikeTrainData.nSpikeTrains, 1]), spikeTrainIndsMat(:), particleIndsMat(:)];
            else
                indsMat = [repmat(trueStateTrajectoryMat(trajectoryMatBinInds), [SpikeTrainData.nSpikeTrains, 1]), spikeTrainIndsMat(:), particleIndsMat(:)];
            end

            spikeCountsMat = accumarray(indsMat, spikesByParticleVec, [nTrueStates, SpikeTrainData.nSpikeTrains, nParticles]);

        else
            spikeCountsMat = zeros([nTrueStates, SpikeTrainData.nSpikeTrains, nParticles]);
        end
    end

    % Position model parameters.
    if ModelSpec.positionModelIndicator > 0,
        % Compute some statistics for position model paramter full
        % conditional distributions: stateEachObsTime,
        % currentTimeParticleInds, nObsEachState, currentPosCountEachState,
        % stateParticleIndsEachObsTime
        
        % Strip data.
        validObsTimeBinInds = PositionData.posTimeVec < updateStepTime;
        PositionData.posTimeVec = PositionData.posTimeVec(validObsTimeBinInds);
        PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds, :);
        
        nObs = length(PositionData.posTimeVec);
    
        % Gives us the index, for each position, of the interval in which it
        % lies, for each particle, in a cell array indexed by particle.
        [~, sampleTimeInds] = cellfun(@(timeVec) histc(PositionData.posTimeVec(validObsTimeBinInds), [timeVec; updateStepTime], 1), jumpTimes, 'UniformOutput', false);
        % Increment state indices by particle for indexing state trajectory
        % matrix. Also transform to column vector.
        sampleTimeInds = cell2mat(cellfun(@(indsVec, iParticle) indsVec + (iParticle - 1) .* maxNTransitions, sampleTimeInds, num2cell((1:nParticles)'), 'UniformOutput', false));
        % State at each position observation.
        stateEachObsTime = reshape(trueStateTrajectoryMat(sampleTimeInds), [nObs, nParticles]);

        % Excluding the initial observation, as per our convention.
        currentStateEachObsTime = stateEachObsTime(2:end, :);

        % For indexing model parameters.
        particleInds = repmat((1:nParticles), [nObs, 1]);
        stateParticleIndsEachObsTime = stateEachObsTime + (particleInds - 1) .* nTrueStates;   
    
        currentTimeParticleInds = particleInds(2:end, :);
        if nObs > 1,
            currentPosCountEachState = accumarray([repmat(PositionData.linearisedPositionTrajectory(2:end), [nParticles, 1]), currentStateEachObsTime(:), currentTimeParticleInds(:)], 1, [PositionData.nValidDiscretePositions, nTrueStates, nParticles]);
            if ModelSpec.modelIndicator == 2,
                previousPosCountEachState = accumarray([repmat(PositionData.linearisedPositionTrajectory(1:(end - 1)), [nParticles, 1]), currentStateEachObsTime(:), currentTimeParticleInds(:)], 1, [PositionData.nValidDiscretePositions, nTrueStates, nParticles]);
            else
                previousPosCountEachState = zeros([0, nTrueStates, nParticles]);
            end
        else
            currentPosCountEachState = zeros([0, nTrueStates, nParticles]);
            previousPosCountEachState = zeros([0, nTrueStates, nParticles]);
        end

        % Total number of observations made whilst in each state, for each
        % particle (not including first observation).
        nObsEachState = sum(currentPosCountEachState, 1);
    end


    if ~ModelSpec.useAugStateModel,
        varargout = {stateOccupancyDurations, transitionCountsMat};
    else
        varargout = {stateOccupancyDurations, transitionCountsMat, firstArrivals};
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        varargout = {varargout{:}, spikeCountsMat};
    end
    if ModelSpec.positionModelIndicator > 0,
        varargout = {varargout{:}, stateEachObsTime, currentTimeParticleInds, stateParticleIndsEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState};
    end
end