data<-read.table("http://users.jyu.fi/~slahola/files/ica_datoja/canc.txt") data<-as.matrix(data) #DEMO2 TEHTÄVÄ 5B mdata<-apply(data,2,mean) Sdata<-cov(data) D<-diag(sqrt(diag(Sdata))) datacent<-sweep(data,2,mdata,"-") #skaalaus zeta<-datacent%*%solve(D) kuvaa<-function(data,n=569){ v<-seq(1,30,by=1) plot(v,data[1,],type="l",ylim=c(min(data),max(data)),ylab="569 käyrää") for(i in 2:n){ points(v,data[i,],type="l") } } kuvaa(zeta) #DEMO2 TEHTÄVÄ 6A G <- eigen(cov(data))$val cumsum(G[1:10])/sum(G) #DEMO2 TEHTÄVÄ 6B G <- eigen(cov(zeta))$val cumsum(G[1:10])/sum(G) #¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤ # 7. tehtävä # a) plottaus<- function(x, S=cov(x)){ plot(x[,1], x[,2],pch=16, col="purple") om.a<- eigen(S)$values om.v<- eigen(S)$vectors mean1<- mean(x[,1]) mean2<- mean(x[,2]) arrows(mean1, mean2, mean1-sqrt(om.a[1])*om.v[1,1], mean2-sqrt(om.a[1])*om.v[2,1], lwd=2) arrows(mean1, mean2, mean1-sqrt(om.a[2])*om.v[1,2], mean2-sqrt(om.a[2])*om.v[2,2], lwd=2) } # b) x<- read.table("http://users.jyu.fi/~slahola/files/ica_datoja/stiffness.dat") plottaus(x) # c) huberM<-function(data) { n <- dim(data)[1] k <- dim(data)[2] q <- pchisq(k+1,k) c <- qchisq(q,k) c2 <- pchisq(c,k+2)+(c/k)*(1-q) bdp <- min(1/c,1-k/c) MAX_ITER <- 250; # Max number of iteration EPS <- 1.0e-6; # Iteration accuracy # STARTING VALUES t <- apply(data,2,mean) C <- cov(data) R<-chol(solve(C)) iter <- 0 d1<-1 d2<-1 while (d1 > EPS & d2 > EPS ) { # SCATTER STEP u<-NULL data2 <- sweep(data,2,t,"-") s <- diag(data2%*%(t(R)%*%R)%*%t(data2)) u[s<=c] <- (1/c2) u[s>c] <- (s[s>c]/(c/c2))^{-1} C <- R%*%(t(data2)%*%sweep(data2,1,u,"*")/n)%*%t(R) R0 <- chol(solve(C)) R <- R0%*%R d1 = max(apply(abs(R0-diag(k)),1,sum)) # LOCATION STEP v <- NULL s <- diag(data2%*%(t(R)%*%R)%*%t(data2)) v[s<=c] <- 1 v[s>c] <- sqrt(c)/sqrt(s[s>c]) h <- apply(sweep(data2,1,v,"*"),2,mean)/mean(v); t <- t + h; d2 <-sqrt(t(h-t)%*%(h-t)) iter <- iter+1 if (iter > MAX_ITER) break } C <- t(R)%*%R C <- solve(C) list (loc=t,cov=C,iter=iter) } x<- as.matrix(x) huber<-huberM(x) S.h<-huber$cov S.h X11() plottaus(x,S.h) # 8. source("http://users.jyu.fi/~slahola/files/ica_demoja/apufunktiot.R") sti <- read.table("http://users.jyu.fi/~slahola/files/ica_datoja/stiffness.dat",header=TRUE) x <- as.matrix(sti) m1 <- apply(x,2,mean) xcent1 <- sweep(x,2,m1,"-") ; xcent1 <- as.matrix(xcent1) S <- cov(x) z1 <- xcent1 %*% solve(mat.sqrt(S)) r1 <- diag(z1 %*% t(z1)) a <- huberM(x) m2 <- a$loc xcent2 <- sweep(x,2,m2,"-") ; xcent2 <- as.matrix(xcent2) z2 <- xcent2 %*% solve(mat.sqrt(a$cov)) r2 <- diag(z2 %*% t(z2)) plot(r1,r2) abline(h = sqrt(qchisq(p=.975,df=2))) abline(v = sqrt(qchisq(p=.975,df=2)))