Derek O'Connor

Top Class, especially with newly-compiled mexw64 files for Windows 7 64-bit

Here is a function that beats those above by a mile:

function S = DiscSampVec2(x,p,ns);
%
[~,idx] = histc(rand(1,ns),[0,cumsum(p)]);
S = x(idx);

This is a slight modification of the function on page 47, Kroese, Taimre, and Botev, Handbook of Monte Carlo Methods, Wiley, 2011.

>> ns=10^6; n=10^3; x=1:n;
>> p = rand(1,n); p = p/sum(p);

>> t2=tic;
>> S = DiscSampVec2(x,p,ns);
>> t2=toc(t2)

      0.15835

It would be useful if the author would explain clearly what P, M, N are. For example, if P is a "discrete probability distribution for the indices of P", then why does it need to be normalized? Also I think the author means "density" or "mass function" rather than "distribution".

Is M the size of the sample and is N the length of P?

I presume the author is trying to do what this simple function does:

function S = DiscITSamp1(x,p,ns);

% Generate a random sample S of size ns
% with replacement from x(1:n) with prob. % density p(1:n), using the discrete
% inverse transform method.
% Time Complexity: O(n*ns).
% Derek O'Connor, 31 July 2011.
% derekroconnor@eircom.net
%
% USE: S = DiscITSamp1([5 7 8 11],
% [0.2 0.3 0.4 0.1],1000);hist(S)

cdf = cumsum(p);
S = zeros(1,ns);
for k = 1:ns
    u = rand;
    i = 1;
    while cdf(i) <= u
        i = i+1;
    end;
    S(k) = x(i);
end

This is three times faster than GENDIST on a 2.3GHz machine, Matlab2008a 64-bit

>> ns=10^6;n = 10^3;x=1:n;p=rand(1,n);p=p/sum(p);
>> tic;S1 = DiscITSamp1(x,p,ns);t=toc;disp([ns t])
       1e+006 5.6115

>> tic;S2 = gendist(p,1,ns);t=toc;disp([ns t])
       1e+006 19.278

Jan,

(1). No, your modification still uses 2 x p memory. So the "clear p" has to stay in p=GRPfys(n) to conserve memory.

(2). Is there a typo in "p=zeros(1, n); p(1) = 1;" ?

p = GRPfysSimon(5) gives [0 0 1 0 0], and moves the '1' randomly in subsequent calls.

All the generators I've written must start with a (any?) permutation p and then do random transpositions on p,
p = Trans(p,i,r), using the three-statement swap: t = p(i); p(i) = p(r); p(r) = t.

The correctness proofs of the algorithms assume this condition.

I do not use Knuth's two-statement short-cut for permutations of [1,2,...,n] and so GRPfys(p) works for any p. For example, if we want to generalize GRPfys to get a permutation of the integers [L,L+1,...,U] then the only change needed in GRPfys is replace
"p=I(n)" with "p = L:U; n = U-L+1" to get p = GRPfysLU(L,U). This works for
p = GRPfysLU(-5,5) and p = GRPfysLU(-6,-2).

Alternatively, "p = GRPfys(U-L+1)+L-1" does the same thing but is error-prone and puts a burden on the user.

(3). I'll put a link to Shuffle in the next update.

(4). GRPfysSimon(3) neatly simulates the Shell or Three-Card-Trick game.

Thanks for you interest,

Derek.

"p(:)=GRPfys(n);" causes an assignment error:
// ??? In an assignment A(:) = B, the number of elements in A and B must be the same.
//
// Error in ==> GRDrej at 17
// p(:) = GRPfys(n);

Again, this shows that the semantics (meaning) of Matlab statements are not well-defined.

Here is another example of "chaotic" semantics and the "perils of vectorization", which I found this afternoon when writing and testing functions for inverting a permutation p:

function q = InvPCmp(p);
%
% Invert the permutation p so that p(q) = q(p) = I
%
% Derek O'Connor 31 Jan 2011.
%
n = length(p);
q = zeros(1,n);
for i = 1:n
   q(p(i)) = i;
end;
return % InvPCmp(n)
%
% Here are two 'slick' one-line alternatives
%
% (1) InvPSrt: [t,q] = sort(p);
%
% (2) InvPVec: q(p) = 1:length(p);
%

Alternative (1) InvPSrt(p), uses the same trick that Matlab's randperm uses.

Alternative (2) InvPVec(p), is what all good Matlab-ers like: a vectorized version of the component-style function InvPCmp(p).

Here are the last lines of timing-memory test results:

InvPCmp: n = 2^29 tg = 195 ti = 70 ti/tg = 0.36 Mem = 8.0 GB

InvPSrt: n = 2^29 tg = 196 ti = 142 ti/tg = 0.72 Mem = 11.5 GB

InvPVec: n = 2^28 tg = 91 ti = 38 ti/tg = 0.42 Mem = 9.5 GB
         n = 2^29 tg = 197 started page swapping Mem > 16.0 GB

You can see that the component version InvPCmp is the fastest and uses only the memory necessary to store two vectors p and q of length n = 2^29.

The sort version InvPSrt takes twice as long and uses almost 50% more memory than InvPCmp.

The vectorized version InvPVec was a complete disaster: it took longer on the smaller vectors and ran out of memory for n = 2^29.

Derek O'Connor