RANDP

Random integers with given probabilities ( (v2.0)

Description

  RANDP - pick random values with relative probability
 
      R = RANDP(PROB,..) returns integers in the range from 1 to
      NUMEL(PROB) with a relative probability, so that the value X is
      present approximately (PROB(X)./sum(PROB)) times in the matrix R.
 
      All values of PROB should be equal to or larger than 0.
 
      RANDP(PROB,N) is an N-by-N matrix, RANDP(PROB,M,N) and
      RANDP(PROB,[M,N]) are M-by-N matrices. RANDP(PROB, M1,M2,M3,...) or
      RANDP(PROB,[M1,M2,M3,...]) generate random arrays.
      RANDP(PROB,SIZE(A)) is the same size as A.
 
      Examples:
        R = randp([1 3 2],1,10000) ;
        % return a row vector with 10000 values with about 16650 two's
        histc(R,1:3) ./ numel(R)
 
        R = randp([1 1 0 0 1],10,1)
        % 10 samples evenly drawn from [1 2 5]
 
        % randomly select 100 elements according to a specific distribution
        V = {'Red','Green','Blue'} ;
        ind = randp([80 10 10],200,1) ;
        R = V(ind) ; % should contain about 160 'Red'
 
      Also see RAND, RANDPERM
               RANDPERMBREAK, RANDINTERVAL, RANDSWAP (MatLab File Exchange)

version 2.0, feb 2009

• MATLAB release: MATLAB 6.5 (R13)

[urs (us) schwarz]

nicely coded snippet
two (unimportant) comments:
1) why not put the validity test for P BEFORE the (potentially consuming) generation of all the random numbers...
2) why not add a pointer to RANDSAMPLE for those who own the stats tbx...
us

[Reza Farrahi Moghaddam]

Very useful when dealing with user-defined PDFs.

[Carlos Baiz]

Works great!

[Francesco Pozzi]

Wow, it really is a nice function, very fast and efficient. It is exactly what I was looking for. Thank you.

[Francesco Pozzi]

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.

[Martin]

What an elegant piece of code! Thanks for this.

[Reza Fathzadeh]