StatRed: Implemented Lab Excercises 21 and 22.

statred/ass1/ex21.py

-#from math import mean
-#import matplotlib.pyplot as plt
-from numpy.random import multivariate_normal
-
-def mu_est(x):
-    return sum(x) / float(len(x))
-
-def cov_est(x):
-    avg = mean(x)
-    n = len(x)
-    # TODO: transpose the second "(i - avg)"
-    return sum([(i - avg)*(i - avg) for i in x]) / float(n-1)
-
-mean = [0,0]
-cov  = [[1,0],[0,100]] # diagonal covariance, points lie on x or y-axis
-N    = 512 # Samples to use.
-
-x,y = multivariate_normal(mean, cov, N).T
-print 'x:', x
-print 'y:', y
-
-print 'mu_est(x):', mu_est(x)
-print 'mu_est(y):', mu_est(y)

statred/ass1/q21-multivariate.py → statred/ass1/q21_multivariate.py

 from pylab import array, eig, diagflat, dot, sqrt, randn, tile, \
-        plot, subplot, axis, show, figure, clf
+        plot, subplot, axis, figure, clf, savefig
 
-mean = array([[3],
+mu = array([[3],
               [4],
               [5],
               [6]])
@@ -15,20 +15,21 @@ cov  = array(
 samples = 1000
 vector_size = 4
 
-figure(16)
-clf()
+def dataset():
+    d, U = eig(cov)
+    L = diagflat(d)
+    A = dot(U, sqrt(L))
+    X = randn(vector_size, samples)
+    return dot(A,X) + tile(mu, samples)
 
-d, U = eig(cov)
-L = diagflat(d)
-A = dot(U, sqrt(L))
-X = randn(vector_size, samples)
-Y = dot(A,X) + tile(mean, samples)
-
-for i in range(1, 5):
-    for j in range(1, 5):
-        if i != j:
-            subplot(4, 4, i+(j-1)*4)
-            plot(Y[i-1], Y[j-1], 'x')
-            axis('equal')
-
-show()
+if __name__ == '__main__':
+    figure(vector_size**2)
+    clf()
+    Y = dataset()
+    for i in range(vector_size):
+        for j in range(vector_size):
+            if i != j:
+                subplot(vector_size, vector_size, (i+1) + j*vector_size)
+                plot(Y[i], Y[j], 'x')
+                axis('equal')
+    savefig('figures/q21.pdf')

statred/ass1/q22_estimate.py

+from q21_multivariate import dataset, samples
+from numpy import array, mean, tile, newaxis, dot
+
+n = 1000
+
+Y = array([mean(dataset(), 1) for i in range(n)]).T
+mu = mean(Y, 1)
+Yzm = Y - tile(mu[:,newaxis], n)
+S = dot(Yzm, Yzm.T) / (n - 1)
+
+print 'S:', S