function sampledMat = samplefromgeneraliseddirichlet( priorDirichletParams, transitionCounts, firstArrivals )
    % Sample vectors from multiple Generalised 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.
    %
    % firstArrivals - 3-D logical or binary matrix of same size as
    % transitionCounts. Values are the first arrival indicators: (i, j)th
    % value (some page) is true if the first appearance of state j follows
    % an appearance of state i, and is false otherwise. For states to be
    % labelled in order of appearance, we assume that the (i, j)th value is
    % false for all i >= j.
    %
    % See Wong (1998) for more details on this sampling procedure.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 06.10.2014 - copied from old version. Allowed for different number of
    % rows and columns, as when sampling transition rates for generator
    % matrices.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if size(transitionCounts) ~= size(firstArrivals),
        error('transitionCounts and firstArrivals must be matrices of the same size and dimensions!')
    end
    dims = size(transitionCounts);
    m = dims(1);
    n = dims(2);
    if ~isequal(size(priorDirichletParams), [1, n]),
        error('priorDirichletParams must be a row vector with the same number of columns as transitionCounts!')
    end
    if any(any(any(priorDirichletParams <= 0))),
        error('priorDirichletParams must be all positive!')
    end
    if any(any(any(~islogical(firstArrivals) | ((firstArrivals ~= 0) & (firstArrivals ~= 1))))),
        error('firstArrivals must be a logical or binary matrix!')
    end
    
    if length(dims) == 3,
        p = dims(3);
    else
        p = 1;
    end

    % First step (form 'beta' parameters in Wong (1998)).
    b = zeros([m, n, p]);
    b(:, n, :) = firstArrivals(:, n, :) + ones([m, 1, p]);
    for i = (n - 1):-1:1,
        b(:, i, :) = bsxfun(@plus, firstArrivals(:, i, :), priorDirichletParams(1, i + 1)) + transitionCounts(:, i + 1, :) - firstArrivals(:, i + 1, :) + b(:, i + 1, :);
    end
    
    % Second step: sampling.
    a = bsxfun(@plus, priorDirichletParams, transitionCounts) - firstArrivals;

    sampledMat = zeros([m, n, p]);
    sampledMat(:, 1, :) = betarnd(a(:, 1, :), b(:, 1, :));
    tally = sampledMat(:, 1, :);
    for i = 2:n,
        sampledMat(:, i, :) = betarnd(a(:, i, :), b(:, i, :)) .* (1 - tally);
        tally = tally + sampledMat(:, i, :);
    end
    sampledMat = bsxfun(@rdivide, sampledMat, tally);
end