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 <- 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" # 0.07957447.7968085 57.597873 34.633949 4.8889367 36.8663898 Rikon 8.4604553 6.808506 5.6378734.8440455 0.0347 0.0486787 0.38998474 38.09756568 5.606663 38.53467 Rikon.8379547 3.4340969 0.98084 0.080895 # # 4.74353e-8.88974e-7 -.40375e-4-4.686558e-5-3.40534e-5 Rikon -3.35945e-5.946434e-5 -.644075e-5.88974e-7-9.94e-7 0.0347 0.0486787 0.38998474 38.09756568 5.606663 38.53467 Rikon.8379547 3.4340969 0.98084 0.080895. 0.6 0. 0.6 0. 0 5 0 5 0.4 0.6 0. 0.6 v v = v v =., vp p vj =. j= 0 5 5 5 5 0 5 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 0.4 0.6 0. 0.6 0 5 5 5 5.0 0.5 Rikon.0 0.5 x (i) yiv x (i) yiv i =,..., n v = n x (i) yiv i= yn 0.6 0.4 0. 0.0 0. 0.4 Fukui Niigata Gifu Fukushima Iwate Nagano Shizuoka Gumma Mie Yamanashi Ishikawa WakayamaOkayama Saitama Kagawa Hyogo Hokkaido Miyazaki Fukuoka 0. 0.0 0. 0.4 0.6 run0087-s run0087-s0 Okinawa 3 0 5 0 5 0 Saitama Hyogo Hokkaido Fukuoka Okinawa Ishikawa Gifu Fukushima Iwate Shizuoka Mie Gumma Yamanashi Kagawa Miyazaki Okayama 0 5 0 5 0 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

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") [] 0 0 0 0 0 0 0 0 0 initial value 5484.690809 iter 0 value 49.04965 iter 0 value 47.466930 final value 47.4385 converged [] 0.03677-0.080609 0.3797844-0.6350376 0.764300-0.679779 [7] -0.03054 0.070694 0.04054889 0.05086 [,] Hokkaido.6588547 -.5054468 Iwate -8.60336604 3.7586794-3.4553368...... -.6897069 0.439747 Miyazaki.4330856 5.64047806 Okinawa 3.8537567 There were 50 or more warnings (use warnings() to see the first 50) 5 y 0 0 0 0 anagawa Okinawa Saitama Hokkaido Fukuoka Hyogo Miyazaki Gumma Wakayama Okayama Shizuoka Yamanashi Mie Kagawa Ishikawa Gifu 0 5 0 5 0 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 0 0 0 0 0 0 0 0 > source("run0089.r") initial value -3.69053 iter 0 value -3997.33688 iter 0 value -406.9440 final value -406.94956 converged [] 0.03675-0.080609 0.3797844-0.6350376 0.764300-0.679779 [7] -0.03054 0.070694 0.04054890 0.05086 [,] Hokkaido.6588547 -.5054468 Iwate -8.60336604 3.7586794-3.4553368...... -.6897069 0.439747 Miyazaki.4330856 5.64047806 Okinawa 3.8537566 There were 50 or more warnings (use warnings() to see the first 50) y run0089-s y: y:.4 0000 0 σ 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

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") # # 9.448707e-8 7.086530e-7-5.333795e-5-7.57e-6-8.88635e-6 Rikon -5.433006e-6.56678e-5-9.838466e-6.8346e-7 -.645638e-6 Rikon initial value -46.000000 iter 0 value -33.76884 final value -33.77458 converged [] 0.95696345-0.306448 0.3699859-0.40483635 0.985463 [6] -0.43040588-0.5958345 0.004907643 0.35036533 0.38005 [,] Hokkaido.786454587-0.69447996 Iwate -.33003459 0.58467509-3.70596460...... -0.843679-0.637708 Miyazaki 0.5880648 0.563643 Okinawa 4.56887557 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 y3 4 0 4 Nagano Fukushima Gifu Iwate Fukui Niigata Miyazaki Shizuoka Kagawa Okayama Gumma Wakayama Yamanashi Mie Ishikawa 0 0 0 0 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") 0.0347063-0.00053006 0.4004043-0.4600704 0.0337843-0.00053006 0.048678734-0.490355.086604-0.77806434 0.400404553-0.490354764 0.3899847-9.843384.764407-0.460070355.086604070-9.843384 38.0975657-6.64777044 0.03378430-0.7780643386.7644-6.6477704 5.606663 [] 0.03694-0.0795377 0.3799395-0.630407 0.7540588-0.6839763 [7] -0.0300084 0.069003 0.0403757 0.058684 [,] Hokkaido.65835480 -.5070706 Iwate -8.60697 3.835-3.457577...... -.6849847 0.43037858 Miyazaki.47853 5.6335558 Okinawa 3.8533693.000000000-0.00567593 0.494957-0.43988 0.0474959-0.00567593.000000000-0.49878 0.797588-0.8955544 0.49495698-0.4987775.0000000-0.7897 0.549397-0.43988084 0.7975877-0.7897.0000000-0.6864788 0.04749586-0.89555440 0.54930-0.686479.00000000 [] 0.957673-0.306346 0.37007066-0.40490395 0.9784563 [6] -0.43035736-0.609944 0.0048948 0.3503554 0.379936774 [,] Hokkaido.7860464-0.6940865 Iwate -.3994564 0.584889550-3.706073...... -0.84440839-0.663670 Miyazaki 0.58766833 0.55797839 Okinawa 4.5744335

