Copyright (c) 2004,2005 Hidetoshi Shimodaira :43:33 shimo X = x x 1p x n1... x np } {{ } p n = x (1) x (n) = [x 1,..

Size: px

Start display at page:

しほこももき
5 years ago
Views:

1 Copyright (c) 2004,2005 Hidetoshi Shimodaira :43:33 shimo X = x x 1p x n1... x np } {{ } p n = x (1) x (n) = [x 1,..., x p ] x (i) x j X X 1 1 n n1 nx R dat <- scale(dat,center=t,scale=f) dat <- scale(dat,scale=f) # run0087.r # dat <- read.table("dat0002.txt") # 10 cat("# \n") dim(dat); names(dat) cat("# \n") mean(dat); apply(dat,2,var) cat("# \n") xx <- scale(dat,scale=f) # cat("# \n") apply(xx,2,mean); apply(xx,2,var) 1

2 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.copy2eps(file="run0087-s0.eps") plot87("zouka","ninzu",xx) dev.copy2eps(file="run0087-s1.eps") plot87("x65sai","tomo",xx) dev.copy2eps(file="run0087-s2.eps") > source("run0087.r",print=t) # [1] [1] "Zouka" "Ninzu" "Kaku" "Tomo" "Tandoku" "X65Sai" "Kfufu" [8] "Ktan" "Konin" "Rikon" # Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon # # Zouka Ninzu Kaku Tomo Tandoku e e e e e-15 X65Sai Kfufu Ktan Konin Rikon e e e e e-17 Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon

3 Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon run0087-s0 Ninzu Yamagata Fukui Toyama Niigata Saga Gifu Fukushima Shiga Akita Tottori Tochigi Ibaraki Iwate Nara Shimane Nagano Shizuoka Gumma Mie Aomori Yamanashi Ishikawa Kumamoto Miyagi Tokushima WakayamaOkayama Saitama Kagawa Aichi Nagasaki Hyogo Chiba Ooita Miyazaki Ehime Yamaguchi Hiroshima Kyoto Fukuoka Kanagawa Osaka Kochi Kagoshima Hokkaido Tokyo Okinawa Tomo Fukui Toyama Niigata Tottori Nagano Ishikawa Gifu Fukushima Saga Iwate Shizuoka Tochigi Shiga Mie Gumma Ibaraki Yamanashi Kagawa Miyazaki Tokushima Okayama Kumamoto Aomori Aichi Kochi Yamaguchi Hiroshima Ooita Miyagi Saitama Ehime Nagasaki Wakayama Chiba Kagoshima Kyoto HyogoNara Hokkaido Fukuoka Okinawa Kanagawa Osaka Tokyo Akita Yamagata Shimane Zouka run0087-s X65Sai run0087-s2 3

4 1.2 v v = 1 v 1 v =., v p p vj 2 = 1. j=1 i v y i X = x x 1p x n1... x np } {{ } p n = x (1) x (n) = [x 1,..., x p ] y i = x (i) v, i = 1,..., n y =. = Xv x (i) y i v x (i) y i v 2 i = 1,..., n v y 1 y n = n x (i) y i v 2 i=1 optim v = 1 (v 1, v 2,..., v p 1 ) v p v p = 1 v 2 1 v 2 p 1 # run0088.r # # dat xx <- scale(dat,scale=f) # vv88 <- function(v) { vp <- sqrt(1-sum(v*v)) # c(v,vp) # } rss88 <- function(v) { # vv <- vv88(v) # 4

5 y <- as.vector(xx %*% vv) # yy <- y %o% vv # sum((yy-xx)^2) } cat(" \n") v0 <- rep(0,ncol(xx)-1) # print(vv88(v0)) a <- optim(v0,rss88,control=list(trace=t,parscale=rep(0.1,9)),method="bfgs") cat(" \n") v1 <- a$par # vv1 <- vv88(v1) y1 <- xx %*% vv1 # print(vv1); print(y1) plot87("x65sai","y1",data.frame(xx,y1)) dev.copy2eps(file="run0088-s1.eps") > source("run0088.r") [1] initial value iter 10 value iter 20 value final value converged [1] [7] [,1] Hokkaido Aomori Iwate Miyagi Akita Kumamoto Ooita Miyazaki Kagoshima Okinawa There were 50 or more warnings (use warnings() to see the first 50) 5

6 y Tokyo anagawa Osaka Okinawa Saitama Chiba Hokkaido Fukuoka Kyoto Hyogo Aichi Nara Hiroshima Kagoshima Miyagi Ehime Miyazaki Nagasaki Yamaguchi Shiga Ibaraki Gumma Ooita Wakayama Kochi Okayama Tochigi Shizuoka Yamanashi Mie Kagawa Aomori Kumamoto Ishikawa Tokushima Gifu Fukushima Nagano Iwate Saga Niigata Tottori Toyama Fukui Akita Shimane Yamaga X65Sai run0088-s1 y1: X65Sai 1.3 v = n x (i) y i v 2 i=1 X yv = tr((x yv ) (X yv )) A, B tr(ab) = tr(ba) = tr(x X) 2y Xv + y y y = Xv = tr(x X) y y y 2 # run0089.r # # dat xx <- scale(dat,scale=f) # rss89 <- function(v) { # vv <- vv88(v) # 6

7 y <- as.vector(xx %*% vv) # sum(y*y) } v0 <- rep(0,ncol(xx)-1) # a <- optim(v0,rss89,control=list(trace=t,parscale=rep(0.1,9),fnscale=-1), method="bfgs") v2 <- a$par # vv2 <- vv88(v2) y2 <- xx %*% vv2 # print(vv2); print(y2) plot(y1,y2); abline(0,1) dev.copy2eps(file="run0089-s1.eps") > source("run0089.r") initial value iter 10 value iter 20 value final value converged [1] [7] [,1] Hokkaido Aomori Iwate Miyagi Akita Kumamoto Ooita Miyazaki Kagoshima Okinawa There were 50 or more warnings (use warnings() to see the first 50) 7

8 y y1 run0089-s1 y1: y2: σ 2 x 1 = 1 n 1 x 1 2,..., σ 2 x p = 1 n 1 x p 2 x j 1 σ xj x j,..., j = 1,..., p R dat <- scale(dat,center=t,scale=t) dat <- scale(dat) # run0090.r # # dat cat("# \n") xx <- scale(dat) # cat("# \n") print(apply(xx,2,mean)); print(apply(xx,2,var)) v0 <- rep(0,ncol(xx)-1) # a <- optim(v0,rss89,control=list(trace=t,parscale=rep(0.1,9),fnscale=-1), method="bfgs") v3 <- a$par # 8

9 vv3 <- vv88(v3) y3 <- xx %*% vv3 # print(vv3); print(y3) plot87("y2","y3",data.frame(y2,y3)) dev.copy2eps(file="run0090-s1.eps") > source("run0090.r") # # Zouka Ninzu Kaku Tomo Tandoku e e e e e-16 X65Sai Kfufu Ktan Konin Rikon e e e e e-16 Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon initial value iter 10 value final value converged [1] [6] [,1] Hokkaido Aomori Iwate Miyagi Akita Kumamoto Ooita Miyazaki Kagoshima Okinawa Warning messages: 1: NaNs produced in: sqrt(1 - sum(v * v)) 2: NaNs produced in: sqrt(1 - sum(v * v)) 3: NaNs produced in: sqrt(1 - sum(v * v)) 4: NaNs produced in: sqrt(1 - sum(v * v)) 5: NaNs produced in: sqrt(1 - sum(v * v)) 9

10 y amagata Tokyo Osaka Okinawa Kanagawa Fukuoka Hokkaido Saitama Chiba Aichi Hyogo Kyoto Hiroshima Miyazaki Miyagi Kagoshima Nara Shizuoka Shiga Ibaraki Tochigi Okayama Gumma Kochi Nagasaki Kagawa Wakayama Yamaguchi Ooita Ehime Yamanashi Mie Ishikawa Kumamoto Aomori Tokushima Nagano Fukushima Gifu Saga Tottori Iwate Toyama Fukui Niigata Akita Shimane y2 run0090-s1 y2: y3: y2 y3 1.5 y = Xv, v = 1 y 2 v X Xv, v v = 1 f(v, λ) = v X Xv λ(v v 1) f v = 2X Xv 2λv = 0, X Xv = λv, v = 1 f λ = v v 1 = 0 X X ( ) v λ v X Xv y 2 = v X Xv = λv v = λ v y 2 10

11 X X 1 n 1 X X λ y 1 n 1 y 2 1 n 1 X 1 Xv = λ, n 1 y 2 = λ 1 n 1 X X 1 n 1 X X # run0091.r # # dat cat(" \n") xx1 <- scale(dat,scale=f) # cv1 <- var(xx1) # print(cv1[1:5,1:5]) cat(" \n") vv4 <- eigen(cv1)$vectors[,1] y4 <- xx1 %*% vv4 # print(vv4); print(y4) plot(y2,y4); abline(0,1) dev.copy2eps(file="run0091-s1.eps") cat(" \n") xx2 <- scale(dat) # cv2 <- var(xx2) # print(cv2[1:5,1:5]) cat(" \n") vv5 <- eigen(cv2)$vectors[,1] y5 <- xx2 %*% vv5 # print(vv5); print(y5) plot(y3,y5); abline(0,1) dev.copy2eps(file="run0091-s2.eps") > source("run0091.r") Zouka Ninzu Kaku Tomo Tandoku Zouka Ninzu Kaku Tomo Tandoku

12 [1] [7] [,1] Hokkaido Aomori Iwate Miyagi Akita Kumamoto Ooita Miyazaki Kagoshima Okinawa Zouka Ninzu Kaku Tomo Tandoku Zouka Ninzu Kaku Tomo Tandoku [1] [6] [,1] Hokkaido Aomori Iwate Miyagi Akita Kumamoto Ooita Miyazaki Kagoshima Okinawa

13 y y y2 run0091-s y3 run0091-s2 xx1: y4: y2 xx2: y5: y3 optim eigen eigen eigen eigen (principal component analysys) PCA (principal component) PC? ( ) X y 1, y 2,..., y p y j = Xv j X 1 n 1 X X λ 1 λ 2 λ p 0 v 1, v 2,..., v p 13

14 V = (v 1,..., v p ) Y = (y 1,..., y p ) Y = XV V V = I p V p x 1, x 2,..., x p p y 1, y 2,..., y p v 1 v 1 v 2 v 1, v 2 v 3 v 1,..., v r 1 v r v j λ j λ j = λ j+1 = = λ j+s 1 s v j, v j+1,..., v j+s 1 1 n 1 y j 2 = λ j λ j y j j k 1 n 1 y jy k = v 1 j n 1 (X X)v k = v j(λ k v k ) = λ k (v jv k ) = 0 1 Y Y = V ( 1 n 1 n 1 X XV ) = V (V Λ) = (V V )Λ = Λ Λ = diag(λ 1,..., λ p ) 1 n 1 X XV = V Λ s v 1, v 2,..., v s = λ λ s λ λ p V = (v 1,..., v p ) V V = I p 1 y n y n 1 p 2 = λ λ p = 1 x n x n 1 p 2 # run0092.r # # dat xx <- scale(dat) # cv <- var(xx) # 14

15 ei <- eigen(cv) # cat(" \n"); print(ei) yy <- xx %*% ei$vectors # cat(" (j=1,2,3)\n"); print(yy[1:5,1:3]); cat("......");print(yy[43:47,1:3]) cat(" \n"); print(cumsum(ei$values)/sum(ei$values)) plot87(1,2,yy); dev.copy2eps(file="run0092-s12.eps") plot87(3,2,yy); dev.copy2eps(file="run0092-s32.eps") > source("run0092.r") $values [1] [7] $vectors [,1] [,2] [,3] [,4] [,5] [,6] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [,7] [,8] [,9] [,10] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] (j=1,2,3) [,1] [,2] [,3] 15

16 Hokkaido Aomori Iwate Miyagi Akita [,1] [,2] [,3] Kumamoto Ooita Miyazaki Kagoshima Okinawa [1] [8] Kagoshima Kagoshima Shimane Akita Kochi Yamaguchi Ehime Ooita Miyazaki Wakayama Nagasaki Hokkaido Tokushima Kumamoto Kagawa Hiroshima Okayama Kyoto Fukuoka Tottori Aomori Hyogo Iwate NaganoYamanashi Mie Saga Nara Niigata FukushimaIshikawa Gumma Toyama Miyagi amagata Fukui Gifu Shizuoka Chiba Tochigi Ibaraki Aichi Saitama Shiga Tokyo Osaka Kanagawa Okinawa Tokyo Kochi Yamaguchi Ehime Ooita Miyazaki Wakayama Hokkaido Nagasaki Shimane Tokushima Kumamoto Hiroshima Kagawa Kyoto Fukuoka AkitaOkayama Aomori Tottori Osaka Hyogo Iwate Yamanashi Nagano Mie Saga Nara Ishikawa Fukushima Niigata Gumma Kanagawa Miyagi Toyama Yamagata Fukui Shizuoka ChibaGifu Tochigi Aichi Ibaraki Saitama Okinawa Shiga run0092-s run0092-s r r r v 1,..., v r y j = Xv j, V r = [v 1,..., v r ], V rv r = I r y ij = x (i) v j, i = 1,..., n, j = 1,..., r 16

17 n r = x (i) y ij v j 2 r = i=1 j=1 n x (i) (I p V r V r) 2 i=1 = tr(x(i p V r V r) 2 X ) = tr(xx XV r V rx ) = tr(x X) tr(v rx XV r ) n p n r = x 2 ij i=1 j=1 i=1 j=1 y 2 ij tr(v rx XV r ), V rv r = I r r r Λ r r f(v r, Λ) = v ix Xv i λ ii (v iv i 1) 2 i=1 i=1 = tr (V rx XV r Λ(V rv r I r )) f v i = 2X Xv i 2 r λ ij v j, j=1 r i=1 r λ ij v iv j j>i f V r = 2X XV r 2V r Λ Λ r r Q V r V r Q Q ΛQ = diag(λ 1,..., λ r ) X Xv i = λ i v i, i = 1,..., r X X v 1,..., v r = tr(x X) (λ λ r ) λ 1,..., λ r r v 1,..., v r r r 2.3 z j = y j λj, j = 1,..., p 17

18 Z = [z 1,..., z p ] = x (i), i = 1,..., n z (i) Z = Y Λ 1/2 z (1). z (n) Λ 1/2 = diag(λ 1/2 1,..., λ 1/2 p ) Z 1 n 1 Z Z = I p 1 n 1 Z Z = Λ 1/2 ( 1 Y Y )Λ 1/2 = Λ 1/2 ΛΛ 1/2 = I n 1 p x j z k 1 n 1 x jz k B B = 1 n 1 X Z, B = [b 1,..., b p ] = x j, j = 1,..., p b (j) 1 n 1 Z Z = I p B X = ZB x j = Z(b (j) ) x j z 1,..., z p b (j) (i) (ii) (i) n z 1 z 2 (ii) p b 1 b 2 p X b (1). b (p) X = ZB = z 1 b z p b p r X z 1 b z r b r r = 2 X 18

19 # run0093.r # # dat xx <- scale(dat) # cv <- var(xx) # ei <- eigen(cv) # yy <- xx %*% ei$vectors # lam2 <- diag(1/sqrt(ei$values)) # Lambda^{-1/2} zz <- yy %*% lam2 # n <- nrow(xx) bb <- crossprod(xx,zz)/(n-1) # =t(xx) %*% zz /(n-1) cat(" Y (i=1:5, j=1:3)\n"); print(yy[1:5,1:3]); cat(" \n"); print(cumsum(ei$values)/sum(ei$values)) cat(" Z (i=1:5, j=1:3)\n"); print(zz[1:5,1:3]); cat(" B (j=1:3)\n"); print(bb[,1:3]); plot87(1,2,zz); dev.copy2eps(file="run0093-z12.eps") plot87(3,2,zz); dev.copy2eps(file="run0093-z32.eps") plot87(1,2,bb); dev.copy2eps(file="run0093-b12.eps") plot87(3,2,bb); dev.copy2eps(file="run0093-b32.eps") plot(xx,zz %*% t(bb)); abline(0,1); dev.copy2eps(file="run0093-zzbb.eps") plot(xx,zz[,1:3] %*% t(bb[,1:3])); abline(0,1); dev.copy2eps(file="run0093-zzbb3.eps") > source("run0093.r") Y (i=1:5, j=1:3) [,1] [,2] [,3] Hokkaido Aomori Iwate Miyagi Akita [1] [8] Z (i=1:5, j=1:3) 19

20 [,1] [,2] [,3] Hokkaido Aomori Iwate Miyagi Akita B (j=1:3) [,1] [,2] [,3] Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Rikon > xx[1:5,1:5] Zouka Ninzu Kaku Tomo Tandoku Hokkaido Aomori Iwate Miyagi Akita > (zz %*% t(bb))[1:5,1:5] Zouka Ninzu Kaku Tomo Tandoku Hokkaido Aomori Iwate Miyagi Akita > (zz[,1:3] %*% t(bb[,1:3]))[1:5,1:5] Zouka Ninzu Kaku Tomo Tandoku Hokkaido Aomori Iwate Miyagi Akita

21 Kagoshima Kagoshima Kochi Kochi Shimane Akita Yamaguchi Ehime Ooita Miyazaki Wakayama Nagasaki Hokkaido Tokushima Kumamoto Kagawa Hiroshima Okayama Kyoto Fukuoka Tottori Aomori Hyogo Iwate NaganoYamanashi Mie Saga Nara Niigata FukushimaIshikawa Gumma Toyama Miyagi amagata Fukui Gifu Shizuoka Chiba Tochigi Ibaraki Aichi Saitama Shiga Tokyo Osaka Kanagawa Okinawa Tokyo Yamaguchi Ehime Ooita Miyazaki Wakayama Hokkaido Nagasaki Shimane Tokushima Kumamoto Hiroshima Kagawa Kyoto Fukuoka AkitaOkayama Aomori Tottori Osaka Hyogo Iwate Yamanashi Nagano Mie Saga Nara Ishikawa Fukushima Niigata Gumma Kanagawa Miyagi Toyama Yamagata Fukui Shizuoka ChibaGifu Tochigi Aichi Ibaraki Saitama Okinawa Shiga run0093-z run0093-z X65Sai Tomo Kfufu Ktan Tandoku Rikon Kaku Konin andoku Ktan Kfufu X65Sai Rikon Tomo Konin Kaku Ninzu Zouka Ninzu Zouka run0093-b run0093-b32 3 zz %*% t(bb) zz[, 1:3] %*% t(bb[, 1:3]) xx run0093-zzbb xx run0093-zzbb3 21

22 z12,z32: Z b12,b32: B zzbb: =X =ZB zzbb3: =X = z 1 b z r b r r = # run0095.r # mybiplot <- function(x,y,choices=1:2,scale=c(1,1), col.arg=c(1,2),cex.arg=c(1,1),magnify=1, xadj.arg=c(0.5,0.5),yadj.arg=c(0.5,0.5), xnames=null,ynames=null) { if(length(choices)!= 2) stop("choices must be length 2") if(length(scale)!= 2) stop("scale must be length 2") # x <- x[,choices] %*% diag(scale) # y <- y[,choices] %*% diag(1/scale) x <- x[,choices] * rep(scale,rep(dim(x)[1],2)) y <- y[,choices] * rep(1/scale,rep(dim(y)[1],2)) if(is.null(xnames)) nx <- dimnames(x)[[1]] else nx <- as.character(xnames) if(is.null(ynames)) ny <- dimnames(y)[[1]] else ny <- as.character(ynames) if(is.null(dimnames(x)[[2]])) nd <- paste("pc",choices) else nd <- dimnames(x)[[2]] 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[1],ylab=nd[2]) ly <- pretty(rx/a) ly[abs(ly) < 1e-10] <- 0 axis(3,at = ly*a,labels = ly) axis(4,at = ly*a,labels = ly) text(yy,ny,col=col.arg[2],cex=cex.arg[2],adj=yadj.arg) arrows(0,0,yy[,1]*0.8,yy[,2]*0.8,col=col.arg[2],length=0.1) text(x,nx,col=col.arg[1],cex=cex.arg[1],adj=yadj.arg) par(oldpar) invisible(list(x=x,y=y)) } 22

23 mybiplot(zz,bb); dev.copy2eps(file="run0095-b12.eps") mybiplot(zz,bb,choices=3:2,scale=c(-1,1)); dev.copy2eps(file="run0095-b32.eps") PC Kagoshima Kochi Ktan Kfufu Yamaguchi Ehime Ooita Wakayama Miyazaki Nagasaki Hokkaido Shimane Tandoku Tokushima Kumamoto X65Sai Kagawa Hiroshima Rikon Tokyo Akita Okayama Kyoto Fukuoka Kaku Tottori Aomori Hyogo Osaka Iwate Nagano Yamanashi Mie Saga Nara Tomo Niigata Fukushima Ishikawa Gumma Toyama Kanagawa Miyagi amagata Fukui Gifu Shizuoka Chiba Konin Tochigi Ibaraki Aichi Ninzu Saitama ZoukaOkinawa Shiga PC Kfufu Ktan Kagoshima Yamaguchi Ehime Tandoku Ooita Wakayama Miyazaki Nagasaki Hokkaido Shimane Rikon Tokushima X65Sai Kaku Kumamoto Kagawa Okayama Hiroshima Akita Fukuoka Kyoto Hyogo Osaka Aomori Tottori Mie Yamanashi Nagano Iwate Saga Nara Gumma Fukushima Niigata Ishikawa Toyama Kanagawa Tomo Yamagata Miyagi GifuChiba Shizuoka Fukui Ibaraki Tochigi Aichi Okinawa Saitama Shiga Konin Ninzu Zouka Kochi Tokyo PC 1 run0095-b PC 3 run0095-b32 PC1 vs,,,, 65,, PC2?,,,, PC3?, 2.5 princomp # run0096.r # princomp # dat cat(" \n") a <- princomp(dat) print(summary(a)) biplot(a); dev.copy2eps(file="run0096-b1.eps") cat(" \n") a <- princomp(dat,cor=t) 23

24 print(summary(a)) biplot(a); dev.copy2eps(file="run0096-b2.eps") > source("run0096.r") Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation Proportion of Variance Cumulative Proportion Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation e-02 Proportion of Variance e-05 Cumulative Proportion e-01 Comp.10 Standard deviation e-02 Proportion of Variance e-06 Cumulative Proportion e+00 Importance of components: Comp.1 Comp.2 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.10 Standard deviation Proportion of Variance Cumulative Proportion

25 Comp Tomo X65Sai Saitama Nara Gifu Gumma Ibaraki Chiba Okinawa Shiga Kaku Shizuoka Tochigi Mie Aichi Toyama Hyogo Fukui Saga NaganoYamanashi Kagawa Wakayama Miyazaki Kanagawa Niigata Ninzu Zouka FukushimaOkayama Rikon Konin Ishikawa Hiroshima Tottori Aomori Osaka Yamagata Tokushima Kfufu Nagasaki Kumamoto Ehime Akita Yamaguchi Ktan Iwate Ooita Hokkaido Miyagi Shimane Fukuoka Kyoto Kagoshima Kochi Tandoku Tokyo Comp Kagoshima Kochi Ktan Kfufu Yamaguchi Ehime Ooita Wakayama Miyazaki Nagasaki Hokkaido Shimane Tandoku Tokushima Kumamoto X65Sai Kagawa Hiroshima Tokyo Akita Okayama Kyoto Fukuoka Rikon Kaku Tottori Aomori Hyogo Osaka Iwate Nagano Yamanashi Mie Saga Nara Tomo Niigata Fukushima Ishikawa Gumma Toyama Kanagawa Miyagi amagata Fukui Gifu Shizuoka Chiba Tochigi Konin Ibaraki Aichi Ninzu Saitama ZoukaOkinawa Shiga Comp.1 Comp.1 run0096-b1 run0096-b2 princomp() princomp(,cor=t) summary() 2.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", "Hiroshima", "Lotte", "Nichiham", "Seibu", "Kintetsu", "Orix", "Daiei","Taiyo") names(team) <- c(" ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " "," ") item <- c("daritsu","choudaritsu","shutsuruiritsu","shubiritsu", "Touruiboushiritsu","Shouritsu","Bougyoritsu") # na <- substr(item,1,nchar(item)-5) # "ritsu" names(na) <- item year <- 2000:2003 # par(mfrow=c(2,2),pty="s") # 2 x 2 x <- list(); pc <- list(); 25

26 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(year2000=c(-1,1),year2001=c(-1,-1), Year2002=c(-1,1),Year2003=c(-1,1)) # for(i in year) { j <- paste("year",i,sep=""); pc[[j]]$scores[,1:2] <- pc[[j]]$scores[,1:2] %*% diag(ex[[j]]); pc[[j]]$loadings[,1:2] <- pc[[j]]$loadings[,1:2] %*% diag(ex[[j]]); biplot(pc[[j]],main=j)} # dev.copy2eps(file="run0097.eps") # EPS par(mfrow=c(1,1),pty="s") 26

27 Year Year Comp Yakult Touruiboushi Kyojin Hansh Shubi Seibu Shou Yokohama Chunichi Daiei Chouda Da Hiroshima Shutsurui Orix Lotte hiham Kintetsu Bougyo Comp Chunichi anshin Yokohama Touruiboushi Yakult Shubi Daiei Seibu Kyojin Orix Shou Da Shutsu Chou Hiroshima Nichiham Lotte Bougyo Kintetsu Comp.1 Comp.1 Year Year Comp Orix okohama Chunichi Yakult Hanshin Touruiboushi Kyojin Lotte Nichiham Bougyo Hiroshima Shubi Daiei Kintetsu Seib Da Choud Shutsu Sho Comp Chunichi Shubi Hiroshima Yokohama Kyojin HanshinSho Nichiham Lotte Yakult Seibu Kintetsu Touruiboush Shutsur Da Daie Bougyo Chouda Orix Comp.1 Comp.1 run (Da) (Chouda) (Shutsurui) (Shubi) (Touruiboshi) (Shou) (Bougyo) 2000 = 2001 = 2002 = 2003 = ( 27

28 2.7 (SVD) (singular value decomposition) X = x x 1p x n1... x np } {{ } p n = x (1) x (n) = [x 1,..., x p ] X = UDV = d 1 u 1 v d p u p v p U = [u 1,..., u p ], V = [v 1,..., v p ] n p p p d 1 0 D =..., d 1 d p 0 0 d p X Σ = 1 n 1 X X = 1 n 1 V D2 V v 1,..., v p λ 1 = 1 n 1 d2 1,..., λ p = 1 n 1 d2 p y j = Xv j, j = 1,..., p Y = [y 1,..., y p ] = XV = UD Λ = 1 n 1 D2 Z = [z 1,..., z p ] = Y Λ 1/2 = n 1 U B = 1 n 1 X Z = 1 n 1 V D 28

29 # 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[1:5,1:5]\n") print(xx[1:5,1:5]) cat("(a$u %*% diag(a$d) %*% t(a$v))[1:5,1:5]\n") print((a$u %*% diag(a$d) %*% t(a$v))[1:5,1:5]) cat(" cumsum(a$d^2)/sum(a$d^2)\n") print(cumsum(a$d^2)/sum(a$d^2)) n <- nrow(xx) zz <- sqrt(n-1)*a$u # bb <- (1/sqrt(n-1))* a$v %*% diag(a$d) # mybiplot(zz,bb); dev.copy2eps(file="run0098-b12.eps") > source("run0098.r") [1] "d" "u" "v" d 10 u v xx[1:5,1:5] Zouka Ninzu Kaku Tomo Tandoku Hokkaido Aomori Iwate Miyagi Akita (a$u %*% diag(a$d) %*% t(a$v))[1:5,1:5] Zouka Ninzu Kaku Tomo Tandoku Hokkaido Aomori

30 Iwate Miyagi Akita cumsum(a$d^2)/sum(a$d^2) [1] [8] Kagoshima PC Kochi Ktan Kfufu Yamaguchi Ehime Ooita Wakayama Miyazaki Nagasaki Hokkaido Shimane Tandoku Tokushima Kumamoto X65Sai Kagawa Hiroshima Rikon Tokyo Akita Okayama Kyoto Fukuoka Kaku Tottori Aomori Hyogo Osaka Iwate Nagano Yamanashi Mie Saga Nara TomoNiigataFukushima Ishikawa Gumma Toyama Miyagi Kanagawa amagata Fukui Gifu Shizuoka Chiba Tochigi Konin Ibaraki Aichi Ninzu Saitama ZoukaOkinawa Shiga PC 1 run0098-b lambda: 1 n 1 X X z: Z b: B mybiplot (run0095.r)

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

$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$ Copyright (c) 004,005 Hidetoshi Shimodaira.. 3.. x... xp X = xn... xnp {{ p n x () = = [x,..., xp] x (n) x (i) xj X X n n nx 005-0-9 09:43:33 shimo R dat