import time
import numpy

def uniform_random_deviate( size,  seed=None ):
    """
    Uses multiplicative congruential method.
    Return size uniform pseudorandom deviates in [0,1) as a numpy array
    Also see numpy.random.uniform
    mean = 0.5, variance = 1/12
    """
    if seed == None:
        seed = int((time.time() % 1) * 2**32)
    a = 16807
    M = 2**32-1
    x = numpy.zeros(size)
    x[0] = float((a*seed) % M)
    for i in xrange(1,size):
        x[i] = float((a*x[i-1]) % M)
    x /= M
    return x

def gaussian_random_deviate( size, N = 12 ):
    """
    Approximates Gaussian deviates by the sum of N 
    uniform random deviates
    Return size deviates of Standard (mean,var=1) Gaussian Deviation
    Also see numpy.random.normal
    """
    x = numpy.zeros(size)
    for i in xrange(size):
        x[i] = numpy.sum(uniform_random_deviate(N))
    x -= float(N)/2. # mean is zero
    x /= numpy.sqrt(N/12)  # variance is 1
    return x


if __name__ == "__main__":
    import pylab
    g = gaussian_random_deviate(100)
    print 'mean =', g.mean()
    print 'variance = ', g.var()
    hist ,be, histgraph = pylab.hist(g,20, normed=True, color='gray', ec = 'white')
    x = numpy.arange(-4,4,.1)
    def gaussian(x) : 
        return 1/numpy.sqrt(2*numpy.pi) * numpy.exp(-x**2/2)
    pylab.plot(x,gaussian(x),color='black')
    pylab.xlabel('x')
    pylab.ylabel('p(x)')
    pylab.legend(('Standard Gaussian','Deviates'),loc='upper left')
    pylab.savefig('gaussian_deviate.png' )