- 空間に多様体がある
- 多様体をあっちこっち観察すると、多様体の存在する点が観測される
- その複数の点の位置データから、多様体が空間上で何次元なのかを知るにはどうしたらよいのだろう?
- 「多次元最小全域木」を作ってやると、多様体の次元までは、「多次元最小全域木」の体積は大きくなっていくが、多様体の次元を超えると、体積が増えなくなる
- これを利用して、多様体の次元がわかりそうだ
- ちなみに、「増えなくなる」というのは、k次元最小全域木の体積とk+1次元最小全域木の体積の比が1になるということだから、k次元最小全域木とk+1次元最小全域木の体積比を「k+1次元化高さ」と呼ぶことにすれば、「k+1次元化高さ」が1になったら、多様体の次元はk次元ということになるのだろう
SimplexVolume<-function(x,Factorial=TRUE){
n<-length(x[,1])
d<-apply(x,2,FUN="diff")
if(Factorial){
ret<-log(abs(det(d))) - lfactorial(n-1)
}else{
ret<-log(abs(det(d)))
}
return(ret)
}
SimplexVolume2<-function(x,Factorial=TRUE){
n<-length(x[,1])
tmpX<-t(x[2:n,])-x[1,]
if(n==2){
ret<-log(sqrt(sum(tmpX^2)))
}else{
tmpX<-SubDfSVD(t(tmpX))
ret<-SimplexVolume(rbind(rep(0,length(tmpX[1,])),tmpX),Factorial=Factorial)
}
return(ret)
}
SubDfSVD<-function(RR){
svd.out<-svd(t(RR))
RR%*%svd.out[[2]]
}
library(sets)
library(rgl)
Nv<-10
df<-9
X<-matrix(rnorm(df*Nv),ncol=df)
nrm<-sqrt(apply(X^2,1,sum))
X<-X/nrm
X[,1]<-abs(X[,1])
X[,2]<-abs(X[,2])
EX<-matrix(runif(Nv*df),ncol=df)
EX<-CategoryVector(Nv)
EX<-X+matrix(rnorm(ncol(X)*nrow(X)),ncol=ncol(X))*1
df<-ncol(X)
library(GPArotation)
Nv<-10
df<-10
X<-matrix(0,Nv,df)
kkk<-4
for(i in 1:kkk){
X[,i]<-runif(Nv)
}
plot3d(X[,1:3])
RM<-Random.Start(df)
X<-X%*%RM
Nv<-20
df<-10
df1<-4
X<-matrix(rnorm(df1*Nv),ncol=df1)
nrm<-sqrt(apply(X^2,1,sum))
X<-X/nrm
plot3d(X)
X<-cbind(X,matrix(0,Nv,df-df1))
RM<-Random.Start(df)
X<-X%*%RM
Ngs<-list()
Factorial<-FALSE
GList<-list()
GVolume<-list()
GVolumeSum<-list()
k<-1
GList[[k]]<-list()
GVolume[[k]]<-list()
for(i in 1:Nv){
GList[[k]][[i]]<-as.set(c(i))
GVolume[[k]][[i]]<-0
}
ks<-1:df
for(kk in 1:length(ks)){
k<-ks[kk]
Ngs[[k]]<-length(GList[[k]])
currentGList<-GList[[k]]
currentNgs<-length(currentGList)
Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
for(i in 1:(currentNgs-1)){
for(j in (i+1):currentNgs){
dif1<-currentGList[[i]]-currentGList[[j]]
dif2<-currentGList[[j]]-currentGList[[i]]
if(length(dif1)==1 & length(dif2)==1){
Em[i,j]<-1
vs<-set_union(currentGList[[i]],currentGList[[j]])
Ws[i,j]<-Ws[j,i]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
}
}
}
cnt<-1
GList[[k+1]]<-list()
GVolume[[k+1]]<-list()
Candidates<-rep(1,Ngs[[k]])
NonCandidates<-rep(0,Ngs[[k]])
for(i in 1:(Ngs[[k]]-1)){
currentNgsX<-currentNgs
currentNgs<-length(currentGList)
NewCandidates<-rep(0,currentNgs)
NewNonCandidates<-rep(1,currentNgs)
NewCandidates[1:currentNgsX]<-Candidates
NewNonCandidates[1:currentNgsX]<-NonCandidates
Candidates<-NewCandidates
NonCandidates<-NewNonCandidates
if(i!=1){
WsX<-Ws
Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
Ws[1:nrow(WsX),1:ncol(WsX)]<-WsX
if(currentNgs>currentNgsX){
for(ii in which(Candidates==1)){
for(jj in max(currentNgsX+1,ii+1):currentNgs){
dif1<-currentGList[[ii]]-currentGList[[jj]]
dif2<-currentGList[[jj]]-currentGList[[ii]]
if(length(dif1)==1 & length(dif2)==1){
Em[ii,jj]<-1
vs<-set_union(currentGList[[ii]],currentGList[[jj]])
Ws[ii,jj]<-Ws[jj,ii]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
}
}
}
}
}else{
Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
for(ii in 1:(currentNgs-1)){
for(jj in (ii+1):currentNgs){
dif1<-currentGList[[ii]]-currentGList[[jj]]
dif2<-currentGList[[jj]]-currentGList[[ii]]
if(length(dif1)==1 & length(dif2)==1){
Em[ii,jj]<-1
vs<-set_union(currentGList[[ii]],currentGList[[jj]])
Ws[ii,jj]<-Ws[jj,ii]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
}
}
}
}
if(i==1){
tmpVolume<-min(Ws[upper.tri(Ws)])
GVolume[[k]][[i]]<-exp(tmpVolume)
selected<-which(Ws==tmpVolume,arr.ind=TRUE)
if(length(selected)>2)selected<-selected[1,]
Candidates[selected]<-0
NonCandidates[selected]<-1
}else{
tmpWs<-Ws[which(Candidates==1),which(NonCandidates==1)]
tmpVolume<-min(tmpWs)
tmpselected<-which(tmpWs==tmpVolume,arr.ind=TRUE)
if(length(tmpselected)>2)tmpselected<-tmpselected[1,]
GVolume[[k]][[i]]<-exp(tmpVolume)
if(i!=(Ngs[[k]]-1)){
selected<-c(which(Candidates==1)[tmpselected[1]],which(NonCandidates==1)[tmpselected[2]])
}else{
selected<-c(which(Candidates==1),which(NonCandidates==1)[tmpselected[1]])
}
Candidates[selected[1]]<-0
NonCandidates[selected[1]]<-1
}
GList[[k+1]][[cnt]]<-set_union(currentGList[[selected[1]]],currentGList[[selected[2]]])
tmpset<-set_intersection(currentGList[[selected[1]]],currentGList[[selected[2]]])
currentGListLen<-length(currentGList)
tmpcnt<-1
for(j in tmpset){
currentGList[[currentGListLen+tmpcnt]]<-GList[[k+1]][[cnt]]-j
tmpcnt<-tmpcnt+1
}
cnt<-cnt+1
}
GVolumeSum[[k]]<-sum(unlist(GVolume[[k]]))
}
GList
plot(X[,1],X[,2])
for(i in length(GList):2){
for(j in 1:length(GList[[i]])){
scmb<-set_combn(GList[[i]][[j]],2)
for(j2 in scmb){
v<-c()
for(j3 in j2){
v<-c(v,j3)
}
segments(X[v[1],1],X[v[1],2],X[v[2],1],X[v[2],2],col=i)
}
}
}
plot3d(X[,1],X[,2],X[,3])
for(i in length(GList):2){
for(j in 1:length(GList[[i]])){
scmb<-set_combn(GList[[i]][[j]],2)
for(j2 in scmb){
v<-c()
for(j3 in j2){
v<-c(v,j3)
}
segments3d(X[v,],col=i)
}
}
}
GVolumeSum
incr<-unlist(GVolumeSum)[2:length(GVolumeSum)]/unlist(GVolumeSum)[1:(length(GVolumeSum)-1)]
par(mfcol=c(2,2))
plot(incr)
plot(unlist(GVolumeSum))
plot(GVolumeSum[[1]]^(1:(length(GVolumeSum)-1)))
plot(unlist(GVolume))
par(mfcol=c(1,1))