ただのメモ
角座標・多次元球・カイ自乗分布・非心カイ自乗分布・αエラー・パワーなど
library(sphere) ThetaAngle<-function(cos,sin){ t<-acos(cos) t[which(sin<0)]<-2*pi-t[which(sin<0)] t } Simplex<-function(Ts,a,r){ A<-Ts%*%t(a) R<-Ts%*%r L<-sqrt(A^2+R^2) Cs<-A/L Ss<-R/L Thetas<-ThetaAngle(Cs,Ss) #Thetas[which(Ss<0)]<-2*pi -Thetas[which(Ss<0)] TPointsX<-1/L * Cs TPointsY<-1/L * Ss xlim<-ylim<-c(-10,10) plot(TPointsX,TPointsY,xlim=xlim,ylim=ylim) selected<-c(which(A==max(A))) current<-selected[1] VsX<-c(1/max(A)) VsY<-c(0) loop<-TRUE while(loop){ NegativeStart<-FALSE if(length(selected)==1 & Thetas[current]>pi){ NegativeStart<-TRUE Mores<-c(which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)), which(Thetas>=0 & Thetas <(Thetas[current]-pi)) ) }else{ NegativeStart<-FALSE Mores<-which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)) } #if(Thetas[current]<pi){ #Mores<-which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)) #}else{ #Mores<-c(which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)), which(Thetas>=0 & Thetas <(Thetas[current]-pi)) ) #} tmpL<-sqrt((1/L[Mores]*Cs[current]-1/L[current]*Cs[Mores])^2 + (1/L[Mores]*Ss[current]-1/L[current]*Ss[Mores])^2 ) tmpSs<-(-1/L[current]*Cs[Mores]+1/L[Mores]*Cs[current]) / tmpL tmpCs<-(1/L[current]*Ss[Mores]-1/L[Mores]*Ss[current]) / tmpL tmpThetas<-ThetaAngle(tmpCs,tmpSs) if(Thetas[current]+pi/2-2*pi>0){ tmpThetas[which(tmpThetas<Thetas[current]+pi/2-2*pi)]<-tmpThetas[which(tmpThetas<Thetas[current]+pi/2-2*pi)]+2*pi } #print(tmpThetas) #print(NegativeStart) #if(!NegativeStart){ #tmpThetas[which(tmpThetas<Thetas[current])]<-2*pi-tmpThetas[which(tmpThetas<Thetas[current])] #tmpThetas[which(tmpThetas>pi]<-2*pi-tmpThetas[which(tmpThetas>pi] #} #if(Thetas[current]>pi){ #tmpThetas<-2*pi -tmpThetas #} #tmpThetas[which(tmpSs<0)]<-2*pi - tmpThetas[which(tmpSs<0)] #tmpMoreThetas<-tmpThetas[which(tmpThetas>Thetas[current])] if(length(tmpThetas)>0){ tmpid<-which(tmpThetas==min(tmpThetas)) tmpselected<-Mores[tmpid] selected<-c(selected,tmpselected) rx<-1/L[current]/cos(tmpThetas[tmpid]-Thetas[current]) VsX<-c(VsX,rx*cos(tmpThetas[tmpid])) VsY<-c(VsY,rx*sin(tmpThetas[tmpid])) current<-tmpselected }else{ loop<-FALSE } if(current==selected[1]){ loop<-FALSE } #print(selected) #plot(1/L[selected]*Cs[selected],1/L[selected]*Ss[selected]) for(i in 1:length(selected)){ abline(1/(L[selected[i]]*Ss[selected[i]]),-Cs[selected[i]]/Ss[selected[i]],col="red") } } #print(Thetas[selected]) crossX<-crossY<-rep(0,length(selected)+1) for(i in 1:length(selected)){ x1<-1/L[selected[i]]*Cs[selected[i]] y1<-1/L[selected[i]]*Ss[selected[i]] x2<-1/L[selected[i+1]]*Cs[selected[i+1]] y2<-1/L[selected[i+1]]*Ss[selected[i+1]] crossX[i]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2) crossY[i]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2) } x1<-1/L[selected[length(selected)]]*Cs[selected[length(selected)]] y1<-1/L[selected[length(selected)]]*Ss[selected[length(selected)]] x2<-1/L[selected[1]]*Cs[selected[1]] y2<-1/L[selected[1]]*Ss[selected[1]] crossX[length(crossX)]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2) crossY[length(crossX)]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2) list(L=L,Cs=Cs,Ss=Ss,selected=selected,VsX=crossX,VsY=crossY) } #Simplex2 #最短辺を選ぶ #ついで、現在の法線ベクトルから反時計周りにpiの法線ベクトルを持つ辺を対象とし #交点と、現在の法線ベクトルとのなす余弦が最大になるような辺を次の辺として採用する #最短辺に戻ったら終わり Simplex2<-function(Ts,a,r){ A<-Ts%*%t(a) R<-Ts%*%r L<-sqrt(A^2+R^2) Cs<-A/L Ss<-R/L Thetas<-ThetaAngle(Cs,Ss) #Thetas[which(Ss<0)]<-2*pi -Thetas[which(Ss<0)] TPointsX<-1/L * Cs TPointsY<-1/L * Ss xlim<-ylim<-c(-10,10) plot(TPointsX,TPointsY,xlim=xlim,ylim=ylim) #selected<-c(which(A==max(A))) selected<-which(L==max(L))[1] StartTheta<-Thetas[selected[1]] Theta2<-(Thetas-StartTheta)%%(2*pi) #Theta2[selected[1]]<-2*pi current<-selected[1] #VsX<-c(1/max(A)) #VsY<-c(0) VsX<-c() VsY<-c() loop<-TRUE #first<-TRUE while(loop){ #if(!first){ #Theta2[selected[1]]<-2*pi #} #if(first){ #first<-FALSE #} currentTheta2<-Theta2[current] #print(current) plusPiTheta2<-currentTheta2+pi EndTheta2<-min(2*pi,plusPiTheta2) #print(EndTheta2) if(EndTheta2==(2*pi)){ # 最初の法線からpi以上回ったら、最初の法線も考慮の対象に入れるため Theta2[selected[1]]<-(currentTheta2+EndTheta2)/2 } Mores<-which(Theta2>currentTheta2 & Theta2<EndTheta2) tmpL<-sqrt((1/L[Mores]*Cs[current]-1/L[current]*Cs[Mores])^2 + (1/L[Mores]*Ss[current]-1/L[current]*Ss[Mores])^2 ) tmpSs<-(-1/L[current]*Cs[Mores]+1/L[Mores]*Cs[current]) / tmpL tmpCs<-(1/L[current]*Ss[Mores]-1/L[Mores]*Ss[current]) / tmpL if(length(tmpCs)>0){ ttcc<-tmpCs * Cs[current] + tmpSs * Ss[current] tmpid<-which(ttcc==max(ttcc))[1] #tmpid<-which(tmpThetas==min(tmpThetas)) tmpselected<-Mores[tmpid] selected<-c(selected,tmpselected) current<-tmpselected }else{ loop<-FALSE } if(current==selected[1]){ loop<-FALSE } #print(selected) #plot(1/L[selected]*Cs[selected],1/L[selected]*Ss[selected]) #for(i in 1:length(selected)){ #abline(1/(L[selected[i]]*Ss[selected[i]]),-Cs[selected[i]]/Ss[selected[i]],col="red") #} } #print(Thetas[selected]) crossX<-crossY<-rep(0,length(selected)-1) for(i in 1:(length(selected)-1)){ x1<-1/L[selected[i]]*Cs[selected[i]] y1<-1/L[selected[i]]*Ss[selected[i]] #if(i<length(selected)){ x2<-1/L[selected[i+1]]*Cs[selected[i+1]] y2<-1/L[selected[i+1]]*Ss[selected[i+1]] #}else{ #x2<-1/L[selected[1]]*Cs[selected[1]] #y2<-1/L[selected[1]]*Ss[selected[1]] #} crossX[i]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2) crossY[i]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2) } list(L=L,Cs=Cs,Ss=Ss,selected=selected,VsX=crossX,VsY=crossY) } # 2次元凸包をある点から眺めた角度で考える # 少しずつ角度を変えてそのときの、原点からの距離と # 視点からの距離とを出し # そのリストを返す DistSimplex<-function(vs,p,n){# vs:(x,y)coordinates, p:new center, n: 各辺をn等分(角として)する newvs<-t(t(vs)-p) r<-sqrt(apply(newvs^2,1,sum)) #outx<-outy<-rep(0,n*length(vs[,1])) outx<-outy<-matrix(0,length(vs[,1]),n) counter<-0 deltatheta<-rep(0,length(vs[,1])) for(i in 1:length(vs[,1])){ #for(i in 1:18){ st<-i end<-i+1 if(i==length(vs[,1])){ end<-1 } if(r[st]!=0 & r[end]!=0){ cosst<-newvs[st,1]/r[st] sinst<-newvs[st,2]/r[st] cosend<-newvs[end,1]/r[end] sinend<-newvs[end,2]/r[end] tst<-ThetaAngle(cosst,sinst) tend<-ThetaAngle(cosend,sinend) difft<-tend-tst if(difft>=0 & difft <=pi){ }else if(difft > pi & difft<=2*pi){ difft<-(-1)*(2*pi - difft) }else if(difft<0 & difft >=-pi){ }else if(difft< -pi & difft >= -2*pi){ difft<-2*pi+difft } deltatheta[i]<-difft #if(difft<0){ #difft<-difft+2*pi #} #print(difft) ts<-tst+difft*(0:(n-1))/n + difft/(2*n) newt1<-tst-ts newt2<-tend-ts newrs<-r[st]*r[end]/(r[st]*sin(newt1)-r[end]*sin(newt2)) * sin(difft) #if(i!=0){ #par(new=TRUE) #} #plot(newrs*cos(ts),newrs*sin(ts)) #outx[(counter+1):(n*i)]<-newrs*cos(ts) #outy[(counter+1):(n*i)]<-newrs*sin(ts) outx[i,]<-newrs*cos(ts) outy[i,]<-newrs*sin(ts) } counter<-n*i #plot(outx,outy) } return(list(x=outx+p[1],y=outy+p[2],p=p,n=n,rp=sqrt(outx^2+outy^2),deltatheta=deltatheta)) } df<-5 Nt<-8 # 法線ベクトルの集合 Ts<-RandomSphere(n=Nt,df=df) Ts<-rbind(Ts,-Ts) # 真実のモデルの方向の単位ベクトル a<-RandomSphere(n=1,df=df) # 輪切りにする面の数 Ns<-3 # その単位ベクトル rs<-RandomSphere(n=Ns,df=df) ip<-rs %*% t(a) rs<-t(t(rs)-t(a)%*%t(ip)) rs<-rs / sqrt(apply(rs^2,1,sum)) plot(as.data.frame(rs)) # 統計量の大きさとモデルの強さとの比として計算する数(Nr*2+1) Nr<-8 Ratios<-10^((-Nr):Nr) # 統計量とモデルの強さの比あたり、計算する、角度の数 Na<-12 Angles<-seq(from=0,to=1,by=1/(Na-1))*2*pi for(i in 1:length(rs[,1])){ simplexOut<-Simplex2(Ts,a,rs[i,]) #simplexOut<-Simplex(Ts,a,r) #simplexOut xlim<-ylim<-c(min(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss), max(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss)) xlim<-ylim-c(min(simplexOut$VsX,simplexOut$VsY),max(simplexOut$VsX,simplexOut$VsY)) xlim<-ylim<-c(-4,4) plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20) for(i in 1:length(simplexOut$selected)){ abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red") } par(new=TRUE) plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue") } plot(simplexOut$VsX,simplexOut$VsY,pch=19,col="blue") for(i in 1:Ns){ #単位凸包を作る simplexOut<-Simplex2(Ts,a,rs[i,]) #print(UnitSimplex) xlim<-ylim<-c(-10,10) plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20) for(i in 1:length(simplexOut$selected)){ abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red") } par(new=TRUE) plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue") #相対的な位置を変えたとき(Ratiosの値)の関数として表したい #アウトプットは、2交点の座標なのだが・・・虚数も含めて-b \pm sqrt(D) /2a のDとかを使って表す } ################ # A/K ratiosを出す # 等比級数にすることで、計算量を減らす geomSeries<-function(stval,endval,n){ exp(seq(from=log(stval),to=log(endval),length.out=n)) } maxp<-10^(-1) minp<-10^(-15) minchival<-sqrt(qchisq(maxp,df=df,lower.tail=FALSE)) maxchival<-sqrt(qchisq(minp,df=df,lower.tail=FALSE)) n<-10 geomSer<-geomSeries(minchival,maxchival,n) logchisqval<-seq(from=log(minchisqval),to=log(maxchisqval),length.out=10) alphachival<-exp(logchisqval) betachival<-c(0,alphachival) # betachival/alphachival の値は、alphachivalとbetachivalが同じ等比数列なので、比の-(要素数)乗・・・(要素数)乗になる ratio<-alphachival[1]/alphachival[2] ratios<-ratio^((-length(alphachival)+1):(length(alphachival)-1)) outprob<-matrix(0,length(ratios),length(betachival)) for(x in 1:Ns){ simplexOut<-Simplex2(Ts,a,rs[x,]) # 凸包の形は決まった vs<-cbind(simplexOut$VsX,simplexOut$VsY) #area<-rep(0,length(ratios)) for(i in 1:length(ratios)){ p<-c(ratios[i],0) n<-100 # 視点からの等角点の計算 dsout<-DistSimplex(vs,p,n) xlim<-ylim<-c(min(dsout$x,dsout$y,p),max(dsout$x,dsout$y,p)) plot(c(t(dsout$x)),c(t(dsout$y)),type="l",xlim=xlim,ylim=ylim) par(new=TRUE) plot(p[1],p[2],col="red",xlim=xlim,ylim=ylim) #area[i]<-sum((pchisq(dsout$rp^2,df=df))*dsout$deltatheta/(2*pi*dsout$n)) for(j in 1:length(outprob[1,])){ outprob[i,j]<-outprob[i,j]+sum((pchisq((dsout$rp*betachival[j])^2,df=df))*dsout$deltatheta/(2*pi*dsout$n)) } } matplot(outprob,type="l") } outprob<-outprob/Ns ###############