非線形カーネル主成分

#kernel component analysis


library(mlbench)

library(kernlab)

n=200

k=10

dat=mlbench.spirals(n,cycles=1.2,sd=0.16)

X=dat$x;Y=dat$classes

data=data.frame(x1=X[,1],x2=X[,2],class=Y)

plot(data$x1[data$class==1],data$x2[data$class==1],type="p",col=2,xlab="x1",ylab="x2",xlim=c(min(data$x1),max(data$x1)),ylim=c(min(data$x2),max(data$x2)))

par(new=T)

plot(data$x1[data$class==2],data$x2[data$class==2],type="p",col=3,xlab="x1",ylab="x2",xlim=c(min(data$x1),max(data$x1)),ylim=c(min(data$x2),max(data$x2)))

kernel=function(x,y){

return(exp(-sum((x-y)^2)))  

}

K=array(0,dim=c(nrow(X),nrow(X)))

N=nrow(K)

val=0

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

val=val+kernel(X[i,],X[j,])   

}
}  

val=val/(N^2)

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

  value=kernel(X[j,],X[i,])+kernel(X,X[j,])/N+kernel(X,X[i,])/N+val

  K[j,i]=value

}
}

eigen_values=eigen(K)$values

eigen_vectors=eigen(K)$vectors

propotion=eigen_values/sum(eigen_values)

comp_principles=eigen_vectors

for(j in 1:length(eigen_values)){

comp_principles[,j]=sqrt(eigen_values[j])*comp_principles[,j]  

}

comp_principles=comp_principles[,1:2]

plot(comp_principles[,1],comp_principles[,2])

comp_principles=data.frame(comp_principles) %>% mutate(class=kmeans(comp_principles,2)$cluster,y=Y) 

plot(comp_principles$X1[comp_principles$class==1],comp_principles$X2[comp_principles$class==1],type="p",col=1,xlab="x1",ylab="x2",xlim=c(min(comp_principles$X1),max(comp_principles$X1)),ylim=c(min(comp_principles$X2),max(comp_principles$X2)))

par(new=T)

plot(comp_principles$X1[comp_principles$class==2],comp_principles$X2[comp_principles$class==2],type="p",col=2,xlab="x1",ylab="x2",xlim=c(min(comp_principles$X1),max(comp_principles$X1)),ylim=c(min(comp_principles$X2),max(comp_principles$X2)))