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