import numpy as np
import matplotlib.pyplot as plt

Writeup

Covariance of Normal Distributions and Sums of Normal Distributions

x1 = np.random.normal(size=1000000)
x2 = np.random.normal(size=1000000)
x3 = x1+x2
x = np.array([x1, x2, x3])

covar = np.zeros((3,3))

covar_est = np.array([[1,0,1],[0,1,1],[1,1,2]])

for i in range(3):
    for j in range(3):
        covar[i,j] = np.mean((np.mean(x[i])-x[i])*(np.mean(x[j])-x[j]))

print 'analytic covar:\n',covar_est

print '\ncomputed covar:\n',covar

print '\ndifference:\n',covar-covar_est

analytic covar:
[[1 0 1]
 [0 1 1]
 [1 1 2]]

computed covar:
[[  1.00070778e+00   7.26979685e-04   1.00143476e+00]
 [  7.26979685e-04   1.00091655e+00   1.00164353e+00]
 [  1.00143476e+00   1.00164353e+00   2.00307828e+00]]

difference:
[[ 0.00070778  0.00072698  0.00143476]
 [ 0.00072698  0.00091655  0.00164353]
 [ 0.00143476  0.00164353  0.00307828]]

print 'analytic eigenvalues:\n',np.linalg.eigvals(covar_est)

print '\ncomputed eigenvalues:\n',np.linalg.eigvals(covar)

print '\ndifference:\n',np.linalg.eigvals(covar)-np.linalg.eigvals(covar_est)

analytic eigenvalues:
[  3.00000000e+00   1.00000000e+00  -3.36770206e-17]

computed eigenvalues:
[  3.00461744e+00   1.00008517e+00   8.16726535e-17]

difference:
[  4.61743963e-03   8.51654295e-05   1.15349674e-16]

M = np.linalg.eig(covar_est)[1].T
print "M:\n", M
print "\neigeinvalues on diagonals from M transform:\n", M.dot(covar_est.dot(M.T))

M:
[[ -4.08248290e-01  -4.08248290e-01  -8.16496581e-01]
 [ -7.07106781e-01   7.07106781e-01  -2.61214948e-16]
 [ -5.77350269e-01  -5.77350269e-01   5.77350269e-01]]

eigeinvalues on diagonals from M transform:
[[  3.00000000e+00  -4.09578822e-17  -3.17272656e-16]
 [  7.10771519e-17   1.00000000e+00   7.85046229e-17]
 [ -1.90178410e-16   2.98842800e-17  -6.40987562e-17]]

y = M.dot(x)
y.shape

covar_y = np.zeros((3,3))

for i in range(3):
    for j in range(3):
        covar_y[i,j] = np.mean((np.mean(y[i])-y[i])*(np.mean(y[j])-y[j]))

print '\ncomputed covar y:\n', covar_y

computed covar y:
[[  3.00461742e+00  -1.80800566e-04  -4.08704094e-16]
 [ -1.80800566e-04   1.00008518e+00   7.85424680e-17]
 [ -4.08704094e-16   7.85424680e-17   6.57196212e-32]]

values of the new covariance are virtually zero on the off-diagonals and the diagonals correspond to the eigenvalues of x

Uniform Random Variable transforms

s1 = np.random.uniform(size=1000)
s2 = np.random.uniform(size=1000)

plt.figure()
plt.plot(s1, s2, 'x')
plt.show()

A = np.array([[1,2],[3,1]])

print "transform point by A="
print A

X = A.dot(np.array([s1,s2]))

plt.figure()
plt.plot(X[0],X[1], 'x')
plt.show()

transform point by A=
[[1 2]
 [3 1]]