## Tehtävä 2. robusti<-function(x) { sigma1<-cov(x) ## ensimmäinen hajontamatriisi S1<-solve(mat.sqrt(sigma1)) xcent<-sweep(x,2,apply(x,2,mean),"-") z<-xcent%*%t(S1) sigma2<-cov4_t2(z) sigma2 lambda<-diag(eigen(sigma2)$values) ## ominaisarvot gamma<-eigen(sigma2)$vectors ## ominaisvektorit gamma%*%lambda%*%t(gamma) ## sama kuin sigma2 W<-t(gamma)%*%solve(mat.sqrt(cov(x))) W } # laskee huipukkuusmatriisin cov4_t2<-function(x) { d<-dim(x)[2] n<-dim(x)[1] xm<-sweep(x,2,colMeans(x),"-") z<-xm%*%solve(mat.sqrt(cov(x))) r<-sqrt(diag(z%*%t(z))) rx<-sweep(xm,1,r,"/") ## tässä muutos (t(rx)%*%rx/n)*d ## ja tässä } ## Tehtävä 3. # Tavallinen FOBI fobi<-function(x) { S1<-solve(mat.sqrt(cov(x))) xcent<-sweep(x,2,apply(x,2,mean),"-") S2<-cov4(xcent%*%t(S1)) U2<-eigen(S2)$vectors t(U2)%*%S1 } # laskee matriisin neliöjuuren mat.sqrt <- function(V) { apu <- eigen(V) L <- apu$values P <- apu$vectors P%*%(diag(L))^(1/2)%*%t(P) } # laskee huipukkuusmatriisin cov4<-function(x) { d<-dim(x)[2] n<-dim(x)[1] xm<-sweep(x,2,colMeans(x),"-") z<-xm%*%solve(mat.sqrt(cov(x))) r<-sqrt(diag(z%*%t(z))) rx<-sweep(xm,1,r,"*") (t(rx)%*%rx/n)/(d+2) } ## a) data<-read.table("http://users.jyu.fi/~slahola/files/ica_datoja/speech.dat", header=T) data<-as.matrix(data) pairs(data) A<-matrix(runif(16), ncol=4) x<-data%*%t(A) W1<-fobi(x) W1 Xcent<-sweep(x,2,apply(x,2,mean),"-") SestF<-Xcent%*%t(W1) pairs(SestF) x11() W2<-T2(x) W2 Sest2<-Xcent%*%t(W2) pairs(Sest2) ## b) ## Alkuperäinen par(mfrow=c(dim(data)[2],1)) for(i in 1:dim(data)[2]) { plot(data[,i], type="l") } x11() ## Fobi par(mfrow=c(dim(SestF)[2],1)) for(i in 1:dim(SestF)[2]) { plot(SestF[,i], type="l") } x11() ## T2 par(mfrow=c(dim(Sest2)[2],1)) for(i in 1:dim(Sest2)[2]) { plot(Sest2[,i], type="l") } ISR<-function(W,A) { G<-W%*%A delta<-abs(sweep(G,1,apply(abs(G),1,max),"/")) sqrt((1/dim(G)[2])*sum(delta^2)-1) } ISR(W1,A) # 0.4961508 ISR(W2,A) # 0.6050086 ## -> Fobi parempi ################################################################## # ICA-DEMO 6 TEHTÄVÄ 4 a) #FOBI # mat.sqrt -funktio laskee matriisin neliöjuuren library(fastICA) ISR<-function(W,A){ G<-W%*%A delta<-abs(sweep(G,1,apply(abs(G),1,max),"/")) sqrt((1/dim(G)[2])*sum(delta^2)-1) } simul<-function(n,A,rep) { ISR1<-NULL ISR2<-NULL for(i in 1:rep){ d<-dim(x)[2] s<-cbind(rnorm(n),rexp(n),runif(n,-sqrt(3),sqrt(3))) x<-s%*%t(A) W1<-fobi(x) W2<-robusti(x) ISR1<-c(ISR1,(ISR(W1,A))) ISR2<-c(ISR2,(ISR(W2,A))) } list(ISR1=ISR1,ISR2=ISR2) } A<-matrix(runif(9),3,3) result<-simul(500,A,100) x11() boxplot(data.frame(result)) ################################################################## # ICA-DEMO 6 TEHTÄVÄ 4 b) library(fastICA) ISR<-function(W,A){ G<-W%*%A delta<-abs(sweep(G,1,apply(abs(G),1,max),"/")) sqrt((1/dim(G)[2])*sum(delta^2)-1) } simul<-function(n,A,rep) { ISR1<-NULL ISR2<-NULL for(i in 1:rep){ d<-dim(x)[2] s<-cbind(rt(n,df=1),rt(n,df=3),rt(n,df=5)) x<-s%*%t(A) W1<-fobi(x) W2<-robusti(x) ISR1<-c(ISR1,(ISR(W1,A))) ISR2<-c(ISR2,(ISR(W2,A))) } list(ISR1=ISR1,ISR2=ISR2) } A<-matrix(runif(9),3,3) result<-simul(500,A,500) x11() boxplot(data.frame(result)) ######################### sim <- function(n,A,rep){ T1 <- NULL T2 <- NULL d <- dim(A)[1] for(i in 1:rep){ s <- matrix(runif(n*d),n,d) x <- s %*% t(A) ptm <- proc.time() fastICA(x,d,alg.typ="parallel",fun="logcosh",method="C") T1[i] <- (proc.time()-ptm)[3] ptm <- proc.time() fobi(x) T2[i] <- (proc.time()-ptm)[3] } T1 <- NULL T2 <- NULL for(i in 1:length(d)){ A <- matrix(runif(d[i]^2),d[i],d[i]) result <- sim(n,A,10) T1[i] <- median(result$fastICA) T2[i] <- median(result$fobi) } n <- c(50,100,200,500,800) d <- 30 T1 <- NULL T2 <- NULL for(i in 1:length(n)){ A <- matrix(runif(d^2),d,d) result <- sim(n[i],A,20) T1[i] <- median(result$fastICA) T2[i] <- median(result$fobi) } MEG <-read.table("http://users.jyu.fi/~slahola/files/ica_datoja/MEG.dat") MEG <- as.matrix(MEG) #alkuperäiset par(mfrow = c(10,1), mar=c(1,1,1,1)) for(i in 31:40){ plot(MEG[,i], type = "l", main="alkuperäiset") } # tehtävä 7: #matrix based on fourth moments cov2 <- function(x){ n <- dim(x)[1] d <- dim(x)[2] m <- apply(x,2,mean) xcent <- sweep(x,2,m,"-") S1 <- cov(x) z <- xcent %*% solve(mat.sqrt(S1)) r <- NULL for(i in 1:n){ r[i] <- t(z[i,]) %*% z[i,] } rx <- sweep(xcent,1,sqrt(r),"*") (t(rx) %*% rx)/(n*(d+2)) } fobi <- function(x){ m <- apply(x,2,mean) xcent <- sweep(x,2,m,"-") S1 <- cov(x) z <- xcent %*% solve(mat.sqrt(S1))#valkaistut havainnot S2 <- cov2(z) V <- eigen(S2)$vectors W <- t(V) %*% solve(mat.sqrt(S1)) S <- x %*% t(W) list(S = S, W = W) } cov2r <- function(x){ n <- dim(x)[1] d <- dim(x)[2] m <- apply(x,2,mean) xcent <- sweep(x,2,m,"-") S1 <- cov(x) z <- xcent %*% solve(mat.sqrt(S1)) r <- NULL for(i in 1:n){ r[i] <- t(z[i,]) %*% z[i,] } rx <- sweep(xcent,1,sqrt(r),"/") d*(t(rx) %*% rx)/n } fobiM <- function(x){ m <- apply(x,2,mean) xcent <- sweep(x,2,m,"-") S1 <- cov(x) z <- xcent %*% solve(mat.sqrt(S1))#valkaistut havainnot S2 <- cov2r(z) V <- eigen(S2)$vectors W <- t(V) %*% solve(mat.sqrt(S1)) S <- x %*% t(W) list(S = S, W = W) } fob<-fobi(MEG) fobM <- fobiM(MEG) fica <- fastICA(MEG,40,alg.typ="parallel", fun="exp") x11() #fastICA par(mfrow = c(10,1), mar=c(1,1,1,1)) for(i in 1:10){ plot(fica$S[,i], type = "l", main="fastICA") } x11() #FOBI par(mfrow = c(10,1), mar=c(1,1,1,1)) for(i in 1:10){ plot(fob$S[,i], type = "l", main="FOBI") } x11() #robusti par(mfrow = c(10,1), mar=c(1,1,1,1)) for(i in 1:10){ plot(fobM$S[,i], type = "l", main="robusti") }