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
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
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
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()
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()
#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()
# 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()