ただのメモ

角座標・多次元球・カイ自乗分布・非心カイ自乗分布・αエラー・パワーなど

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

###############