function sample = samplefromcategoricaldistribution( samplingDistribution, nSamples, sampleDimension )
    % Obtain a vector of indices of elements sampled from a categorical
    % distribution whose point masses are listed in a vector. Dimensions of
    % input and output can be arbitrary (i.e. vectors in any dimension).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % samplingDistribution - vector of probabilities that must sum to 1.
    % Must not be in log domain.
    %
    % nSamples - size of the sample to obtain - number of times to draw
    % from the distribution.
    %
    % sampleDimension - dimension in which we should return the samples
    % vector.
    %
    % sample - vector in dimension sampleDimension. E.g. if this is 1 and
    % nSamples is 3, we will return a 3 by 1 column vector. This is the
    % case regardless of the dimension of the sampling distribution vector.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 10.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    distributionDims = size(samplingDistribution);
    if (length(distributionDims) > 2) || ((distributionDims(1) > 1) && (distributionDims(2) > 1)),
        error('samplingDistribution must be a vector!')
    end

    nDims = length(distributionDims);
    [~, distributionDimension] = max(distributionDims);
    
    sampleDims = ones([1, max([2, sampleDimension])]);
    sampleDims(sampleDimension) = nSamples;

    randVec = rand(sampleDims);
    
    if sampleDimension == distributionDimension,
        samplingDistribution = shiftdim(samplingDistribution, 1);       
        distributionDimension = mod(distributionDimension - 2, nDims) + 1;
    end
    
    [~, sample] = max(bsxfun(@minus, cumsum(samplingDistribution), randVec) > 0, [], distributionDimension);
end