import random import matplotlib.pyplot as plt import math import numpy as np x = np.array([random.random() for i in xrange(1000)]) y = np.array([random.random() for i in xrange(1000)]) vec = np.concatenate(( x[:,np.newaxis], y[:,np.newaxis]),1) transform = np.matrix([[-.2,1.1],[1.1,-.6]]) trans_vec = np.array((transform* np.matrix(vec).transpose()).transpose()) x_mean = sum(trans_vec[:,0])/ len(trans_vec[:,0]) y_mean = sum(trans_vec[:,1])/ len(trans_vec[:,1]) centered = np.empty_like(vec) centered[:,0] = [i-x_mean for i in trans_vec[:,0] ] centered[:,1] = [i-y_mean for i in trans_vec[:,1] ] C = np.cov(centered, rowvar = 1) eigval, eigenvect = np.linalg.eig(C) eigenvect[:,0] /= eigval[0]**(1.0/len(eigval)) eigenvect[:,0] /= eigval[1]**(1.0/len(eigval)) M = np.transpose(eigenvect[:,abs(eigval)>1e-8]) xt = np.dot(M,centered) Cxt = np.cov(xt) # xt should have a diagonal covariance matrix, with no zero eigenvalues # plot transformed data: print xt plt.scatter(vec[:,0],vec[:,1],c='y') plt.scatter(trans_vec[:,0],trans_vec[:,1], c="r") plt.scatter(centered[:,0],centered[:,1]) plt.scatter(xt[:,0],xt[:,1]) plt.axis([-1,2, -1,2]) plt.show()