function varargout = samplediscretetimespikerate( ModelSpec, spikeRateGammaHyperparams, updateStepTime, nSpikeTrains, nStates, nParticles, stateOccupancyDurations, spikeCountsMat )
    % Compute the posterior sampling distribution parameters for each
    % neuron's firing rate conditional on (true) state, and sample, for
    % multiple streams (particles).
    %
    % This is the discrete time version. In the continuous time version, we
    % will have full spike times in a cell array, rather than the discrete
    % bins indexing rows of discreteSpikeTrainsMat.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logSpikeRate - a 3-D matrix with (true) states indexing rows,
    % neurons indexing columns and streams (particles) indexing pages. We
    % may be sampling for a subset of the total number of streams.
    %
    % logNonspikeRate - for Bernoulli model only. Same size and
    % dimensions as logSpikeRate. When this is not needed (e.g. when
    % using the discrete Poisson model), only one output needs to be
    % extracted from this function when it is called.
    %
    % spikeRateGammaHyperparams - matrix with first two dimensions the same
    % size as logSpikeRate, and 2 pages. Each page contains one of the
    % two parameters of the prior sampling distribution. This is Beta for
    % the Bernoullli model.
    %
    % stateOccupancyDurations - a 2-D matrix with true states indexing
    % rows and streams (particles) indexing columns. Contains the duration
    % in seconds spend in each state for each stream.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 09.04.2013 - created. Copied code for Bernoulli (thresholded spike
    % presence) from old version. Made allowance for discrete Poisson model
    % in future.
    % 10.04.2013 - introduced logNonspikeRate output for Bernoulli model.
    % - changed size of output to have exactly nPages (those streams we are
    % sampling for), and not to fit output matrices into the correct
    % 'asbolute' pages of a preallocated structure. This may help with
    % memory. Also removed the now redundant absolutePageInds.
    % 11.04.2013 - changed structure of posMeans and posCovs to have size
    % nPages particles dimensions rather than 'nAbsolutePages'. This may
    % help with memory. Also removed the now redundant absolutePageInds.
    % - removed posCovs from input, removed references to absolute indices.
    % 01.05.2013 - found and corrected a bug in the Bernoulli model
    % sampling code. The second beta parameter posterior calculation was
    % not being performed correctly.
    % 02.05.2013 - added stateOccupancyDurations to input. This is
    % available (for free) when we are also sampling the position
    % distributions.
    % 24.06.2013 - implemented Poisson counts (discrete time) model
    % sampling.
    % 27.06.2013 - trueStateTrajectory now includes the initial state
    % separately (changed the model structure), i.e. it now has one more
    % time bin than the observations.
    % - added alternative code for computing the statistics when no states
    % have been sampled, as in when we wish to sample these parameters
    % directly from the prior.
    % 17.07.2013 - renamed function
    % samplediscretetimefiringrateparameters.
    % 18.07.2013 - made stateOccupancyDurations a compulsory input, so
    % that we can ensure it is the same across functions that use it.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 6,
        error('Insufficient arguments!')
    end
    if updateStepTime < 0,
        error('Invalid updateStepTime!')
    end
    if updateStepTime > 0,
        if nargin < 8,
            error('Insufficient arguments!')
        end
        if ~isequal(size(stateOccupancyDurations), [nStates, nParticles]),
            error('stateOccupancyDurations must have nStates rows and nParticles columns!');
        end
        if size(spikeCountsMat, 1) ~= nStates,
            error('spikeCountsMat must have nStates rows, nSpikeTrains columns and nParticles pages!');
        end
        if size(spikeCountsMat, 2) ~= nSpikeTrains,
            error('spikeCountsMat must have nStates rows, nSpikeTrains columns and nParticles pages!');
        end
        if size(spikeCountsMat, 3) ~= nParticles,
            error('spikeCountsMat must have nStates rows, nSpikeTrains columns and nParticles pages!');
        end
    end
    if ~isequal(size(spikeRateGammaHyperparams), [1, 2]),
        error('spikeRateGammaHyperparams must have two elements, containing the two prior parameters for every state/neuron firing rate!');
    end
    if any(any(any(spikeRateGammaHyperparams))) <= 0,
        error('spikeRateGammaHyperparams must be all positive!')
    end

    % For discrete-time thresholded spike data: discreteSpikeTrainsMat is a
    % logical matrix.
    if ModelSpec.spikeTrainModelIndicator == 1,
        % updateStepTime being zero means we are sampling from the prior.
        if updateStepTime == 0,
            spikeCountsMat = zeros([nStates, nSpikeTrains, nParticles]);
            nonspikeCountsMat = zeros([nTrueStates, nSpikeTrains, nParticles]);
        else
            nonspikeCountsMat = bsxfun(@minus, permute(stateOccupancyDurations, [1, 3, 2]), spikeCountsMat);
        end
%         if nTimeBins > 0,
%             spikesByStates = bsxfun(@times, discreteSpikeTrainsMat(1:nTimeBins, :), permute(trueStateTrajectory(2:end, :), [1, 3, 2]));
%             spikeCountsMat = histc(spikesByStates, 1:nTrueStates, 1);
%             nonspikeCountsMat = bsxfun(@minus, permute(stateOccupancyDurations, [1, 3, 2]), spikeCountsMat);

        % This is the case when we sample from the prior - before sampling
        % any state trajectories.
%         else
%             nSpikeTrains = size(discreteSpikeTrainsMat, 2);
%             nParticles = size(trueStateTrajectory, 2);
%             spikeCountsMat = zeros([nTrueStates, nSpikeTrains, nParticles]);
%             nonspikeCountsMat = zeros([nTrueStates, nSpikeTrains, nParticles]);
%         end
        logSpikeRate = betarnd(spikeCountsMat + spikeRateGammaHyperparams(1), nonspikeCountsMat + spikeRateGammaHyperparams(2));
        logSpikeRate(logSpikeRate > (1 - 1e-14)) = 1 - 1e-14;
        logSpikeRate(logSpikeRate < (1e-14)) = 1e-14;
        logNonspikeRate = log(1 - logSpikeRate);
        logSpikeRate = log(logSpikeRate);
        
        varargout = {logSpikeRate, logNonspikeRate};
        
    % For discrete-time count spike data: discreteSpikeTrainsMat is a
    % matrix of spike counts in each time bin.
    elseif ModelSpec.spikeTrainModelIndicator == 2,
        % updateStepTime being zero means we are sampling from the prior.
        if updateStepTime == 0,
            spikeCountsMat = zeros([nStates, nSpikeTrains, nParticles]);
            stateOccupancyDurations = zeros([nStates, nParticles]);
        end
%         nSpikeTrains = size(discreteSpikeTrainsMat, 2);
%         nParticles = size(trueStateTrajectory, 2);
%         if nTimeBins > 0,
%             spikeCountsMat = zeros([nTrueStates, nSpikeTrains, nParticles]);
%             for iParticle = 1:nParticles,
%                 spikesVec = reshape(discreteSpikeTrainsMat(1:nTimeBins, :), [nTimeBins * nSpikeTrains, 1]);
%                 indsMat = [repmat(trueStateTrajectory(2:end, iParticle), [nSpikeTrains, 1]), reshape(repmat(1:nSpikeTrains, [nTimeBins, 1]), [nTimeBins * nSpikeTrains, 1])];
%                 spikeCountsMat(:, :, iParticle) = accumarray(indsMat, spikesVec, [nTrueStates, nSpikeTrains]);
%             end
        % This is the case when we sample from the prior - before sampling
        % any state trajectories.
%         else
%             spikeCountsMat = zeros([nTrueStates, nSpikeTrains, nParticles]);
%             stateOccupancyDurations = zeros([nTrueStates, nParticles]);
%         end
%         stateOccupancyDurationsBySpikeTrain = repmat(permute(stateOccupancyDurations, [1, 3, 2]), [1, nSpikeTrains, 1]);
%         scalePosteriorHyperparam = spikeRateGammaHyperparams(2) ./ (permute(stateOccupancyDurations, [1, 3, 2]) .* spikeRateGammaHyperparams(2) + 1);
        scalePosteriorHyperparam = 1 ./ repmat(spikeRateGammaHyperparams(2) + permute(stateOccupancyDurations, [1, 3, 2]), [1, nSpikeTrains, 1]);
        logSpikeRate = log(gamrnd(spikeCountsMat + spikeRateGammaHyperparams(1), scalePosteriorHyperparam));
        varargout = {logSpikeRate};
    end

end