2011-04-04

球は入れ子

R 球多次元球

昨日、「同じ $r^2$ のマーカーの取り方と多次元球・正単体について書いた
検定のパワーとＬＤのことを考えるための導入だ
話を進める前に多次元球とその球面上の減次元球について考えよう
k次元球があり、その外部の点から、球を見ると、太陽を地球から眺めているように見える。遠くから見れば小さい丸に見えるし、近づくと、圧倒される
この太陽の大きさを決めているのは、「視角」。そんな話
さて。
- 次元はk
- 的は、原点Oを中心とした半径１の球面である
- 外側の点Xは、原点からの距離がx, $x>1$ である
- 今、的上の点Pがあって、角POXがtであるという
- Xの座標を $(1,0,0,...,0)$ とし、 $(\cos(t),\sin(t),0,0,...,0)$ とする
- 角PXO s は、的の中心を狙ったときの、ずれの角度である
- $\tan(s)=\frac{\sin{t}}{x-\cos(t)}$ となる
- 的上の点であって、POXがtとなるような点は、 $(\cos(t),....)$ なOXに垂直な面上にある半径 $\sin(t)$ なるk-1次元球面を作る
- このk-1次元球面の面積は(OXから角tの点の「面積」)は半径なるk-1次元球面であるから、それは
  - $2\frac{\pi^{\frac{k-1}{2}}}{\Gamma(\frac{\pi}{2})}$ となる
- k次元球面は、このようにk-1次元球面が連続して並んでいる
- k-1次元球面はk-2次元球面が並んだもの…というように入れ子になっている
- k次元球のある軸から、偏角tが作るk-1次元球の球面積を求める関数を作ろう
- 2次元面にある単位球(単位円)の円周を考える
- 円の中心を原点とする角座標で $0\le t \le 2\pi$ を考える
- 微小角に対応する円周はだから、円周をに関するからの積分と考えれば
  - $\int_0^{2\pi} 1 dt=2\pi$
- 3次元空間にある単位球の表面積を考える
- 球の中心を原点とし、ある方向を一つ定めその方向のベクトルを法線ベクトルとする面で球を輪切りにすることを考える
- 輪切り面は円で交わるが、この交円上の点の位置ベクトルが法線ベクトルとなす角であるような円がの範囲で存在する
  - $\int_0^{\pi} 2\pi\sin(t)dt = 2\pi (-(\cos(\pi)-\cos(0)))=4\pi$
- 4次元空間にある単位球の表面積を考える
- 球の中心を原点とし、ある方向を一つ定めその方向のベクトルを法線ベクトルとする面で４次元球を輪切りにすることを考える
- 輪切り面は３次元球で交わるが、この交球上の点の位置ベクトルが法線ベクトルとなす角であるような3次元球がの範囲で存在する
  - $\int_0^{\pi} 4\pi \sin^2(t) dt = 4\pi \int_0^{\pi} \frac{1}{2}(1-\cos(2t)) dt = 4\pi \frac{1}{2}(\pi-0)=2 \pi^2$

# df次元球の表面積
SphereSurface<-function(df,r=1){
	2*r^df*pi^(df/2)/gamma(df/2)
}
# df次元球の軸と偏角tであるdf-1次元球の表面積
SphereSurfaceAngle<-function(df,t,r=1){
	d<-r*abs(sin(t))
	SphereSurface(df-1,d)
}
SphereSurfaceAngle(3,pi/2)

今、角POX,tに関して考えた
角PXO, sについて考えるとどうなるだろうか
- 的のことを考えずにXから(何かを)まっすぐに投げたとき、そもそも、当たるか当たらないか、が問題になる
- そのうえで、当たるとしたらどこ(tがいくつのところ)にあたるか、はどうなるのだろうか…
２次元空間に単位円周の的があって、円周の外の点からそこに向かって的当てをするときの、あたりやすさを考える
- 円周状の的は、実は直線の張りぼてで、凸多角形だとする
- このようなときの当たりやすさも同様に考えられる
- さて、この２次元空間の凸多角形をより高次な次元にろくろを回すようにぐるぐる回して「ある意味では凸多角形で、ある意味では、丸い」ような的にする
- この当たりやすさは、２次元の凸多角形の的の当たりやすさを評価しておけば、当的位置から的の位置への偏角と次元の係数とを使って、任意次元の「凸多角形的性格を持つ丸い的」への当たりやすさを計算できる
- そのことを計算するソースが以下。ただし、「÷0」とかがあって、振動する部分が制御できていない
- また、これをやると、「多次元空間」を2次元に関する「網羅」のみで済む(かもしれない)点は、素敵なこと

# 指定したKvとTiisに関して
# 指定したKvとTiisに関して
# 指定した(複数の)Ksx(複数の)CheckDistsに関する、交点のリストを
# nperm試行別に出力
# ただし、交点の数は0,1,2のいずれかになっている
# この出力から、任意のdfでのパワーが計算できる

