球面配置

  • 昨日は球の分裂について書いた
  • 分裂の均等性・不均等性を議論するには、分裂後の球に引き継がれた点の分布を球面上で評価しなくてはならない
  • そのために球表面の分布が取りたい
  • 球表面をどんな風に分割するのがよいだろうか、と考える
  • こんなことを考えたこともあった
  • まずは簡単のために円上の度数分布を考えたい
  • 普通の1次元尺度の度数分布では、等幅にbinを取る
  • 多尺度の度数分布では、通常、格子状にbinを取る
  • いずれも、一つのbinは同じ幅・面積を持つ
  • 円について角座標を使ってbinを作るときも、個々のbinの面積は同じにしたい
  • 角座標にしているということは、「原点からの距離」と「方向」とでbinを作りたいということ
  • 中央に円をとり、そこから同心円状に帯を作り、帯を角度で均等割りして、結果として、すべてのbinが等しい面積を持たせたい
  • いくつの同心円状の帯に分けるか、帯が一つ外側になるにつれ、どれくらい角による分割数を増やすかを決めると、同心円の半径と分割の角が定まる
  • 後は、binに入る点を数え上げるだけ
  • この後、3次元球表面に拡張しなくては…
X <- matrix(runif(100000*2)*2-1,ncol=2)
X <- matrix(c(rnorm(10000,0.5,0.3),rnorm(10000,0.2,0.8)),ncol=2)
r <- sqrt(apply(X^2,1,sum))
s <- which(r<=1)
X <- X[s,]

plot(X)

t <- 4
k <- 6
my.S <- function(t,k){
	(t^(k+1)-1)/(t-1)
}
r <- sqrt(my.S(t,0:k)/my.S(t,k))

plot(r)

thetas <- list()
for(i in 1:(k+1)){
	thetas[[i]] <- seq(from = 0, to = 2*pi, length =t^(i-1))
}

Sphere.Histogram.2 <- function(X,t,k){
	r <- c(0,sqrt(my.S(t,0:k)/my.S(t,k)))
	#print(r)
	Hist <- list()
	
	thetas <- list()
	for(i in 1:k){
		thetas[[i]] <- seq(from = 0, to = 2*pi, length =t^(i-1)+1)
		Hist[[i]] <- rep(0,length(thetas[[i]]-1))
	}
	#print(thetas)
	Cmpx <- complex(real = X[,1],imaginary = X[,2])
	Md <- Mod(Cmpx)
	Ag <- Arg(Cmpx)
	neg <- which(Ag < 0)
	Ag[neg] <- 2*pi +Ag[neg]
	for(i in 1:length(Hist)){
		s <- which(Md >= r[i] & Md < r[i+1])
		a <- Ag[s] %/% (thetas[[i]][2]) +1
		#print(i)
		#print(range(Ag[s]))
		#print(a)
		Hist[[i]] <- tabulate(a,nbins = t^(i-1)+1)
		Hist[[i]] <- Hist[[i]][-length(Hist[[i]])]
	}
	return(list(Hist=Hist,X=X,Md=Md,Ag=Ag,t=t,k=k,r=r,thetas=thetas))
}
t <- 4
k <- 4
sp.h <- Sphere.Histogram.2(X,t,k)

x <- y <- c()
freq <- c()
for(i in 1:(length(sp.h$r)-2)){
	for(j in 1:(length(sp.h$thetas[[i]])-1)){
		x <- c(x, sp.h$r[i]*cos(sp.h$thetas[[i]][j]))
		y <- c(y, sp.h$r[i]*sin(sp.h$thetas[[i]][j]))
		freq <- c(freq,sp.h$Hist[[i]][j])
	}
}

plot(x,y,col=gray((max(freq)-freq)/max(freq)))
  • 3次元球面に拡張しよう
    • (0,0,1)を北極として、緯度を下げて南極まで進む。緯度に関する輪切りによって、球表面を分割する。各輪切りは、1,t,t^2,..,t^k,t^k,t^{k-1},t^{k-2},...,t,1に角によって分割する。このように分割したら、すべての区画が同面積にするためには…
    • まず、北極から南極まで\piの緯度角度で輪切りにするときのその角\thetaが作る球表面の円を縁として北極を中心とする部分球面の面積が1-\cos \thetaであることを利用する
    • その点だけ変えて、あとは2次元のときとほぼ同じ(ただし、赤道を越えた南極部分は北極部分と対称で、分割数は南極に向かって減り始めるし、輪切りを定める緯度の増分も北半球・南半球で対称である
n.pt <- 1000

X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
apply(X^2,1,sum)

#X <- X[which(X[,3]>=0),]

R <- sqrt(X[,1]^2+X[,2]^2)

the <- asin(R)

t <- seq(from=0,to=pi/2,length=100)

L <- rep(0,length(t))
for(i in 1:length(t)){
	L[i] <- length(which(R <= sin(t[i])))
}

plot(t,L)

X.renew <- X
rs <- thetas <- phis <- rep(0,length(X.renew[,1]))
for(i in 1:length(rs)){
	tmp <- my.sphere.coord(X.renew[i,])
	rs[i] <- tmp$r
	thetas[i] <- tmp$theta
	phis[i] <- tmp$phi
}
sss <- which(thetas <= pi/2)
plot(thetas[sss]*cos(phis[sss]),thetas[sss]*sin(phis[sss]))


Sphere.Histogram.3 <- function(X,t,k){
	phis <- acos((t^(k)-t^(0:(k)))/(t^(k)-1))
	diff.phis <- diff(phis)
	phis <- c(phis[-length(phis)],cumsum(c(phis[length(phis)],diff.phis[length(diff.phis):1])))
	#print(r)
	Hist <- list()
	
	thetas <- list()
	for(i in 1:k){
		thetas[[i]] <- seq(from = 0, to = 2*pi, length =t^(i-1)+1)
		Hist[[i]] <- rep(0,length(thetas[[i]]-1))
	}
	tmpks <- k:1
	for(i in 1:length(tmpks)){
		thetas[[k+i]] <- thetas[[k+1-i]]
		Hist[[k+i]] <- Hist[[k+1-i]]
	}
	#print(thetas)
	Cmpx <- complex(real = X[,1],imaginary = X[,2])
	Md <- acos(X[,3])
	#Md <- Mod(Cmpx)
	Ag <- Arg(Cmpx)
	neg <- which(Ag < 0)
	Ag[neg] <- 2*pi +Ag[neg]
	for(i in 1:length(Hist)){
		s <- which(Md >= phis[i] & Md < phis[i+1])
		a <- Ag[s] %/% (thetas[[i]][2]) +1
		#print(i)
		#print(range(Ag[s]))
		#print(a)
		if(i<=k){
			Hist[[i]] <- tabulate(a,nbins = t^(i-1)+1)
		}else{
			Hist[[i]] <- tabulate(a,nbins = t^(k-1-i+k+1)+1)
		}
		Hist[[i]] <- Hist[[i]][-length(Hist[[i]])]
	}
	return(list(Hist=Hist,X=X,Md=Md,Ag=Ag,t=t,k=k,phis=phis,thetas=thetas))
}
n.pt <- 100000

X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
apply(X^2,1,sum)

#X <- X[which(X[,3]>=0),]

R <- sqrt(X[,1]^2+X[,2]^2)

the <- asin(R)

t <- seq(from=0,to=pi/2,length=100)

L <- rep(0,length(t))
for(i in 1:length(t)){
	L[i] <- length(which(R <= sin(t[i])))
}

plot(t,L)


sp.out3 <- Sphere.Histogram.3(X,2,3)
sp.out3$Hist