等分するベクトル2

  • 4日目。一昨昨日、一昨日と昨日(記事はこちらこちらこちら)の続き
  • \sum_{i=1}^{df+1} x_i^2=1を満足する均一分布を発生させる
  • そのうち、\sum_{i=1}^{df+1} \sqrt{N_i} x_i=0(周辺度数制約空間で原点を通っている)を満足させる点のみを考えたい
  • df+1次元球面上の点から、周辺度数制約空間に垂線を下ろす
  • 球面の点をF=\{f_1,...,f_{df+1}\}、垂線の足をF'=\{f_i',...,f_{df+1}'\}とすると
    • F-F' = k (\sqrt{N_1}.\sqrt{N_2},...\sqrt{N_{df+1}})となって、かつ、
    • \sum_{i=1}^{df+1} \sqrt{N_i} f_i'=0を満足する
  • これをkについて解くことができて、F'の座標は定まる
  • 今、F'はdf+1次元単位球と周辺度数制約空間との交わりの空間で、あるが、単位球の内部の点となっているので、これを、原点から外側に伸ばしてやって、単位球表面の点に置きなおす
  • このようにすることで、単位球面上であって、制約空間に含まれる部分空間に均一分布の点を発生させることができている(はず)。それを描くRのソース
df<-4
np<-10000
rs<-matrix(rnorm((df+1)*np),nrow=np)
rssq<-rs^2
ds<-sqrt(apply(rssq,1,sum))
rs<-rs/ds

weight<-round(runif(df+1)*100,0)
#weight<-rep(1,df+1)
#weight<-c(1,0,0)

tmp<-t(t(sqrt(weight))%*%t(rs))/sum(weight)
rs2<-rs-tmp%*%t(sqrt(weight))
#plot(as.data.frame(rs2))


#rs2<-rs-(t(weight)%*%t(rs))/sum(weight)

rssq<-rs2^2
ds<-sqrt(apply(rssq,1,sum))
rs2<-rs2/ds

plot(as.data.frame(rs2[,1:min(df,5)]),cex=0.1)