SpherePower2DBasic<-function(Kv,Ks,CheckDists,Tiis,nperm=1000){
	Ks[which(Ks==0)]<-0.0001
	CheckDists[which(CheckDists==0)]<-0.0001
	Kss<-matrix(rep(Kv,length(Ks)),nrow=length(Ks),byrow=TRUE)*Ks
	Ntest<-length(Tiis[,1])
	#rss<-RandomSphere(df=df,n=nperm)
	rt<-runif(nperm)*2*pi
	rt<-sort(rt)
	#wt<-SphereSurfaceAngle(df,rt)
	#wt<-wt/sum(wt)
	rss<-cbind(cos(rt),sin(rt))
	#PowerOut<-matrix(0,length(Ks),length(CheckDists))
	PowerOut<-array(0,c(length(Ks),length(CheckDists),nperm,2))
	for(ik in 1:length(Ks)){
		Pii<-Kss[ik,]
		AccumProbs <- matrix(0, nperm, length(CheckDists))
		for (i in 1:nperm) {
			rs <- matrix(rss[i,],nrow=1)
			tmptheta <- acos(cosxy(t(Pii), rs))
			C <- nrm(Pii) * cos(tmptheta)
			D <- nrm(Pii) * sin(tmptheta)
			As <- t(Tiis %*% t(rs))
			Bs <- t(Tiis %*% Pii)
			current <- which(Bs == max(Bs))[1]
			SerialTopT <- c(current)
			currentR <- 0
			SerialR <- c(currentR)
			SerialStat <- c(Bs[current])
			CheckList <- 1:Ntest
			CheckList <- CheckList[which(CheckList != current)]
			while (length(CheckList) > 0) {
				tmpks <- -(Bs[current] - Bs[CheckList])/(As[current] - 
				As[CheckList])
				tmpList <- which(tmpks > currentR)
				if (length(tmpList) > 0) {
					current <- CheckList[which(tmpks == min(tmpks[tmpList]))][1] 
					currentR <- min(tmpks[tmpList]) 
					SerialTopT <- c(SerialTopT, current) 
					SerialR <- c(SerialR, currentR) 
					SerialStat <- c(SerialStat, As[current] * currentR + Bs[current]) 
				}	
				CheckList <- CheckList[tmpList] 
				CheckList <- CheckList[which(CheckList != current)] 
			}
			intMyX<-intMyX3 <- matrix(0,length(CheckDists),length(SerialTopT))
			AnswerPs<-rep(0,length(CheckDists))
			for(j in 1:length(CheckDists)){
				tmpSerialStat<-SerialStat-CheckDists[j]
				if(tmpSerialStat[1]>=0){
					AnswerPs[j]<-1
				}
				if(length(tmpSerialStat)>1){
					for(jj in 1:(length(tmpSerialStat)-1)){
						if(tmpSerialStat[jj]==0){
							intMyX[j,jj]<-1
							intMyX3[j,jj]<-sign(As[SerialTopT[jj]])
						}else{
							if(tmpSerialStat[jj]*tmpSerialStat[jj+1]<0){
								intMyX[j,jj]<-1
								intMyX3[j,jj]<-sign(As[SerialTopT[jj]])
							}
						}
					}
				}
				
				if(tmpSerialStat[length(tmpSerialStat)]==0){
					intMyX[j,length(tmpSerialStat)]<-1
					intMyX3[j,length(tmpSerialStat)]<-sign(As[SerialTopT[length(tmpSerialStat)]])
				}else{
					if(tmpSerialStat[length(tmpSerialStat)]*As[SerialTopT[length(tmpSerialStat)]]<0){
						intMyX[j,length(tmpSerialStat)]<-1
						intMyX3[j,length(tmpSerialStat)]<-sign(As[SerialTopT[length(tmpSerialStat)]])
					}
				}
			}


			zeros <- which(apply(abs(intMyX), 1, sum) == 0) 
			Crossings <- t(outer(CheckDists, Bs[SerialTopT], FUN = "-"))/As[SerialTopT]
			# ここが変！
			#print(Crossings)
			#print(intMyX3)
			for(cc in 1:length(Crossings[1,])){
				tmpid<-which(intMyX3[cc,]!=0)
				#tmpid<-order(abs(Crossings[,cc]))
				if(length(tmpid)>2){
						print(tmpid)
				}
			
				#print(tmpid[1:2])
				if(length(tmpid)>0){
					PowerOut[ik,cc,i,1:length(tmpid)]<-Crossings[tmpid,cc]
					#PowerOut[ik,cc,i,1:min(length(tmpid),2)]<-Crossings[tmpid[1:min(length(tmpid),2)],cc]
					#print("---")
					#print(PowerOut[ik,cc,i,1:min(length(tmpid),2)])
				}
			}
			
			#ChiProbs <- pchisq(Crossings^2, df, lower.tail = FALSE)
			#AnswerPs <- AnswerPs + apply(t(intMyX3) * ChiProbs, 2, sum)
			#AnswerPs[zeros]<-AnswerPs[zeros] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 1)/2
			#AnswerPs[which(AnswerPs<0)]<-1+AnswerPs[which(AnswerPs<0)]
			#AnswerPs[which(AnswerPs>1)]<-AnswerPs[which(AnswerPs>1)]-1
			#AccumProbs[i, ] <- AnswerPs
		}
		#PowOut <- apply(AccumProbs, 2, sum)/nperm
		#PowOut <- apply(AccumProbs*wt, 2, sum)
		#PowerOut[ik,]<-PowOut
		#PowerOut[ik,,]<-PowOut
	}

	list(Kv=Kv,Ks=Ks,CheckDists=CheckDists,Tiis=Tiis,nperm=nperm,crosses=PowerOut,rt=rt)
}

Kv<-c(1,0)
Ks<-seq(from=3,to=4,length.out=6)
Ks<-c(3.1)
CheckDists<-seq(from=1,to=7,length.out=6)
#CheckDists<-c(3)
library(MCMCpack)

Tiis<-rdirichlet(10,rep(1,2))
Tiis<-rbind(Tiis,-Tiis)
nperm<-1000
ssppbb<-SpherePower2DBasic(Kv=Kv,Ks=Ks,CheckDists=CheckDists,Tiis=Tiis,nperm=nperm)

 matplot(pchisq(ssppbb$crosses[1,2,,]^2,2,lower.tail=FALSE),type="l")
 matplot(ssppbb$crosses[1,2,,],type="l")