StatRed: completed question 21.

statred/ass1/q21-multivariate.py

+from pylab import array, eig, diagflat, dot, sqrt, randn, tile
+from pylab import plot, axis, show, figure, clf
+
+mean = array([[3],
+              [4],
+              [5],
+              [6]])
+
+cov  = array(
+    [[ 3.01602775,  1.02746769, -3.60224613, -2.08792829],
+     [ 1.02746769,  5.65146472, -3.98616664,  0.48723704],
+     [-3.60224613, -3.98616664, 13.04508284, -1.59255406],
+     [-2.08792829,  0.48723704, -1.59255406,  8.28742469]])
+
+samples = 1000
+vector_size = 4
+
+figure(1)
+clf()
+
+d, U = eig(cov)
+L = diagflat(d)
+A = dot(U, sqrt(L))
+X = randn(vector_size, samples)
+Y = dot(A,X) + tile(mean, samples)
+
+#import matplotlib.pyplot as plt
+for i in range(4):
+    for j in range(4):
+        if i != j:
+            plot(Y[i], Y[j], 'x')
+
+axis('equal')
+#plt.plot(x,y,'x');
+#plt.axis('equal');
+#plt.show()
+show()
+
+#import numpy as np
+#x,y = np.random.multivariate_normal(mean, cov, 5000).T
+#print 'x:', x
+#print 'y:', y
+#plt.plot(x,y,'x'); plt.axis('equal'); plt.show()
+