function logTransitionMat = samplehmmtransitionmatrix( ModelSpec, transitionCountsMat, firstArrivals, priorDirichletParams, nTrueStates, nPages, nAugmentedStates )
    % This function samples transition matrices for a hidden state process,
    % potentially for multiple streams (e.g. particles) simultaneously.
    % This could be as part of a Gibbs sweep, in which case we are sampling
    % from a full conditional over transition probabilities given sampled
    % state trajectories (one for each stream). Otherwise we are sampling
    % from a prior distribution.
    % We allow for the state space to be 'regular' or 'augmented'. In the
    % former case, the sampling distribution is Dirichlet for each row of
    % each transtion matrix. For the latter it is Generalised Dirichlet.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logTransitionMat - preallocated structure. Should be log
    % probabilities with normalised rows. First two dimensions hold square
    % transition matrices.
    %
    % trueStateTrajectory - a T by H matrix where T is the maximum number
    % of discrete time bins for which any of the concurrent sampling
    % streams (i.e. particles) has state samples, and H is the number of
    % those streams (particles). In each column is a sampled trajectory of
    % (regular, or true - not augmented) states.
    %
    % priorDirichletParams is a matrix of the prior Dirichlet parameters.
    % Must have number of columns equal to the dimension of the
    % distributions (same as the number of columns of transitionCounts and
    % firstArrivals) but can have any number of rows or pages.
    %
    % nTrueStates - for 'regular' states this is the size of the state
    % space; for augmented states it is the number of 'true' states.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 09.04.2013 - created. Copied some code from old version. Allowed for
    % 'regular' or 'augmented' transition matrix sampling.
    % 10.04.2013 - wrote data checks and initialisations.
    % 16.04.2013 - removed logTransitionMat from input, removed references
    % to absolute indices.
    % 17.04.2013 - corrected procedure so that when nTimeBins == 0,
    % transitionCounts is all zeros, as expected.
    % 01.07.2013 - found and corrected a bug: transitionIndexingMat's third
    % column was not always correct due to the use of floor and numerical
    % truncation. We now construct the third column using repmat and
    % reshape.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 3,
        error('Insufficient arguments!')
    end
    if ModelSpec.useAugStateModel,
        if nargin < 7,
            nAugmentedStates = nTrueStates * (nTrueStates + 1) / 2;
        end
    end
    if nargin < 4,
        nTrueStates = size(priorDirichletParams, 1);
    else
        if size(priorDirichletParams, 2) ~= nTrueStates,
            error('priorDirichletParams must have nTrueStates columns!')
        end
    end
    if any(priorDirichletParams < 0),
        error('priorDirichletParams must be all positive!')
    end
    if nargin < 5,
        nPages = size(transitionCountsMat, 3);
    else
        if size(transitionCountsMat, 3) ~= nPages,
            error('transitionCountsMat must have nPages columns!')
        end
    end

    if ModelSpec.useAugStateModel,
        logTransitionMat = samplefromgeneraliseddirichlet( priorDirichletParams, transitionCountsMat, firstArrivals );
        
        % Convert to augmented form.
        logTransitionMat = transformregularmatrixtoaugmented( logTransitionMat, nTrueStates, nPages, nAugmentedStates );
    else
        % If we are only using 'regular' states.
        logTransitionMat = samplefromdirichlet( priorDirichletParams, transitionCountsMat, nTrueStates, nTrueStates, nPages );
    end
    
    % Convert to log domain.
    logTransitionMat = log(logTransitionMat);
end