数量化４類（多変量解析入門２、河口至商先生）

多変量解析入門２、河口至商著の数量化４類をプロットしてみました。

```{r}

library(dplyr)

#数量化４類

e=array(0,dim=c(14,14))

e[1,]=c(0,0.86,0.42,0.42,0.18,0.06,0.07,0.04,0.02,0.07,0.09,0.12,0.13,0.16)

e[2,]=c(0,0,0.5,0.44,0.22,0.09,0.07,0.07,0.02,0.04,0.07,0.11,0.13,0.14)
e[3,]=c(0,0,0,0.81,0.47,0.17,0.1,0.08,0.02,0.01,0.02,0.01,0.05,0.03)
e[4,]=c(0,0,0,0,0.54,0.25,0.1,0.09,0.02,0.01,0,0.01,0.02,0.04)
e[5,]=c(0,0,0,0,0,0.61,0.31,0.26,0.07,0.02,0.02,0.01,0.02,0)
e[6,]=c(0,0,0,0,0,0,0.62,0.45,0.14,0.08,0.02,0.02,0.02,0.01)
e[7,]=c(0,0,0,0,0,0,0,0.73,0.22,0.14,0.05,0.02,0.02,0)
e[8,]=c(0,0,0,0,0,0,0,0,0.33,0.19,0.04,0.03,0.02,0.02)
e[9,]=c(0,0,0,0,0,0,0,0,0,0.58,0.37,0.27,0.2,0.23)
e[10,]=c(0,0,0,0,0,0,0,0,0,0,0.74,0.5,0.41,0.28)
e[11,]=c(0,0,0,0,0,0,0,0,0,0,0,0.76,0.62,0.55)
e[12,]=c(0,0,0,0,0,0,0,0,0,0,0,0,0.85,0.68)
e[13,]=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0.76)

e=e+t(e)

H=array(0,dim=c(nrow(e),ncol(e)))

for(j in 1:nrow(H)){
for(i in 1:ncol(H)){

if(j!=i){  

H[j,i]=e[i,j]+e[j,i]

}    
}  
}

alpha=-max(H)

H=H+alpha;diag(H)=0

for(j in 1:ncol(H)){

diag(H)[j]=-sum(H[j,])  

}

X=eigen(H)$vectors[,1];Y=eigen(H)$vectors[,2]

lambda=eigen(H)$values[1];mu=eigen(H)$values[2]

Q=lambda*sum(X^2)+mu*sum(Y^2)

print(t(X)%*%H%*%X+t(Y)%*%H%*%Y)

plot(-X,Y,type="p",xlab="X",ylab="Y",col=2)

lambda2_n1=eigen(H)$values[length(eigen(H)$values)-1]

lambda_hats=eigen(H)$values[1:(length(eigen(H)$values)-2)]

propotion=(lambda_hats-lambda2_n1)/(sum(lambda_hats)-(length(eigen(H)$values)-2)*lambda2_n1)

Z=cbind(X,Y)



```

```{r}

#k-mean method(multi dimension)

gap_stats=function(X){
  
 X=data.frame(X) 

 classes=sort(unique(X$c))
 
 D=c()
 
 for(k in classes){

 dat=X %>% filter(c==k)
 
 mat=as.matrix(dat[,colnames(dat) %in% c("x.X","x.Y")])
 
 D_k=c()
 
 for(l in 1:nrow(mat)){
   
 mat_one=mat[l,]
 
 diff=sum((mat-mat_one)^2)
 
 D_k=c(D_k,diff)
   
 }
 
 D=c(D,sum(D_k))
 
}

 n_data=X %>% group_by(c) %>% summarise(n=n())
 
 n_data=n_data[order(n_data$c,decreasing = F),]
 
 Wf=sum(D/(2*n_data$n))
 
 times=1000
   
W_values=c()
  
for(m in 1:times){
  
  W_vec=c()
  
for(K in classes){  
   
 D=c()
 
 dat=X %>% filter(c==K)
 
 dat_num=data.frame(num=c(1:nrow(dat)),dat)
 
 nums=matrix(sample(c(1:nrow(dat)),nrow(dat)*2,replace = T),ncol=2)
 
 nums_data=data.frame(num=nums[,1])
 
 dat_sub1=dat_num[,1:(2)]
 
 colnames(dat_sub1)=c("num","x.X")
 
 dat_sub1=left_join(nums_data,dat_sub1,"num")
 
 nums_data=data.frame(num=nums[,2])
 
 dat_sub2=dat_num[,c(1,3)]
 
 colnames(dat_sub2)=c("num","x.Y")
 
 dat_sub2=left_join(nums_data,dat_sub2,"num")
 
 
 dats=cbind(dat_sub1,dat_sub2,"num")
 
 mat=as.matrix(dats[,colnames(dats) %in% c("x.X","x.Y")])
 
 D_k=c()
 
 for(l in 1:nrow(mat)){
   
 mat_one=mat[l,]
 
 diff=sum((mat-mat_one)^2)
 
 D_k=c(D_k,diff)
   
 }
 
 W=sum(D_k/(2*n_data$n[n_data$c==K]))
 
 W_vec=c(W_vec,W)
   
 }
  
W_val=sum(W_vec)  

W_values=c(W_values,W_val)
  
}
    
sd_k=sum(log(W_values)-mean(log(W_values)))

sd_k[abs(sd_k)<10^(-4)]=0
 
 gap=c(mean(log(W_values)),log(Wf),ifelse(sd_k>0,sqrt(sd_k),0))
 
 return(gap)
 
}


```

```{r}
#k-mean method(multi dimension)

library(dplyr)

n=nrow(Z);dim=ncol(Z)

x=Z*100

cluster_num=4

C=x[sample(1:n,cluster_num),]

Data=data.frame(num=rep(0,cluster_num))

for(j in 1:dim){

  C_D=data.frame(C=C[,j])

  colnames(C_D)=c(as.character(paste0("c.",j)))

  Data=cbind(Data,C_D)

}

Data=Data %>% mutate(L=0)


ite=10

for(k in 1:ite){

 X=data.frame(x=x,c=0,L=0)

for(j in 1:nrow(X)){

  x_one=x[j,]

  c_x=data.frame(num=1:nrow(C),c=C,L=0) 

  for(l in 1:nrow(C)){

  c_x$L[l]=sum((x_one-C[l,])^2)

  }

  X$c[j]=c_x$num[c_x$L==min(c_x$L)]

  X$L[j]=min(c_x$L)

}


XX=data.frame(c=1:cluster_num)

for(j in 1:dim){

  x_c=data.frame(c=X$c,x=X[,j])

  x_c=x_c %>% group_by(c) %>% summarise(center_of_gravity=mean(x))

  colnames(x_c)=c("c",as.character(paste0("c.",j)))

  XX=left_join(XX,x_c,by="c")

} 

x_c=X %>% group_by(c) %>% summarise(L=sum(L))

XX=left_join(XX,x_c,by="c")

C=XX[,2:(1+ncol(C))]

data_c=data.frame(num=1:cluster_num)

for(j in 1:dim){

  C_data=data.frame(num=1:cluster_num,c=C[,j])

  colnames(C_data)=c("num",as.character(paste0("c.",j)))

  data_c=left_join(data_c,C_data,by="num")

}

data=data_c %>% mutate(L=sum(XX$L))

data$num=k

Data=rbind(Data,data)

print(paste0("C=",C))

print(paste0("L=",sum(XX$L)))

}

gap_stats(X)


```