library(ggplot2)
n<- 300
x0 <- (n/5):n
x<- x0/n*2 + cos(2*3.14 * x0/n * 3) + rnorm(length(x0),0,0.1)
y<- sin(2*3.14 * x0/n * 3) + rnorm(length(x0),0,0.1)
z <- x0/n
options(rgl.useNULL = TRUE)
options(rgl.printRglwidget = TRUE)
require(rgl)
## Loading required package: rgl
grid1<- seq(-1,3,0.1)
grid2<- seq(-1.5,1.5,0.1)
vis_x = rep(grid1,length(grid2))
vis_y = rep(grid2,each=length(grid1))
vis_z <- rep(0,length(vis_y))
for(ii in 1:length(x0)){
tmp1<- which.min( abs(grid1-x[ii]) )
tmp2<- which.min( abs(grid2-y[ii]) )
vis_z[(tmp2-1) * length(grid1) + tmp1] <- z[ii]
}
df<-data.frame(x=vis_x ,y=vis_y, z = vis_z,color = 1)
with(df, {
plot3d(
x = x, y = y, z = z, type = "s",col=2,size=0.5
)
})
#######################
n0<-50
x<-c(n0:1 ,rep(0,n0),rep(0,n0),rep(0,n0),
1:n0, rep(n0,n0)) +10 + rnorm(4*n0,0,1)
## Warning in c(n0:1, rep(0, n0), rep(0, n0), rep(0, n0), 1:n0, rep(n0, n0)) + :
## longer object length is not a multiple of shorter object length
y<-c(rep(0,n0),1:n0 ,n0:1 ,rep(0,n0),
1:n0, n0:1) +10 + rnorm(4*n0,0,1)
## Warning in c(rep(0, n0), 1:n0, n0:1, rep(0, n0), 1:n0, n0:1) + 10 + rnorm(4 * :
## longer object length is not a multiple of shorter object length
z<-c(rep(0,n0),rep(0,n0),1:n0 ,n0:1,
1:n0, n0:1) + 10+ rnorm(4*n0,0,1)
## Warning in c(rep(0, n0), rep(0, n0), 1:n0, n0:1, 1:n0, n0:1) + 10 + rnorm(4 * :
## longer object length is not a multiple of shorter object length
options(rgl.useNULL = TRUE)
options(rgl.printRglwidget = TRUE)
require(rgl)
grid1<- seq(0,n0,1)
grid2<- seq(0,n0,1)
vis_x = rep(grid1,length(grid2))
vis_y = rep(grid2,each=length(grid1))
vis_z <- rep(0,length(vis_y))
for(ii in 1:length(x0)){
tmp1<- which.min( abs(grid1-x[ii]) )
tmp2<- which.min( abs(grid2-y[ii]) )
vis_z[(tmp2-1) * length(grid1) + tmp1] <- z[ii]
}
df<-data.frame(x=vis_x ,y=vis_y, z = vis_z,color = 1)
with(df, {
plot3d(
x = x, y = y, z = z, type = "s",col=2,size=0.5
)
})
library(plot3D)
scatter3D(x, y, z, type="b")
x<-cbind(x,y,z)
n<-dim(x)[1]
MDS
### MDS
###########
MDS<-function(x){
n<- dim(x)[1]
q<-2
H<-diag(1,n)-rep(1,n)%*%t(rep(1,n))/n
# B<-t(xx)%*%H
# B<-t(B)%*%B
cross<-x%*%t(x)
num<-length(cross[,1])
rowm<-matrix(rep(diag(cross),num),num,num,byrow = TRUE)
colm<-matrix(rep(diag(cross),num),num,num,byrow = FALSE)
TT<-rowm+colm-2*cross
# TT<-matrix(0,n,n)
# for(i in 1:n) {
# for(j in 1:n){
# TT[i,j]<-sum(abs((x[i,]-x[j,])^3))^(1/3)
# }
# }
B<-t(H)%*%TT%*%H
svd1<-svd(B)
par(mfrow=c(1,1))
xx<-svd1$u[,1:2]%*%diag(sqrt(svd1$d[1:2]))
#plot(xx,main = "MDS")
xx <- as.data.frame(xx)
colnames(xx) <- c("x","y")
if(n<250){
p<- ggplot(xx, aes(x=x, y=y)) +
geom_point(aes(color= 1:n - n/2 >0, alpha=abs(1:n - n/2))) +
scale_color_manual(values=c("blue","red")) + theme_bw() + theme(legend.position="none") +
annotate("text",x=xx[25,1],y=xx[25,2],label="1-50")+
annotate("text",x=xx[50*1+25,1],y=xx[50*1+25,2],label="51-100")+
annotate("text",x=xx[50*2+25,1],y=xx[50*2+25,2],label="101-150")+
annotate("text",x=xx[50*3+25,1],y=xx[50*3+25,2],label="151-200")+
annotate("text",x=xx[50*4+25,1],y=xx[50*4+25,2],label="201-250")
}else{
p<- ggplot(xx, aes(x=x, y=y)) +
geom_point(aes(color= 1:n - n/2 >0, alpha=abs(1:n - n/2))) +
scale_color_manual(values=c("blue","red")) + theme_bw() + theme(legend.position="none") +
annotate("text",x=xx[25,1],y=xx[25,2],label="1-50")+
annotate("text",x=xx[50*1+25,1],y=xx[50*1+25,2],label="51-100")+
annotate("text",x=xx[50*2+25,1],y=xx[50*2+25,2],label="101-150")+
annotate("text",x=xx[50*3+25,1],y=xx[50*3+25,2],label="151-200")+
annotate("text",x=xx[50*4+25,1],y=xx[50*4+25,2],label="201-250")+
annotate("text",x=xx[50*4+25,1],y=xx[50*4+25,2],label="251-300")+
annotate("text",x=xx[50*5+25,1],y=xx[50*5+25,2],label="301-350")
}
p
}
result12<- MDS(x)
result1<- MDS(x[1:200 ,])
result2<- MDS(x[201:300,])
result12
result1
## Warning: Removed 1 rows containing missing values (`geom_text()`).
Local Tagent Spave Alignment (Local Linear Estimation)
## Warning: Removed 1 rows containing missing values (`geom_text()`).
## Warning: Removed 1 rows containing missing values (`geom_text()`).
