How to generate t random numbers using Gauss

Here is a fast way to generate student t random matrix with df degrees of freedom:

proc rndt(r,c,df);  /* r:rows, c:cols, df:df */
    local x,z;
    x = rndn(df,r*c);
    x = sumc(x.*x);
    x = reshape(x,r,c);
    z = rndn(r,c);
    retp(z./sqrt(x/df));
endp;
/* Example: y = rndt(100,5,2); */

Note.

  1. This is a procedure (almost) everybody knows. I just present it here for those who don't want to bother programming (like myself).
  2. This procedure is fast because it does not operate element-wise.
  3. You can make the procedure shorter by condensing multiple lines into one.
    The current proc is written longer than necessary in order to enhance readability.
  4. You can easily modify it to use rndns instead of rndn.
  5. The theory behind this procedure is this: Let X be chi-square(df) and Z be N(0,1). Then Z/(X/df)1/2 is distributed as t with df degrees of freedom.
  6. Please report any errors to me by email (c[dot]han[at]auckland.ac.nz). I would greatly appreciate it.