y <- as.vector(xx %*% vv) # yy <- y %o% vv # sum((yy-xx)^) cat("\n") v0 <- rep(0,ncol(xx)-) # print(vv88(v0)) a <- optim(v0,rss88,control=list(trace=t

Transcription

1 Copyright (c) 004,005 Hidetoshi Shimodaira x... xp X = xn... xnp {{ p n x () = = [x,..., xp] x (n) x (i) xj X X n n nx :43:33 shimo R dat <- scale(dat,center=t,scale=f) dat <- scale(dat,scale=f) # run0087.r # dat <- read.table("dat000.txt") # 0 cat("# \n") dim(dat); names(dat) cat("# \n") mean(dat); apply(dat,,var) cat("# \n") xx <- scale(dat,scale=f) # cat("# \n") apply(xx,,mean); apply(xx,,var) plot87 <- function(x,y,dat) { plot(dat[,x],dat[,y],type="n",xlab=x,ylab=y) text(dat[,x],dat[,y],rownames(dat)) invisible(cbind(dat[,x],dat[,y])) pairs(xx) dev.copyeps(file="run0087-s0.eps") plot87("","",xx) dev.copyeps(file="run0087-s.eps") plot87("","",xx) dev.copyeps(file="run0087-s.eps") > source("run0087.r",print=t) # [] 47 0 [] "" "" "" "" "" "" "" [8] "" "" "Rikon" # Rikon Rikon # # e e e e e-5 Rikon e e e e e Rikon v v = v v =., vp p vj =. j= i v yi x... xp X = xn... xnp n x () = = [x,..., xp] x (n) 0 5 {{ p y yi = x (i) v, i =,..., n y =. = Xv Rikon x (i) yiv x (i) yiv i =,..., n v = n x (i) yiv i= yn Fukui Niigata Gifu Fukushima Iwate Nagano Shizuoka Gumma Mie Yamanashi Ishikawa WakayamaOkayama Saitama Kagawa Hyogo Hokkaido Miyazaki Fukuoka run0087-s run0087-s0 Okinawa Saitama Hyogo Hokkaido Fukuoka Okinawa Ishikawa Gifu Fukushima Iwate Shizuoka Mie Gumma Yamanashi Kagawa Miyazaki Okayama Nagano Wakayama run0087-s Fukui Niigata optim v = (v, v,..., vp ) vp # run0088.r # # dat xx <- scale(dat,scale=f) # vv88 <- function(v) { vp = v vp vp <- sqrt(-sum(v*v)) # c(v,vp) # rss88 <- function(v) { # vv <- vv88(v) # 4

2 y <- as.vector(xx %*% vv) # yy <- y %o% vv # sum((yy-xx)^) cat("\n") v0 <- rep(0,ncol(xx)-) # print(vv88(v0)) a <- optim(v0,rss88,control=list(trace=t,parscale=rep(0.,9)),method="bfgs") cat("\n") v <- a$par # vv <- vv88(v) y <- xx %*% vv # print(vv); print(y) plot87("","y",data.frame(xx,y)) dev.copyeps(file="run0088-s.eps") > source("run0088.r") [] initial value iter 0 value iter 0 value final value converged [] [7] [,] Hokkaido Iwate Miyazaki Okinawa There were 50 or more warnings (use warnings() to see the first 50) 5 y anagawa Okinawa Saitama Hokkaido Fukuoka Hyogo Miyazaki Gumma Wakayama Okayama Shizuoka Yamanashi Mie Kagawa Ishikawa Gifu run0088-s Fukushima Nagano Iwate Niigata Fukui Yamaga y:.3 v = n x (i) yiv X yv i= = tr((x yv ) (X yv )) A, B tr(ab) = tr(ba) y = Xv = tr(x X) y Xv + y y = tr(x X) y y y # run0089.r # # dat xx <- scale(dat,scale=f) # rss89 <- function(v) { # vv <- vv88(v) # 6 y <- as.vector(xx %*% vv) # sum(y*y) v0 <- rep(0,ncol(xx)-) # a <- optim(v0,rss89,control=list(trace=t,parscale=rep(0.,9),fnscale=-), method="bfgs") v <- a$par # vv <- vv88(v) y <- xx %*% vv # print(vv); print(y) plot(y,y); abline(0,) dev.copyeps(file="run0089-s.eps") y > source("run0089.r") initial value iter 0 value iter 0 value final value converged [] [7] [,] Hokkaido Iwate Miyazaki Okinawa There were 50 or more warnings (use warnings() to see the first 50) y run0089-s y: y: σ x = n x,..., σ xp = n xp xj xj,..., j =,..., p σxj R dat <- scale(dat,center=t,scale=t) dat <- scale(dat) # run0090.r # # dat cat("# \n") xx <- scale(dat) # cat("# \n") print(apply(xx,,mean)); print(apply(xx,,var)) v0 <- rep(0,ncol(xx)-) # a <- optim(v0,rss89,control=list(trace=t,parscale=rep(0.,9),fnscale=-), method="bfgs") v3 <- a$par # 7 8

3 vv3 <- vv88(v3) y3 <- xx %*% vv3 # print(vv3); print(y3) plot87("y","y3",data.frame(y,y3)) dev.copyeps(file="run0090-s.eps") > source("run0090.r") # # e e e e e-6 Rikon e e e e e-6 Rikon initial value iter 0 value final value converged [] [6] [,] Hokkaido Iwate Miyazaki Okinawa Warning messages: : NaNs produced in: sqrt( - sum(v * v)) : NaNs produced in: sqrt( - sum(v * v)) 3: NaNs produced in: sqrt( - sum(v * v)) 4: NaNs produced in: sqrt( - sum(v * v)) 5: NaNs produced in: sqrt( - sum(v * v)) 9 y Nagano Fukushima Gifu Iwate Fukui Niigata Miyazaki Shizuoka Kagawa Okayama Gumma Wakayama Yamanashi Mie Ishikawa y run0090-s Fukuoka Hokkaido Saitama Hyogo Okinawa y: y3: y y3.5 y = Xv, v = y v X Xv, v v = f(v, λ) = v X Xv λ(v v ) f v = X Xv λv = 0, X Xv = λv, v = f λ = v v = 0 X X () v λ v X Xv y = v X Xv = λv v = λ v y 0 X X n X X λ y n y n X Xv = λ, n y = λ n X X n X X # run009.r # # dat cat(" \n") xx <- scale(dat,scale=f) # cv <- var(xx) # print(cv[:5,:5]) cat(" \n") vv4 <- eigen(cv)$vectors[,] y4 <- xx %*% vv4 # print(vv4); print(y4) plot(y,y4); abline(0,) dev.copyeps(file="run009-s.eps") cat(" \n") xx <- scale(dat) # cv <- var(xx) # print(cv[:5,:5]) cat(" \n") vv5 <- eigen(cv)$vectors[,] y5 <- xx %*% vv5 # print(vv5); print(y5) plot(y3,y5); abline(0,) dev.copyeps(file="run009-s.eps") > source("run009.r") [] [7] [,] Hokkaido Iwate Miyazaki Okinawa [] [6] [,] Hokkaido Iwate Miyazaki Okinawa

4 V = (v,..., vp) Y = (y,..., y p ) y y Y = XV V V = Ip V p x, x,..., xp p y, y,..., y p v v v v, v v3 v,..., vr vr y y3 run009-s run009-s xx: y4: y xx: y5: y3 optim eigen eigen eigen eigen vj λj λj = λj+ = = λj+s s vj, vj+,..., vj+s n y j = λj λj y j j k n y jy k = v j n (X X)vk = v j(λkvk) = λk(v jvk) = 0 n Y Y = V ( n X XV ) = V (V Λ) = (V V )Λ = Λ. (principal component analysys) PCA (principal component) PC? ( ) X y, y,..., y p y j = Xvj X n X X λ λ λp 0 v, v,..., vp 3 Λ = diag(λ,..., λp) n X XV = V Λ s v, v,..., vs λ + + λs = λ + + λp V = (v,..., vp) V V = Ip # run009.r n y + n y p = λ + + λp = n x + n xp # # dat xx <- scale(dat) # cv <- var(xx) # 4 ei <- eigen(cv) # cat("\n"); print(ei) yy <- xx %*% ei$vectors # cat(" (j=,,3)\n"); print(yy[:5,:3]); cat("......");print(yy[43:47,:3]) cat("\n"); print(cumsum(ei$values)/sum(ei$values)) plot87(,,yy); dev.copyeps(file="run009-s.eps") plot87(3,,yy); dev.copyeps(file="run009-s3.eps") > source("run009.r") $values [] [7] $vectors [,] [,] [,3] [,4] [,5] [,6] [,] [,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [0,] [,7] [,8] [,9] [,0] [,] [,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [0,] Hokkaido Iwate [,] [,] [,3] Miyazaki Okinawa [] [8] Miyazaki Wakayama Hokkaido Kagawa Okayama Fukuoka Hyogo Iwate NaganoYamanashi Mie Niigata FukushimaIshikawa Gumma Fukui Gifu Shizuoka Saitama Okinawa run009-s run009-s r r r v,..., vr y j = Xvj, V r = [v,..., vr], V rv r = Ir 0 4 Miyazaki Wakayama Hokkaido Kagawa Fukuoka Okayama Hyogo Iwate Yamanashi Nagano Mie Ishikawa Fukushima Niigata Gumma Fukui Shizuoka Gifu Saitama Okinawa (j=,,3) yij = x (i) vj, i =,..., n, j =,..., r [,] [,] [,3] 5 6

5 n r = x (i) yijv j i= j= n = x (i) (Ip V rv r) i= = tr(x(ip V rv r) X ) = tr(xx XV rv rx ) = tr(x X) tr(v rx XV r) n p n r = x ij yij i= j= i= j= r tr(v rx XV r), V rv r = Ir r r Λ r r r r f(v r, Λ) = v ix Xvi λii(v ivi ) λijv ivj i= i= i= j>i = tr (V rx XV r Λ(V rv r Ir)) f r = X f Xvi λijvj, = X XV r V rλ vi V j= r Λ r r Q Q ΛQ = diag(λ,..., λr) V r V rq X Xvi = λivi, i =,..., r X X v,..., vr = tr(x X) (λ + + λr) λ,..., λr r v,..., vr r r.3 zj = y j, j =,..., p λj z () Z = [z,..., zp] =. z (n) x (i), i =,..., n z (i) Z = Y Λ / Λ / = diag(λ /,..., λ / p ) Z n Z Z = Ip n Z Z = Λ / ( n Y Y )Λ / = Λ / ΛΛ / = Ip xj zk n x jzk B B = n X Z, B = [b,..., bp] = xj, j =,..., p b (j) n Z Z = Ip B X = ZB b (). b (p) xj = Z(b (j) ) xj z,..., zp b (j) (i) (ii) (i) n z z (ii) p b b p X X = ZB = zb + + zpb p r X zb + + zrb r r = X 7 8 # run0093.r # # dat xx <- scale(dat) # cv <- var(xx) # ei <- eigen(cv) # yy <- xx %*% ei$vectors # lam <- diag(/sqrt(ei$values)) # Lambda^{-/ zz <- yy %*% lam # n <- nrow(xx) bb <- crossprod(xx,zz)/(n-) # =t(xx) %*% zz /(n-) cat(" Y (i=:5, j=:3)\n"); print(yy[:5,:3]); cat("\n"); print(cumsum(ei$values)/sum(ei$values)) cat(" Z (i=:5, j=:3)\n"); print(zz[:5,:3]); cat(" B (j=:3)\n"); print(bb[,:3]); plot87(,,zz); dev.copyeps(file="run0093-z.eps") plot87(3,,zz); dev.copyeps(file="run0093-z3.eps") plot87(,,bb); dev.copyeps(file="run0093-b.eps") plot87(3,,bb); dev.copyeps(file="run0093-b3.eps") plot(xx,zz %*% t(bb)); abline(0,); dev.copyeps(file="run0093-zzbb.eps") plot(xx,zz[,:3] %*% t(bb[,:3])); abline(0,); dev.copyeps(file="run0093-zzbb3.eps") > source("run0093.r") Y (i=:5, j=:3) [,] [,] [,3] Hokkaido Iwate [] [8] Z (i=:5, j=:3) 9 [,] [,] [,3] Hokkaido Iwate B (j=:3) [,] [,] [,3] Rikon > xx[:5,:5] Hokkaido Iwate > (zz %*% t(bb))[:5,:5] Hokkaido Iwate > (zz[,:3] %*% t(bb[,:3]))[:5,:5] Hokkaido Iwate

6 0 z,z3: Z b,b3: B zz %*% t(bb) Miyazaki Wakayama Hokkaido Kagawa Okayama Fukuoka Hyogo Iwate NaganoYamanashi Mie Niigata FukushimaIshikawa Gumma Fukui Gifu Shizuoka Saitama Okinawa run0093-z Rikon run0093-b xx run0093-zzbb zz[, :3] %*% t(bb[, :3]) Miyazaki Wakayama Hokkaido Kagawa Fukuoka Okayama Iwate andoku 3 Hyogo Yamanashi Nagano Mie Ishikawa Fukushima Niigata Gumma Fukui Shizuoka Gifu Saitama Okinawa run0093-z3 Rikon run0093-b xx run0093-zzbb3 zzbb: =X =ZB zzbb3: =X = zb + + zrb r r = 3.4 # run0095.r # mybiplot <- function(x,y,choices=:,scale=c(,), col.arg=c(,),cex.arg=c(,),magnify=, xadj.arg=c(0.5,0.5),yadj.arg=c(0.5,0.5), xnames=null,ynames=null) { if(length(choices)!= ) stop("choices must be length ") if(length(scale)!= ) stop("scale must be length ") # x <- x[,choices] %*% diag(scale) # y <- y[,choices] %*% diag(/scale) x <- x[,choices] * rep(scale,rep(dim(x)[],)) y <- y[,choices] * rep(/scale,rep(dim(y)[],)) if(is.null(xnames)) nx <- dimnames(x)[[]] else nx <- as.character(xnames) if(is.null(ynames)) ny <- dimnames(y)[[]] else ny <- as.character(ynames) if(is.null(dimnames(x)[[]])) nd <- paste("pc",choices) else nd <- dimnames(x)[[]] rx <- range(x); ry <- range(y) oldpar <- par(pty="s") a <- min(rx/ry); yy <- y*a plot(x,xlim=rx*magnify,ylim=rx*magnify,type="n",xlab=nd[],ylab=nd[]) ly <- pretty(rx/a) ly[abs(ly) < e-0] <- 0 axis(3,at = ly*a,labels = ly) axis(4,at = ly*a,labels = ly) text(yy,ny,col=col.arg[],cex=cex.arg[],adj=yadj.arg) arrows(0,0,yy[,]*0.8,yy[,]*0.8,col=col.arg[],length=0.) text(x,nx,col=col.arg[],cex=cex.arg[],adj=yadj.arg) par(oldpar) invisible(list(x=x,y=y)) mybiplot(zz,bb); dev.copyeps(file="run0095-b.eps") mybiplot(zz,bb,choices=3:,scale=c(-,)); dev.copyeps(file="run0095-b3.eps") PC Wakayama Miyazaki Kagawa Okayama Iwate Nagano Yamanashi Mie Niigata Fukushima Ishikawa Gumma Fukui Gifu Shizuoka Rikon Fukuoka 0 PC Hokkaido Hyogo run0095-b Saitama Okinawa PC Wakayama Miyazaki Hokkaido Rikon Kagawa Okayama Fukuoka Hyogo Mie Yamanashi Nagano Iwate Gumma Fukushima Niigata Ishikawa Gifu Shizuoka Fukui Okinawa Saitama PC 3 run0095-b3 PC vs,,,, 65,, PC?,,,, PC3?,.5 princomp # run0096.r # princomp # dat cat("\n") a <- princomp(dat) print(summary(a)) biplot(a); dev.copyeps(file="run0096-b.eps") cat(" \n") a <- princomp(dat,cor=t) print(summary(a)) biplot(a); dev.copyeps(file="run0096-b.eps") > source("run0096.r") Importance of components: Comp. Comp. Comp.3 Comp.4 Comp.5 Standard deviation Proportion of Variance Cumulative Proportion Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation e-0 Proportion of Variance e-05 Cumulative Proportion e-0 Comp.0 Standard deviation.538e-0 Proportion of Variance.96080e-06 Cumulative Proportion e+00 Importance of components: Comp. Comp. Comp.3 Comp.4 Comp.5 Standard deviation Proportion of Variance Cumulative Proportion Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation Proportion of Variance Cumulative Proportion Comp.0 Standard deviation Proportion of Variance Cumulative Proportion

7 Comp Saitama Gifu Gumma Okinawa Shizuoka Mie Hyogo Fukui Wakayama Miyazaki NaganoYamanashi Kagawa Wakayama Miyazaki Hokkaido Niigata FukushimaOkayama Rikon Ishikawa Kagawa Okayama Fukuoka Rikon Iwate Hokkaido Fukuoka Hyogo Iwate Nagano Yamanashi Mie Niigata Fukushima Ishikawa Gumma Fukui Gifu Shizuoka Saitama Okinawa Comp. Comp. run0096-b run0096-b princomp() princomp(,cor=t) Comp for(i in year) { j <- paste("year",i,sep=""); x[[j]] <- k <- x0[x0$year == i,c("team",item)]; x[[j]]$team <- NULL; rownames(x[[j]]) <- team[as.character(k$team)]; colnames(x[[j]]) <- na[colnames(x[[j]])]; pc[[j]] <- princomp(x[[j]],cor=t); biplot(pc[[j]],main=paste("year =",i)) # ex <- list(year000=c(-,),year00=c(-,-), Year00=c(-,),Year003=c(-,)) # for(i in year) { j <- paste("year",i,sep=""); pc[[j]]$scores[,:] <- pc[[j]]$scores[,:] %*% diag(ex[[j]]); pc[[j]]$loadings[,:] <- pc[[j]]$loadings[,:] %*% diag(ex[[j]]); biplot(pc[[j]],main=j) # dev.copyeps(file="run0097.eps") # EPS par(mfrow=c(,),pty="s") summary().6 # run0097.r # dageki <- read.table("teamdageki.txt",header=t,sep="\t") # toushu <- read.table("teamtoushu.txt",header=t,sep="\t") # x0 <- data.frame(dageki,toushu) # team <- c("kyojin", "Yakult", "Yokohama", "Chunichi", "Hanshin", "", "Lotte", "Nichiham", "Seibu", "Kintetsu", "Orix", "Daiei","Taiyo") names(team) <- c(" ", " ", " ", " ", " ", " ", "", " ", " ", " ", "", " "," ") item <- c("daritsu","choudaritsu","shutsuruiritsu","shubiritsu", "Touruiboushiritsu","Shouritsu","Bougyoritsu") # na <- substr(item,,nchar(item)-5) # "ritsu" names(na) <- item year <- 000:003 # par(mfrow=c(,),pty="s") # x x <- list(); pc <- list(); 5 6 Comp Year Yakult Touruiboushi Kyojin Hansh Shubi Seibu Shou Yokohama Chunichi Daiei Chouda Da Shutsurui Orix Lotte hiham Kintetsu Bougyo Comp. 0 4 Comp Year anshin Chunichi Yokohama Touruiboushi Yakult Shubi Seibu Orix Comp. Daiei Kyojin Shou Da Shutsu Chou Nichiham Lotte Bougyo Kintetsu (SVD) (singular value decomposition) x... xp X = xn... xnp {{ p X = UDV n x () = = [x,..., xp] x (n) = duv + + dpupv p Comp Year Orix okohama Nichiham Lotte Bougyo Chunichi Yakult Hanshin Touruiboushi Kyojin Shubi Daiei Kintetsu Seib Da Choud Shutsu Sho 4 0 Comp Year Chunichi Shubi Yokohama Kyojin HanshinSho Nichiham Lotte Yakult Seibu Kintetsu Touruiboush Shutsur Da Daie Bougyo Chouda Orix U = [u,..., up], V = [v,..., vp] n p p p d 0 D =..., d dp 0 0 dp X Σ = n X X = V n D V Comp. Comp. run (Da) (Chouda) (Shutsurui) (Shubi) (Touruiboshi) (Shou) (Bougyo) 000 =00 =00 = 003 = ( v,..., vp y j = Xvj, j =,..., p λ = n d,..., λp = n d p Y = [y,..., y p ] = XV = UD Λ = n D Z = [z,..., zp] = Y Λ / = n U B = n X Z = n V D 7 8

8 # run0098.r # # dat xx <- scale(dat) # a <- svd(xx) # rownames(a$u) <- rownames(xx) rownames(a$v) <- colnames(xx) cat("\n") print(names(a)) # cat("d ", length(a$d),"\n") cat("u ", dim(a$u),"\n") cat("v ", dim(a$v),"\n") cat("xx[:5,:5]\n") print(xx[:5,:5]) cat("(a$u %*% diag(a$d) %*% t(a$v))[:5,:5]\n") print((a$u %*% diag(a$d) %*% t(a$v))[:5,:5]) cat(" cumsum(a$d^)/sum(a$d^)\n") print(cumsum(a$d^)/sum(a$d^)) n <- nrow(xx) zz <- sqrt(n-)*a$u # bb <- (/sqrt(n-))* a$v %*% diag(a$d) # mybiplot(zz,bb); dev.copyeps(file="run0098-b.eps") > source("run0098.r") [] "d" "u" "v" d 0 u 47 0 v 0 0 xx[:5,:5] Hokkaido Iwate (a$u %*% diag(a$d) %*% t(a$v))[:5,:5] Hokkaido Iwate cumsum(a$d^)/sum(a$d^) [] [8] PC Wakayama Miyazaki Hokkaido Kagawa Rikon Okayama Fukuoka Iwate Nagano Yamanashi Mie NiigataFukushima Ishikawa Gumma Fukui Gifu Shizuoka PC Hyogo run0098-b Saitama Okinawa lambda: n X X z: Z b: B mybiplot (run0095.r)