SphereSurface<-function(df,r=1){
2*r^df*pi^(df/2)/gamma(df/2)
}
SphereSurfaceAngle<-function(df,t,r=1){
d<-r*abs(sin(t))
ret<-0
if(df==1){
ret<-rep(2,length(t))
}else{
ret<-SphereSurface(df-1,d)
}
ret
}
SphereSurfaceAngle(3,pi/2)
SpherePowerND<-
function (df,Kv, Ks, CheckDists, Tiis, nperm = 1000)
{
dfin<-df
df<-2
Ks[which(Ks == 0)] <- 1e-04
Kss <- matrix(rep(Kv, length(Ks)), nrow = length(Ks), byrow = TRUE) *
Ks
Ntest <- length(Tiis[, 1])
tmprss<-seq(from=0,to=1,length.out=nperm)*2*pi
rss <- cbind(cos(tmprss),sin(tmprss))
ws<-SphereSurfaceAngle(dfin-1,tmprss)
ws<-ws/sum(ws)
PowerOut <- matrix(0, length(Ks), length(CheckDists))
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]
ChiProbs <- pchisq(Crossings^2, dfin, 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*ws[i]
}
PowOut <- apply(AccumProbs, 2, sum)
PowerOut[ik, ] <- PowOut
}
list(df = dfin, Kv = Kv, Ks = Ks, CheckDists = CheckDists,
Tiis = Tiis, nperm = nperm, p.out = PowerOut)
}
rs<-seq(from=0,to=1,length.out=101)
us<-rs*pi/2
Ks<-5
CheckDists<-5
nperm<-1000
dfs<-c(2,3,4,5,10)
power2markerdfs<-matrix(0,length(dfs),length(us))
for(x in 1:length(dfs)){
df<-dfs[x]
ppout<-rep(0,length(us))
Tiis<-matrix(0,2,2)
Tiis[1,1]<-1
Kv<-c(1,0)
for(i in 1:length(us)){
Tiis[2,1]<-cos(us[i])
Tiis[2,2]<-sin(us[i])
tmp<-rbind(Tiis,-Tiis)
spout<-SpherePowerND(df,Kv,Ks,CheckDists,Tiis,nperm)
ppout[i]<-c(spout$p.out)
}
power2markerdfs[x,]<-ppout
}
matplot(t(power2markerdfs2),type="l")
power2markerdfs2<-matrix(0,length(dfs),length(us))
for(x in 1:length(dfs)){
df<-dfs[x]
ppout<-rep(0,length(us))
Tiis<-matrix(0,3,2)
Tiis[1,1]<-1
Kv<-c(1,0)
for(i in 1:length(us)){
Tiis[2,1]<-cos(us[i])
Tiis[2,2]<-sin(us[i])
Tiis[3,1]<-cos(us[i])
Tiis[3,2]<--sin(us[i])
tmp<-rbind(Tiis,-Tiis)
spout<-SpherePowerND(df,Kv,Ks,CheckDists,Tiis,nperm)
ppout[i]<-c(spout$p.out)
}
power2markerdfs2[x,]<-ppout
}
matplot(t(power2markerdfs2),type="l")