0 0 0 0 V = (v,..., vp) Y = (y,..., y p ) y4 0 0 0 0 y5 4 0 4 Y = XV V V = Ip V p x, x,..., xp p y, y,..., y p v v v v, v v3 v,..., vr vr 4 0 4 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 [] 5.08005 3.4444357 0.9743003 0.96066 0.837380 0.0549588 [7] 0.069483 0.040804337 0.00468594 0.00993747 $vectors [,] [,] [,3] [,4] [,5] [,6] [,] 0.957673-0.36968439 0.957658-0.04636798-0.4594746-0.443088 [,] -0.306346-0.3506560 0.773744 0.5344748-0.09954994-0.758407 [3,] 0.37007066 0.085896 0.676335-0.604907 0.3638 0.05754 [4,] -0.40490395-0.479065-0.04845855-0.678763-0.490574 0.5048437 [5,] 0.9784563 0.5505-0.5957756-0.08998077-0.089455-0.55038 [6,] -0.43035736 0.30559 0.07340 0.435679-0.04940-0.097970 [7,] -0.609944 0.5008804 0.3459-0.3499549-0.36304 0.84378 [8,] 0.0048948 0.5395684 0.00684 0.90570-0.483309-0.454770 [9,] 0.3503554-0.86706-0.6433590-0.4944893-0.4557700 0.60 [0,] 0.379936774 0.44848 0.6775 0.780537-0.64483 0.40888 [,7] [,8] [,9] [,0] [,] 0.3704683-0.35547558-0.53773-0.05668 [,] -0.88778-0.008068 0.3505999 0.67868353 [3,] 0.087373 0.30075804-0.375333 0.39834356 [4,] 0.4660554 0.8566073-0.04848 0.03539 [5,] 0.789638-0.0690537-0.9488 0.57303579 [6,] -0.306888-0.5575865-0.759700-0.00368 [7,] -0.06373-0.6874668 0.508848 0.5644046 [8,] -0.0707873 0.5030490 0.38954-0.483809 [9,] -0.6599845 0.070973 0.093357 0.0570439 [0,] 0.0659877-0.3397676 0.0990653 0.05469995 Hokkaido.78605.8946503-0.38884-0.694086-0.05387806-0.357883 Iwate -.399456-0.3095053 -.07090 0.5848895 -.4790939 -.78795-3.70607 0.6534589-0.37304...... [,] [,] [,3] -0.844408.054473 0.38475-0.66367.40609 0.644958 Miyazaki 0.5876683.8787.484 0.557978 4.70755 0.5059443 Okinawa 4.57443 -.458857.5793804 [] 0.50800 0.834454 0.996597 0.95988 0.977654 0.988034 0.99557 [8] 0.9993 0.9997006.0000000 0 4 Miyazaki Wakayama Hokkaido Kagawa Okayama Fukuoka Hyogo Iwate NaganoYamanashi Mie Niigata FukushimaIshikawa Gumma Fukui Gifu Shizuoka Saitama Okinawa 4 0 4 run009-s run009-s3 0.93. 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 4 3 0 3 (j=,,3) yij = x (i) vj, i =,..., n, j =,..., r [,] [,] [,3] 5 6

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.78605.8946503-0.38884-0.694086-0.05387806-0.357883 Iwate -.399456-0.3095053 -.07090 0.5848895 -.4790939 -.78795-3.70607 0.6534589-0.37304 [] 0.50800 0.834454 0.996597 0.95988 0.977654 0.988034 0.99557 [8] 0.9993 0.9997006.0000000 Z (i=:5, j=:3) 9 [,] [,] [,3] Hokkaido.358886.056708-0.70735-0.306703-0.099097-0.369087 Iwate -.033548-0.78853 -.086574 0.59454-0.840865 -.7533860 -.643956 0.368997-0.378380 B (j=:3) [,] [,] [,3] 0.666757-0.665689 0.93093-0.778903-0.63564 0.46836 0.7463899 0.467895 0.6697-0.97849-0.65043-0.04777883 0.66003344 0.454-0.5873737-0.9536675 0.359896 0.056684-0.67658 0.9040993 0.96570 0.0070 0.959684 0.085900 0.7898009-0.507977-0.603078 Rikon 0.8565034 0.457 0.03808 > xx[:5,:5] Hokkaido -0.97998 -.70785546 0.76497 -.30683.809468-0.5530447 0.864054-0.677646-0.040046-0.047404 Iwate -0.8307487 0.5583559 -.45070 0.6783979-0.06030 0.557775 0.044658 -.736808-0.44605438 0.936733 -.8304833 0.9094973 -.504387 0.79658969-0.935397 > (zz %*% t(bb))[:5,:5] Hokkaido -0.97998 -.70785546 0.76497 -.30683.809468-0.5530447 0.864054-0.677646-0.040046-0.047404 Iwate -0.8307487 0.5583559 -.45070 0.6783979-0.06030 0.557775 0.044658 -.736808-0.44605438 0.936733 -.8304833 0.9094973 -.504387 0.79658969-0.935397 > (zz[,:3] %*% t(bb[,:3]))[:5,:5] Hokkaido 0.07979833 -.6076974 0.887948 -.396986.470885-0.546645 0.668536-0.4638890 0.3057-0.0086345 Iwate -0.784354 0.650946 -.4830853.040857-0.679 0.38354-0.0446957 -.096395 0.0646603 0.8803336 -.407007 0.8776076 -.37086.44666-0.6980683 0

