GeoRatio_graph: Building GeoRatio Graph
In Yanhao29/GeoRatio: Nonlinear dimension reduction and clustering

Description Usage Arguments Value Examples

Takes in data, minimum nbhd size, maximum nbhd size, intrinsic dimension, variance threshold and construct a graph for estimating GeoRatio measure

1 2	GeoRatio_graph(X, k_max = k_M, k_min = k_m, d, thresh = 0.9, gmode = "max", distance = NaN)

`X`	Data.
`k_max`	maximum neighborhood size
`k_min`	minimum neighborhood size
`d`	intrinsic dimension
`thresh`	threshold of proportion of variance explained by d leading eigenvalues in each neighborhood
`gmode`	graph edge choosing
`distance`	if data passed in is a distance matrix

a list of GeoRatio dissimilarity, geodesic distance, nbhd adjacency matrix, weighted graph, nbhd size

############################## example I: V-shape
## package for pam() kmeans clustering
library(cluster)
library(igraph) ## used by GeoRatio_graph

## load data
data(V_shape)

## intrinsic dimension of the data
trueDim = 1

############## correlation dimension estimation
## number of epsilons
num = 50

## pair-wise Euclidean distance
dis = dist(V_shape,method = "euclidean", diag = FALSE, upper = TRUE)

## correlation dimension estimation
est = correDimEst(num,dis)

## plot correlation dimension estimation
par(mfrow=c(1,1))
plot(est$x[2:(length(est$x)-1)],est$deri,type="l",ylim=c(0,10),xlab="log(epsilon)",ylab="est. correlation dimension",
     main="Estimated correlation dimension")
abline(trueDim,0,lty=3,col="gray60",lwd=2)

##### decide k_M (plateau log(epsi) from -4.8 to -2.8)
epsi_s = -4.5
epsi_e = -2.5
epsi_seq = seq(epsi_s,epsi_e,length=5)

dis_L2 = as.matrix(dist(V_shape))

for(i in 1:length(epsi_seq)){
  epsi = epsi_seq[i]
  dist_thre = exp(epsi)
  epsi_nbhd = apply(dis_L2,1,function(x) sum(x<dist_thre))
  avg_epsi_nbhd = mean(epsi_nbhd)
  print(epsi)
  print(avg_epsi_nbhd)
}

## parameters for GeoRatio_graph()
K = 2
k_M = 15
k_m = 3
pca_thre = 0.9

set.seed(100)

V_GR = GeoRatio_graph(V_shape,k_max=k_M,k_min=k_m,trueDim,thresh=pca_thre,gmode="max")

dis_L2 = dist(V_shape)
dis_G = as.dist(V_GR$dis_G)
dis_GR = as.dist(V_GR$dis_GR)

hc_L2 = hclust(dis_L2,method="average")
hc_G = hclust(dis_G,method="average")
hc_GR = hclust(dis_GR,method="average")

indi_KML2 = pam(dis_L2,K,diss=T)$clustering
indi_KMG = pam(dis_G,K,diss=T)$clustering
indi_KMGR = pam(dis_GR,K,diss=T)$clustering
indi_HL2 = cutree(hc_L2,k=K)
indi_HG = cutree(hc_G,k=K)
indi_HGR = cutree(hc_GR,k=K)

lpca_V = kplanes(V_shape,K,trueDim,iter.max=100,thresh=1e-5)
indi_lpca = lpca_V$indicator

par(mfrow=c(3,3))
plot(V_shape,col=indi_HGR,main="HC-GR")
plot(V_shape,col=indi_HG,main="HC-G")
plot(V_shape,col=indi_HL2,main="HC-L2")
plot(V_shape,col=indi_KMGR,main="KMC-GR")
plot(V_shape,col=indi_KMG,main="KMC-G")
plot(V_shape,col=indi_KML2,main="KMC-L2")
plot(V_shape,col=indi_lpca,main="LPCA")

############################## example II: Swiss roll
## package for 3d plot
library(rgl)
## package for pam() kmeans clustering
library(cluster)

library(igraph) ## used by GeoRatio_graph

## load data
data(SwissRoll)

## intrinsic dimension of the data
trueDim = 2

############## correlation dimension estimation
## number of epsilons
num = 50

## pair-wise Euclidean distance
dis = dist(SwissRoll,method = "euclidean", diag = FALSE, upper = TRUE)

## correlation dimension estimation
est = correDimEst(num,dis)

## plot correlation dimension estimation
par(mfrow=c(1,1))
plot(est$x[2:(length(est$x)-1)],est$deri,type="l",ylim=c(0,10),xlab="log(epsilon)",ylab="est. correlation dimension",
     main="Estimated correlation dimension")
abline(trueDim,0,lty=3,col="gray60",lwd=2)

##### decide k_M (plateau log(epsi) from -4.8 to -2.8)
epsi_s = -1.2
epsi_e = 1.8
epsi_seq = seq(epsi_s,epsi_e,length=5)

dis_L2 = as.matrix(dist(SwissRoll))

