Problem 13.3

In [118]:
from __future__ import division
import numpy as np
import random
import matplotlib.pyplot as plt

C = [[1, 0, 1], [0, 1, 1], [1, 1, 2]]
w,v = np.linalg.eig(C)
print "eigenvals: "
print w
print ""
print "eigenvects:"
print v

eigenvals: 
[  3.00000000e+00   1.00000000e+00  -3.36770206e-17]

eigenvects:
[[ -4.08248290e-01  -7.07106781e-01  -5.77350269e-01]
 [ -4.08248290e-01   7.07106781e-01  -5.77350269e-01]
 [ -8.16496581e-01  -2.61214948e-16   5.77350269e-01]]
In [125]:
X = []
numPts = 1000;

for i in range(numPts):
    x0 = random.gauss(0, 1)
    x1 = random.gauss(0, 1)
    X.append([x0, x1, x0+x1])

X = np.array(X)
C = np.array([
    [np.dot(X[:,0], X[:,0]),np.dot(X[:,0], X[:,1]), np.dot(X[:,0], X[:,2])], 
    [np.dot(X[:,1], X[:,0]),np.dot(X[:,1], X[:,1]), np.dot(X[:,1], X[:,2])], 
    [np.dot(X[:,2], X[:,0]),np.dot(X[:,2], X[:,1]), np.dot(X[:,2], X[:,2])]
])/numPts
    
print "covar:"
print C
print ""
w,v = np.linalg.eig(C)
print "eigenvals: "
print w
print ""
print "eigenvects:"
print v

covar:
[[ 1.03719509 -0.01237133  1.02482376]
 [-0.01237133  0.98232752  0.96995619]
 [ 1.02482376  0.96995619  1.99477995]]

eigenvals: 
[  2.99331535e+00   1.02098721e+00   4.08910584e-16]

eigenvects:
[[-0.42517007  0.69706318  0.57735027]
 [-0.39108938 -0.71673967  0.57735027]
 [-0.81625946 -0.0196765  -0.57735027]]
In [138]:
M = (v.T)
M = M[:,:2]#remove zero eigenvector
y = np.dot(X, M)#diagonalize

C2 = np.array([
    [np.dot(y[:,0], y[:,0]),np.dot(y[:,0], y[:,1])], 
    [np.dot(y[:,1], y[:,0]),np.dot(y[:,1], y[:,1])]
])/numPts

print "covar:"
print C2
print ""
w2,v2 = np.linalg.eig(C2)
print "eigenvals: "
print w2
print ""
print "eigenvects:"
print v2

covar:
[[ 1.61464883 -0.14777508]
 [-0.14777508  0.05571196]]

eigenvals: 
[ 1.6285331   0.04182769]

eigenvects:
[[ 0.9956152   0.09354345]
 [-0.09354345  0.9956152 ]]

Problem 13.4

In [139]:
from __future__ import division
import numpy as np
import random
import matplotlib.pyplot as plt

numPts = 5000;

ptsX = []
ptsY = []

for i in range(numPts):
    ptsX.append(random.random())
    ptsY.append(random.random())

plt.scatter(ptsX, ptsY, marker=".")
plt.axes().set_aspect('equal', 'datalim')
plt.show()
In [140]:
transform = [[1,2], [3,1]]

for i in range(numPts):
    x = ptsX[i]
    y = ptsY[i]
    ptsX[i]  = transform[0][0]*x + transform[0][1]*y
    ptsY[i]  = transform[1][0]*x + transform[1][1]*y
    
plt.scatter(ptsX, ptsY, marker=".")
plt.axes().set_aspect('equal', 'datalim')
plt.show()
In [144]:
#nice explanation of PCA:
#https://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues

#make zero mean
npPtsX = np.array(ptsX)
npPtsY = np.array(ptsY)
npPtsX -= np.mean(npPtsX)
npPtsY -= np.mean(npPtsY)

#calc covariance
C = np.array([
    [np.dot(npPtsX, npPtsX),np.dot(npPtsX, npPtsY)], 
    [np.dot(npPtsX, npPtsY), np.dot(npPtsY,npPtsY)]]
)/numPts

w,v = np.linalg.eig(C)

M = (v/np.sqrt(w)).T#divide by sqrt(w) for unit variance
x = np.array([npPtsX, npPtsY])
y = np.dot(M, x)#diagonalize

C = np.array([
    [np.dot(y[0], y[0]),np.dot(y[0], y[1])], 
    [np.dot(y[0], y[1]), np.dot(y[1], y[1])]]
)/numPts

print "covariance:"
print C

plt.scatter(y[0], y[1], marker=".")
plt.axes().set_aspect('equal', 'datalim')
plt.show()

covariance:
[[  1.00000000e+00   1.50635060e-16]
 [  1.50635060e-16   1.00000000e+00]]
In [157]:
# Find the independent components of x with the log cosh contrast function, and plot.
#f(x) = log(cosh(x))
#from wolfram alpha:
#df = tanh(_x)
#ddf = (1/cosh(x))**2

def df(_x):
    return np.tanh(_x)
def ddf(_x):
    return (1/np.cosh(_x))**2

def expectation(data):
    return np.array([np.mean(data[0]), np.mean(data[1])])

num = 0;
def stepICA(w):
    global num
    w_next = expectation(y*df(np.dot(w, y))) - w*expectation(ddf(np.dot(w, y)))#13.37
    w_next /= np.linalg.norm(w_next)#13.38 - sq is wrong according to https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf
    num+=1;
    if (np.linalg.norm(w_next+w)<1e-8 or np.linalg.norm(w_next-w)<1e-8 or num > 100): 
        print str(num) + " iterations"
        return w_next
    return stepICA(w_next)

W = stepICA(np.array([random.random(), random.random()]))

print ""
print "w:"
print W
plt.scatter(y[0], y[1], marker=".")
plt.plot([0, W[0]], [0, W[1]], 'r')
plt.plot([0, -W[1]], [0, W[0]], 'r')
plt.axes().set_aspect('equal', 'datalim')
plt.show()

12 iterations

w:
[ 0.55784119 -0.82994771]