0 z,z3: Z b,b3: B zz %*% t(bb) 0 3 4 0.5 0.0 0.5.0 0 Miyazaki Wakayama Hokkaido Kagawa Okayama Fukuoka Hyogo Iwate NaganoYamanashi Mie Niigata FukushimaIshikawa Gumma Fukui Gifu Shizuoka Saitama Okinawa run0093-z Rikon.0 0.5 0.0 0.5 run0093-b 0 3 4 xx run0093-zzbb zz[, :3] %*% t(bb[, :3]) 0 0.5 0.0 0.5.0 3 0 3 4 Miyazaki Wakayama Hokkaido Kagawa Fukuoka Okayama Iwate 4 3 0 andoku 3 Hyogo Yamanashi Nagano Mie Ishikawa Fukushima Niigata Gumma Fukui Shizuoka Gifu Saitama Okinawa run0093-z3 Rikon 0.6 0.4 0. 0.0 0. 0.4 0.6 3 run0093-b3 0 3 4 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 0 0.5 0 0.5 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 0.5 0 0.5 PC 0 3 4 0.5 0 0.5 Wakayama Miyazaki Hokkaido Rikon Kagawa Okayama Fukuoka Hyogo Mie Yamanashi Nagano Iwate Gumma Fukushima Niigata Ishikawa Gifu Shizuoka Fukui Okinawa Saitama 0 3 4 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) 0.5 0 0.5 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 9.44846 3.93803 3.58603.6736088 0.48040903 Proportion of Variance 0.7390 0.375 0.09856 0.03986 0.00976405 Cumulative Proportion 0.7390 0.863077 0.97996 0.996958 0.9988954 Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation 0.540695 0.09430 0.3493693 6.3377e-0 Proportion of Variance 0.00054066 0.0003809477 0.0005034 3.439684e-05 Cumulative Proportion 0.999433705 0.9998338 0.999963643 9.999980e-0 Comp.0 Standard deviation.538e-0 Proportion of Variance.96080e-06 Cumulative Proportion.000000e+00 Importance of components: Comp. Comp. Comp.3 Comp.4 Comp.5 Standard deviation.5433.8006789 0.985973 0.544653 0.486305 Proportion of Variance 0.5080 0.34444 0.09743 0.09606 0.083738 Cumulative Proportion 0.5080 0.834454 0.996597 0.95988 0.9776549 Comp.6 Comp.7 Comp.8 Comp.9 Standard deviation 0.347963 0.6359483 0.0000833 0.0684484764 Proportion of Variance 0.005496 0.0069483 0.004080434 0.000468594 Cumulative Proportion 0.9880345 0.995567 0.999306 0.999700653 Comp.0 Standard deviation 0.05475456 Proportion of Variance 0.000993747 Cumulative Proportion.0000000000 3 4

