Build your own random number generator

From Computational Statistics Course Wiki
Jump to navigation Jump to search

The class activity is to make up your own random number generator. The goal is to learn how hard this actually is! If you Google and copy an existing generator from the web, you'll get zero credit unless you can give a detailed explanation of the mathematics behind it. So, I think you're better off to actually try to invent your own. Even though it's unlikely to be very good, you can still win if it is even tolerable (or better than the others that people invent).

Your generator should return i.i.d. uniform random deviates X in the range 0 < X < 1. You should test it by generating a vector of one million deviates and running it through one of the test routines below.

Here are the test routines that you should use to test your generators:

PYTHON:

 import numpy as np
 import scipy.stats
 
 def test_generator(x):
     ''' Input:      x - an array of purported independent Uniform[0, 1]
                     deviates
         Returns:    a mysterious p-value.
     '''
     ntable = 100
     counts, _, _ = np.histogram2d(x[:-1], x[1:],
                                   bins=ntable,
                                   range=[(0, 1), (0, 1)],
                                  )
     chisq, pval = scipy.stats.chisquare(counts.flatten())
     return pval

(If you are using the IPython class server, you can just say 'from our_random import test' to have access to the function.)


MATLAB:

 function pval = testran(x)
   ntable = 100;
   n = numel(x);
   sub1 = 1+fix(ntable*x(1:n-1));
   sub2 = 1+fix(ntable*x(2:n));
   subs = [sub1 sub2];
   counts = accumarray(subs,1,[ntable ntable]);
   ave = (n-1)/(ntable^2);
   chisq = sum( (counts(:)-ave).^2 / ave);
   pval = 1 - chi2cdf(chisq,ntable^2-1);
 end