In [60]:

    

from numpy.random import randn

rescale = lambda x: np.asarray(x) / np.sqrt(np.dot(x, x))
N = 1000
sigma_1 = 1.
sigma_2 = 1.

v1 = rescale([1, 1])
v2 = rescale([1, -1])

xs = [sigma_1 * randn() * v1 + sigma_2 * randn() * v2 
     for i in xrange(N)]


    

In [61]:

    

X = 1. / N * np.sum([x[:, None] * x[None, :] for x in xs], axis=0)
evals, evecs = np.linalg.eigh(X)
sort = np.argsort(evals)
evals, evecs = evals[sort], evecs[sort]
print(evals)
print(evecs)


    

    




[ 0.95651022  1.06432058]
[[-0.72579879 -0.68790705]
 [-0.68790705  0.72579879]]

In [62]:

    

def plot_vec(x, **kwargs):
    pl.plot([0, x[0]], [0, x[1]], **kwargs)
    
pl.figure(figsize=(8,8))
pl.scatter([x[0] for x in xs], [x[1] for x in xs])

plot_vec(evecs[1], color='r')
plot_vec(evecs[0], color='b')

plot_vec(v1, color='r', ls=':')
plot_vec(v2, color='b', ls=':')


    

In [ ]: