data =  data - data.mean(axis=0)

cov = np.dot(data.T, data) / data.shape[0] # get covariance matrix
eigs, eigv = np.linalg.eigh(cov)
normed_eigs = eigs[::-1] / np.sum(eigs)
eigs_sum = np.cumsum(normed_eigs)
n_components = np.argmax(eigs_sum > 0.90) # get components with total variance of 0.9
P = eigv * np.sqrt(1.0/(eigs + 0.01))
P = P[:, -n_components:]
transformed = np.dot(data, P) # get pca-whitened data

eigv = eigv[:, -n_components:]
np.dot(transformed,np.sqrt(eigs[:, np.newaxis]) * eigv) + mean