#Rory Clune
#NMM Problem 7.2 (c)
from numpy import*

x = array([1.0,0.01])#seed
y = empty([2])
M = zeros([3,2])#moments
C = zeros([3,2])#cumulants

iter = 10000; n = 0
while (n<iter):
    n+=1    
    #Generate the next set of random numbers (linear congruential), between 1 and 134456
	#(division done later to bring x into the range 0 < x < =1 )
    x = ((8121*x + 28411)%134456)
    #Apply transform
    y[0] = sqrt(-2*log(x[0]/134456))*sin(x[1])
    y[1] = sqrt(-2*log(x[0]/134456))*cos(x[1])
    M[0,:] += y #Expected value of y
    M[1,:] += y**2 #Expected value of y^2
    M[2,:] += y**3 #Expected value of y^3   
M/=iter	
C[0] = M[0,:]
C[1] = M[1,:]-M[0,:]**2
C[2] = M[2,:]-2*M[0,:]**3-3*M[0,:]*M[1,:]
print "Cumulants: "; print C