#
# ica.r
# (c) Neil Gershenfeld  4/10/03
# 2D ICA example
#
npts = 1000
eps = 1e-6
#
# generate data
#
s = rbind(matrix(runif(npts,min=0,max=1),1,npts),
          matrix(runif(npts,min=0,max=1),1,npts))
par(mfrow=c(3,2),mai=c(.2,.1,.2,.1))
plot(s[1,],s[2,],xlab="data",ylab="",axes=FALSE,
   xlim=c(-.1,1.2),ylim=c(-.1,1.2),pch=16,cex=.9)
title("data")
abline(h=0,v=0)
#
# mix data
#
A = matrix(c(1,3,2,1),2,2)
x = A %*% s
plot(x[1,],x[2,],xlab="mixed",ylab="",axes=FALSE,
   xlim=c(-.3,4),ylim=c(-.3,4),pch=16,cex=.9)
title("mixed")
abline(h=0,v=0)
#
# subtract mean
#
x = x - apply(x,1,mean)
plot(x[1,],x[2,],xlab="",ylab="",axes=FALSE,
   xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),pch=16,cex=.9)
title("zero mean")
abline(h=0,v=0)
#
# diagonalize
#
C = cov(t(x))
M = t(eigen(C)$vectors)
x = M %*% x
plot(x[1,],x[2,],xlab="",ylab="",axes=FALSE,
   xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),pch=16,cex=.9)
title("diagonalize")
abline(h=0,v=0)
#
# normalize
#
v = eigen(C)$values
x = diag(1/sqrt(v)) %*% x
plot(x[1,],x[2,],xlab="",ylab="",axes=FALSE,
   xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),pch=16,cex=.9)
title("normalize")
abline(h=0,v=0)
#
# ICA iteration for first component
#
f = function(x) (log(cosh(x)))
df = function(x) (tanh(x))
ddf = function(x) (1-tanh(x)^2)
w1 = matrix(runif(2,min=0,max=1),2,1)
w1 = w1 / sqrt(sum(w1 * w1))
err = 1
i = 0
while (err > eps) {
   i = i+1
   w1new = apply(x * t(array(df(t(w1) %*% x),dim=c(npts,2))),1,mean)
      - mean(ddf(t(w1) %*% x)) * w1
   w1new = w1new / sqrt(sum(w1new * w1new))
   err = sum(abs(w1-w1new))
   w1 = w1new
   }
cat(i,"iterations\n")
#
# since we're in 2D, orthogonalize to find other vector
#
w2 = c(w1[2],-w1[1])
W = rbind(w1,w2)
y = W %*% x
plot(y[1,],y[2,],xlab="",ylab="",axes=FALSE,
   xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),pch=16,cex=.9)
title("ICA")
abline(h=0,v=0)
dev.print(postscript,file="ica.ps")