0.4 0. 0.0 0. Comp. 0.4 0. 0.0 0. 40 30 0 0 0 0 0 5 0 5 Saitama 40 30 0 0 0 0 0 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 0.3 0. 0. 0.0 0. 0. 0.3 0.4 Comp. Comp. run0096-b run0096-b princomp() princomp(,cor=t) Comp. 0.3 0. 0. 0.0 0. 0. 0.3 0.4 5 0 5 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. 0.4 0.0 0.4 Year000 0 4 Yakult Touruiboushi Kyojin Hansh Shubi Seibu Shou Yokohama Chunichi Daiei Chouda Da Shutsurui Orix Lotte hiham Kintetsu Bougyo 0.4 0.0 0.4 Comp. 0 4 Comp. 0.6 0. 0. 0.6 Year00 3 0 3 anshin Chunichi Yokohama Touruiboushi Yakult Shubi Seibu Orix 0.6 0. 0. 0.6 Comp. Daiei Kyojin Shou Da Shutsu Chou Nichiham Lotte Bougyo Kintetsu 3 0 3.7 (SVD) (singular value decomposition) x... xp X = xn... xnp {{ p X = UDV n x () = = [x,..., xp] x (n) = duv + + dpupv p Comp. 0.6 0. 0. Year00 4 0 Orix okohama Nichiham Lotte Bougyo Chunichi Yakult Hanshin Touruiboushi Kyojin Shubi Daiei Kintetsu 0.6 0. 0. Seib Da Choud Shutsu Sho 4 0 Comp. 0.6 0. 0. 0.6 Year003 4 0 3 Chunichi Shubi Yokohama Kyojin HanshinSho Nichiham Lotte Yakult Seibu Kintetsu Touruiboush Shutsur Da Daie Bougyo Chouda Orix 0.6 0. 0. 0.6 4 0 3 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. run0097 000 003 (Da) (Chouda) (Shutsurui) (Shubi) (Touruiboshi) (Shou) (Bougyo) 000 =00 =00 = 003 = (http://kamakura.cool.ne.jp/kojikiro/index.htm) 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

# 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 -0.97998 -.70785546 0.76497 -.30683.809468-0.5530447 0.864054-0.677646-0.040046-0.047404 Iwate -0.8307487 0.5583559 -.45070 0.6783979-0.06030 0.557775 0.044658 -.736808-0.44605438 0.936733 -.8304833 0.9094973 -.504387 0.79658969-0.935397 (a$u %*% diag(a$d) %*% t(a$v))[:5,:5] Hokkaido -0.97998 -.70785546 0.76497 -.30683.809468-0.5530447 0.864054-0.677646-0.040046-0.047404 9 Iwate -0.8307487 0.5583559 -.45070 0.6783979-0.06030 0.557775 0.044658 -.736808-0.44605438 0.936733 -.8304833 0.9094973 -.504387 0.79658969-0.935397 cumsum(a$d^)/sum(a$d^) [] 0.50800 0.834454 0.996597 0.95988 0.977654 0.988034 0.99557 [8] 0.9993 0.9997006.0000000 PC 0 0.5 0 0.5 3 Wakayama Miyazaki Hokkaido Kagawa Rikon Okayama Fukuoka Iwate Nagano Yamanashi Mie NiigataFukushima Ishikawa Gumma Fukui Gifu Shizuoka 3. 9-0 PC Hyogo run0098-b Saitama Okinawa 0.5 0 0.5 lambda: n X X z: Z b: B 3. 9- mybiplot (run0095.r) 3.3 9-3 30