for(i in 1:length(epsi_seq)){
  epsi = epsi_seq[i]
  dist_thre = exp(epsi)
  epsi_nbhd = apply(dis_L2,1,function(x) sum(x<dist_thre))
  avg_epsi_nbhd = mean(epsi_nbhd)
  print(epsi)
  print(avg_epsi_nbhd)
}

## number of clusters
K = 8
k_M = 15
k_m = 4
pca_thre = 0.9

set.seed(100)

SwissRoll_GR = GeoRatio_graph(SwissRoll,k_max=k_M,k_min=k_m,trueDim,thresh=pca_thre,gmode="max")

dis_L2 = dist(SwissRoll)
dis_G = as.dist(SwissRoll_GR$dis_G)
dis_GR = as.dist(SwissRoll_GR$dis_GR)

hc_L2 = hclust(dis_L2,method="average")
hc_G = hclust(dis_G,method="average")
hc_GR = hclust(dis_GR,method="average")

indi_KML2 = pam(SwissRoll,K)$clustering
indi_KMG = pam(dis_G,K,diss=T)$clustering
indi_KMGR = pam(dis_GR,K,diss=T)$clustering
indi_HL2 = cutree(hc_L2,k=K)
indi_HG = cutree(hc_G,k=K)
indi_HGR = cutree(hc_GR,k=K)

lpca_V = kplanes(SwissRoll,K,trueDim,iter.max=100,thresh=1e-5)
indi_lpca = lpca_V$indicator

open3d()
plot3d(SwissRoll,col=indi_HGR,main="HC-GR")
open3d()
plot3d(SwissRoll,col=indi_HG,main="HC-G")
open3d()
plot3d(SwissRoll,col=indi_HL2,main="HC-L2")
open3d()
plot3d(SwissRoll,col=indi_KMGR,main="KMC-GR")
open3d()
plot3d(SwissRoll,col=indi_KMG,main="KMC-G")
open3d()
plot3d(SwissRoll,col=indi_KML2,main="KMC-L2")
open3d()
plot3d(SwissRoll,col=indi_lpca,main="LPCA")

############################## example III: Open box
## package for 3d plot
library(rgl)
## package for pam() kmeans clustering
library(cluster)
library(igraph) ## used by GeoRatio_graph

## load data
data(OpenBox)

## intrinsic dimension of the data
trueDim = 2

############## correlation dimension estimation
## number of epsilons
num = 50

## pair-wise Euclidean distance
dis = dist(OpenBox,method = "euclidean", diag = FALSE, upper = TRUE)

## correlation dimension estimation
est = correDimEst(num,dis)

## plot correlation dimension estimation
par(mfrow=c(1,1))
plot(est$x[2:(length(est$x)-1)],est$deri,type="l",ylim=c(0,10),xlab="log(epsilon)",ylab="est. correlation dimension",
     main="Estimated correlation dimension")
abline(trueDim,0,lty=3,col="gray60",lwd=2)

##### decide k_M (plateau log(epsi) from -4.8 to -2.8)
epsi_s = -3.5
epsi_e = 0
epsi_seq = seq(epsi_s,epsi_e,length=5)

dis_L2 = as.matrix(dist(OpenBox))

for(i in 1:length(epsi_seq)){
  epsi = epsi_seq[i]
  dist_thre = exp(epsi)
  epsi_nbhd = apply(dis_L2,1,function(x) sum(x<dist_thre))
  avg_epsi_nbhd = mean(epsi_nbhd)
  print(epsi)
  print(avg_epsi_nbhd)
}

## number of clusters
K = 6
k_M = 20
k_m = 4
pca_thre = 0.9

set.seed(100)

OpenBox_GR = GeoRatio_graph(OpenBox,k_max=k_M,k_min=k_m,trueDim,thresh=pca_thre,gmode="max")

dis_L2 = dist(OpenBox)
dis_G = as.dist(OpenBox_GR$dis_G)
dis_GR = as.dist(OpenBox_GR$dis_GR)

hc_L2 = hclust(dis_L2,method="average")
hc_G = hclust(dis_G,method="average")
hc_GR = hclust(dis_GR,method="average")

indi_KML2 = pam(dis_L2,K,diss=T)$clustering
indi_KMG = pam(dis_G,K,diss=T)$clustering
indi_KMGR = pam(dis_GR,K,diss=T)$clustering
indi_HL2 = cutree(hc_L2,k=K)
indi_HG = cutree(hc_G,k=K)
indi_HGR = cutree(hc_GR,k=K)

lpca_V = kplanes(OpenBox,K,trueDim,iter.max=100,thresh=1e-5)
indi_lpca = lpca_V$indicator

open3d()
plot3d(OpenBox,col=indi_HGR,main="HC-GR",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_HG,main="HC-G",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_HL2,main="HC-L2",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_KMGR,main="KMC-GR",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_KMG,main="KMC-G",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_KML2,main="KMC-L2",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))
open3d()
plot3d(OpenBox,col=indi_lpca,main="LPCA",xlim=c(0,2),ylim=c(0,2),zlim=c(0,2))