genrandsorted(p, n)

%% Generate sorted uniformly distributed random vectors
%
% function x = genrandsorted(p, n)
% Generate n sorted p-vectors of variates from U(0,1), i.e. subject to 
% the constraint that each vector satisfies
%     x_1 <= x_2 <= x_3 <= ... <= x_p
% Implements the very clever algorithm of Bentley and Saxe, 1980.
%
% Parameters:
%   p     = dimension of vectors to generate
%   n     = number of vectors to generate
%
% Returns:
%   x     = [n x p] matrix of sorted random vectors
%
% References:
%
% (1) Generating Sorted Lists of Random Numbers
% Jon Louis Bentley and James B. Saxe
% ACM Transactions on Mathematical Software (TOMS)
% Volume 6, Issue 3, September 1980
%
% (c) Copyright Daniel F. Schmidt and Enes Makalic, 2008

function x = genrandsorted(p, n)

% Generate variates
x = zeros(n, p);
CurMax = ones(n, 1);

for i = p:-1:1
    CurMax = CurMax .* rand(n,1).^(1/i);
    x(:,i) = CurMax;
end