多変量グラフ解析入門　河口至商・水田正弘先生　森北出版

```r
#星座グラフ p.80

library(dplyr)

data=data.frame(num=1:20,x0=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

guzai_pi=array(0,dim=c(nrow(data),ncol(data)-1))

max_x=apply(data[,2:ncol(data)],2,max)

min_x=apply(data[,2:ncol(data)],2,min)

X=data[,2:ncol(data)]


for(i in 1:ncol(guzai_pi)){
for(j in 1:nrow(guzai_pi)){
  
guzai_pi[j,i]=(X[j,i]-min_x[i])/(max_x[i]-min_x[i])  
  
}  
}



guzai_pi=pi*guzai_pi

guzai=guzai_pi

w=1/ncol(guzai_pi)

sita=c();V_i=c()

for(j in 1:nrow(guzai)){
  
sita=c(sita,atan(sum(w*sin(guzai[j,]))/sum(w*cos(guzai[j,]))) ) 

V_i=c(V_i,1-sqrt(sum(w*cos(guzai[j,]))^2+sum(w*sin(guzai[j,]))^2))
  
}

V_c=1-sqrt(sum((1/nrow(guzai))*cos(sita))^2+sum((1/nrow(guzai))*sin(sita))^2)

plot(apply(w*cos(guzai),1,sum),apply(w*sin(guzai),1,sum),type="p",xlab="cos",ylab="sin",main="星座グラフ",xlim=c(-1,1),ylim=c(0,1))

data=data %>% dplyr::mutate(cos=apply(w*cos(guzai),1,sum),sin=apply(w*sin(guzai),1,sum),V_i=V_i)


#cluster(package,kmeans使用)

cluster=6

result=kmeans(cbind(data$cos,data$sin),cluster)

data=data %>% dplyr::mutate(cluster=result$cluster)

for(j in 1:cluster){

if(j!=cluster){
  
data_sub=data %>% dplyr::filter(cluster==j)

plot(data_sub$cos,data_sub$sin,type="p",xlim=c(min(data$cos),max(data$cos)),ylim=c(min(data$sin),max(data$sin)),col=j,xlab="cos",ylab="sin",main="星座グラフ")

par(new=T)
  
}else{
  
data_sub=data %>% dplyr::filter(cluster==j)
  
plot(data_sub$cos,data_sub$sin,type="p",xlim=c(min(data$cos),max(data$cos)),ylim=c(min(data$sin),max(data$sin)),col=j,xlab="cos",ylab="sin",main="星座グラフ")
  
  
} }

```

```r

#正準分解法　p.116

s1980=matrix(c(0.624,0.372,0.504,0.567,0.45,0.488,0.263,0.261,0.243,0.27,0.262,0.259,161,134,153,132,134,135,1832,1745,1671,1797,1746,1753,516,472,482,537,492,509,716,711,705,680,660,773,435,341,404,371,402,421,108,50,83,111,87,93,3.37,4.43,2.95,3.17,3.73,4.18),nrow=6)

s1981=matrix(c(0.554,0.472,0.603,0.491,0.536,0.344,0.274,0.268,0.268,0.255,0.272,0.252,181,151,135,120,114,105,1932,1840,1800,1672,1734,1622,582,522,492,492,503,422,647,753,751,710,608,766,390,370,370,349,379,399,54,62,119,105,71,70,3.66,3.71,2.88,4.3,3.32,4.41),nrow=6)

s1982=matrix(c(0.504,0.577,0.569,0.375,0.533,0.449,0.254,0.266,0.254,0.24,0.262,0.250,139,143,133,95,118,125,1711,1760,1720,1497,1701,1681,499,50,489,418,494,451,704,667,714,687,618,801,393,344,374,359,368,398,99,63,112,87,103,39,3.33,3.27,2.93,3.64,3.44,3.92),nrow=6)

s1983=matrix(c(0.542,0.439,0.59,0.434,0.496,.5,0.269,0.263,0.275,0.266,0.274,0.272,164,160,156,149,169,137,1910,1844,1922,1826,1899,1848,594,548,591,531,573,539,725,706,684,619,600,679,433,453,442,326,428,453,165,63,134,90,86,78,3.65,4.11,3.77,4.52,4.22,4.52),nrow=6)

s1984=matrix(c(0.625,0.598,0.554,0.418,0.434,0.374,0.274,0.282,0.268,0.264,0.264,0.263,167,191,186,101,165,100,1915,2007,1989,1648,1834,1620,614,615,618,498,544,459,717,660,717,695,662,682,443,442,468,428,437,428,140,91,94,49,71,110,3.37,3.82,3.66,4.76,4.46,4.55),nrow=6)

dat=data.frame(rbind(cbind(1980,s1980),cbind(1981,s1981),cbind(1982,s1982),cbind(1983,s1983),cbind(1984,s1984)))
dat=data.frame(rbind(cbind(1980,s1980),cbind(1981,s1981),cbind(1982,s1982),cbind(1983,s1983)))

colnames(dat)=c("year",paste0("V",c(1:ncol(s1980))))

#データの大小感調整(特に0~1だと２乗は小さくなるため)
dat[,2]=dat[,2]*100
dat[,3]=dat[,3]*100
dat[,4]=dat[,4]/10
dat[,5]=dat[,5]/100
dat[,6]=dat[,6]/10
dat[,7]=dat[,7]/10
dat[,8]=dat[,8]/10
dat[,10]=dat[,10]*10

t=2
n=nrow(s1980)
m=ncol(s1980)
years=sort(unique(dat$year))
Y=c(unlist(dat[,-1]))

f=function(x){
a=matrix(x[1:(t*n)],nrow=n);x=x[-c(1:(t*n))]
b=matrix(x[1:c(t*m)],nrow=m);x=x[-c(1:(t*m))]
c=matrix(x,nrow=length(years))

values=c()
for(k in 1:length(years)){
values=c(values,t(b%*%t(a*c(c[k,]))))  
}

return(values)  
}


h=0.01
lam=0.01
ite=3*10^3
sita=rep(1,t*(n+m+length(years)))
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}


plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,f(sita))),max(c(Y,f(sita)))),xlab="index",ylab="",main=paste0("重相関：",cor(Y,f(sita))),col=2)
par(new=T)
plot(c(1:length(Y)),f(sita),xlim=c(1,length(Y)),ylim=c(min(c(Y,f(sita))),max(c(Y,f(sita)))),xlab="index",ylab="",main=paste0("重相関：",cor(Y,f(sita))),col=3)
pred=f(sita)


avec=matrix(sita[1:(t*n)],nrow=n);sita=sita[-c(1:(t*n))]
bvec=matrix(sita[1:c(t*m)],nrow=m);sita=sita[-c(1:(t*m))]
cvec=matrix(sita,nrow=length(years))

plot(avec)
plot(bvec)
plot(cvec)
```