SpherePower2D
フェノタイプ情報のベクトルと、ジェノタイプ情報の行列があるとす
SpherePower2D(Kv,Ks,CheckDists,Tiis,nperm)
SpherePower2D<-
function (Kv, Ks, CheckDists, Tiis, nperm = 1000)
{
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))
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, 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
PowerOut[ik, ] <- PowOut
}
list(df = df, Kv = Kv, Ks = Ks, CheckDists = CheckDists,
Tiis = Tiis, nperm = nperm, p.out = PowerOut)
}
\name{SpherePower2D}
\alias{SpherePower2D}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
SpherePower2D
}
\description{
}
\usage{
SpherePower2D(Kv,Ks,CheckDists,Tiis,nperm)
}
\arguments{
}
%- maybe also 'usage' for other objects documented here.
\details{
%% ~~ If necessary, more details than the description above ~~
}
\value{
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%% ~~who you are~~
}
\note{
%% ~~further notes~~
}
%% ~Make other sections like Warning with \section{Warning }{....} ~
\seealso{
}
\examples{
}
