from numpy import*
from numpy import linalg as LA

n = 10000
x = zeros((n,3))

u = array([1.0,0.01])#  seed for uniformly distributed pseudorandom generation
for i in range (0,n):
    #  Generate pseudorandom numbers from a Gaussian distro with mean 0 and var 1    
    u = ((8121*u + 28411)%134456)        
    x[i,0] = sqrt(-2*log(u[0]/134456))*sin(u[1])#transform uniform to Gaussian
    x[i,1] = sqrt(-2*log(u[0]/134456))*cos(u[1])
    x[i,2] = x[i,0]+x[i,1]
Cx = cov(x, rowvar = 0)
print '\nCovariance matrix of x:\n' + str(Cx)

eig, phi = LA.eig(Cx)
print '\neigenvalues of the covariance matrix of x:\n' + str(eig)

#  M contains the eigenvectors of Cx (corresponding to nonzero eigenvalues) as columns
M = transpose(phi[:,fabs(eig)>1e-8])
#  Transfrom the data points x
y = dot(M,transpose(x))
Cy = cov(y)
#  y should have a diagonal covariance matrix, with no zero eigenvalues
print '\nCovariance matrix of y:\n' + str(Cy)