function samplesMat = samplefromcategoricaldistributions(samplingDistributionsMat, distributionDim, dims)
    % Obtain a vector or vectors of sampled indices from categorical
    % distributions supplied in a matrix. One dimension of the matrix
    % contains the probabilities of the sampling distributions. Other
    % dimensions are arbitrary. Supplied with an N dimensional
    % distributions matrix, the function will return an N dimensional
    % matrix (with a singleton dimension) of samples from the categorical
    % distributions detailed in the removed dimension.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % E.g. if the distributions are in columns of a m by n square matrix,
    % the function will return a row vector of indices between 1 and m from
    % each column's distribution. If the matrix is m by n by p, returned is
    % a 1 by n by p matrix of samples.
    %
    % distributionDim indicates which dimension contains the categorical
    % probabilities. This the matrix must be normalised in this dimension.
    %
    % Note the distributions matrix must not be in the log domain.
    %
    % CAUTION: I have removed the check for normalised distributions
    % because numerical precision makes it hard to check when the
    % dimensions of the sampling matrix is not constrained.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 13.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 3,
        dims = size(samplingDistributionsMat);
    end
%     if ~isequal(sum(samplingDistributionsMat, distributionDim), ones([dims(1:(distributionDim - 1)), 1, dims((distributionDim + 1):end)])),
%         error('samplingDistributionsMat must be normalised in dimension distributionDim!')
%     end

    u = rand([dims(1:(distributionDim - 1)), 1, dims((distributionDim + 1):end)]);
    [~, samplesMat] = max(bsxfun(@minus, cumsum(samplingDistributionsMat, distributionDim), u) > 0, [], distributionDim);
%     samplesMat = find(bsxfun(@minus, cumsum(samplingDistributionsMat, distributionDim), u) > 0, 1);
end