Gianluca Dorini

I've tried this small experiment:

N = 100000; % sample size
w = 1:100; % quite steep distribution
w = w / sum(w); % distribution (pdf)
Y_std = sqrt(w .* (1 - w)); % Standard deviations associated to each probability

%%%%%%%%%%%%%%%%%

Y1 = randp(w, N, 1); % Try Jos function
Y1_w = histc(Y1, [1:100]); % absolute frequencies table
Y1_w = Y1_w / N; % transform into relative frequencies
Y1_w_std = sqrt(Y1_w .* (1 - Y1_w)); % Empirical standard deviations associated to each probability

%%%%%%%%%%%%%%%%%

Y2 = gDiscrPdfRnd(w, N, 1); % Try Gianluca Dorini's function
Y2_w = histc(Y2, [1:100]); % absolute frequencies table
Y2_w = Y2_w / N; % transform into relative frequencies
Y2_w_std = sqrt(Y2_w .* (1 - Y2_w)); % Empirical standard deviations associated to each probability

%%%%%%%%%%%%%%%%%

% Check if empirical standard deviations are coherent with theretical ones
plot(Y_std - Y1_w_std', '.b') % plot empirical differences for Jos function
hold on
plot(Y_std - Y2_w_std', '.r') % plot empirical differences for Gianluca Dorini's function

%%%%%%%%%%%%%%%%%

It looks as though both functions statistically behave very well. They are both excellent functions. Thank you.

Kinda slow for large dimensions due to the for loop.

Input and output are actually a matrix of column vectors, not row vectors as stated in the documentation. (matrix dimensions are usually row x column, not column x row)

Faster version ** (note this doesn't produce the exact same output as the original HyperSphere.m function, see below):
function Y = HyperSphere2(X,r=1)
% compute the cartesian coordinates from the input
% hypersphere angles and radius

% n = number of column vectors
n = size(X,2);

% cosine terms
C = [cos(X); ones(1,n)];
% sine terms
S = [ones(1,n); cumprod(sin(X),1)];
% calculate output
Y = r * C .* S;

=========
This version produces the same output as the original HyperSphere
.m (swapping the sines and cosines, and change the order of the vectors)
function Y = HyperSphere2(X,r=1)
% compute the cartesian coordinates from the input
% hypersphere angles and radius

% dimension of the input vectors
m = size(X,1);
% number of input vectors
n = size(X,2);

% reverse the order of X
X = X(m:-1:1,:);

% sine terms
S = [sin(X); ones(1,n)];
% cosine terms
C = [ones(1,n); cumprod(cos(X),1)];
% calculate output
Y = r * C .* S;

% reverse the order of Y
Y = Y(m+1:-1:1,:);