function sampledMat = samplefromdirichlet( priorDirichletParams, transitionCounts, nStates, nDists, nPages )
    % Sample vectors from multiple Dirichlet sampling distributions.
    % Vectors to be sampled are row vectors. Then different distributions
    % index rows and pages. Rows correspond to different source states of a
    % transition matrix. Pages correspond to different 'streams' or
    % particles we wish to produce samples for. Parameters supplied are the
    % relevant summary statistics. When these are not trivial, we are
    % sampling from a posterior distribution, in, say, a Gibbs sampler
    % sweep.
    % We assume a fixed dimension for all distributions (number of
    % probabilities to sample for each).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % priorDirichletParams is a row vector 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).
    %
    % transitionCounts - 3-D matrix of observed transition counts. The
    % first two dimensions must be of the same size as
    % priorDirichletParams.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 06.10.2014 - copied from earlier version. Now allowing different
    % number of rows (distributions to sample per particle) from columns
    % (states). This is to facilitate sampling of state transitions as part
    % of a generator matrix.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 2,
        error('Insufficient arguments!')
    end
    if nargin < 3,
        nStates = size(transitionCounts, 2);
    end
    % Number of distributions to sample.
    if nargin < 4,
        nDists = size(transitionCounts, 1);
    end
    if nargin < 5,
        nPages = size(transitionCounts, 3);
    else
        if size(transitionCounts, 3) ~= nPages,
            error('nPages does not match the number of pages implied by transitionCounts!')
        end
    end
    if size(priorDirichletParams, 2) ~= nStates,
        error('priorDirichletParams must have the same number of columns as transitionCounts!')
    end
    if any(any(any(priorDirichletParams <= 0))),
        error('priorDirichletParams must be all positive!')
    end
    
    posteriorParameters = bsxfun(@plus, priorDirichletParams, transitionCounts);
    sampledMat = gamrnd(posteriorParameters, ones([nDists, nStates, nPages]));
    sampledMat = bsxfun(@rdivide, sampledMat, sum(sampledMat, 2));
end