Nothing
hcfcirc.perm <- function(u1, u2, rads = TRUE, B = 999) {
if ( !rads ) {
u1 <- u1 * pi/180
u2 <- u2 * pi/180
}
u <- c(u1, u2)
ina <- c( rep(1, length(u1) ), rep(2, length(u2) ) )
ni <- tabulate(ina)
n <- sum(ni)
x1 <- cos(u)
x2 <- sin(u)
Ci <- Rfast::group(x1, ina)
Si <- Rfast::group(x2, ina)
Ri <- sqrt(Ci^2 + Si^2)
V <- sum(Ri)
C <- sum(Ci)
S <- sum(Si)
R <- sqrt(C^2 + S^2)
mu <- atan(S/C) + pi * ( C < 0 )
con <- sum( cos(u - mu) )
k1 <- (1.28 - 0.53 * R^2/n^2) * tan(0.5 * pi * R/n)
if (k1 < 710) {
der <- con - n * besselI(k1, 1, expon.scaled = TRUE)/besselI(k1, 0, expon.scaled = TRUE)
a <- besselI(k1, 0)^2/2 + besselI(k1, 2) * besselI(k1, 0)/2 - besselI(k1, 1)^2
der2 <- n * a/besselI(k1, 0)^2
k2 <- k1 + der/der2
while (abs(k1 - k2) > 1e-7) {
k1 <- k2
der <- con - n * besselI(k1, 1, expon.scaled = TRUE)/besselI(k1, 0, expon.scaled = TRUE)
a <- besselI(k1, 0)^2/2 + besselI(k1, 2) * besselI(k1, 0)/2 - besselI(k1, 1)^2
der2 <- n * a/besselI(k1, 0)^2
k2 <- k1 + der/der2
}
} else k2 <- k1
kapa <- k2
if (kapa > 2) {
Ft <- (n - 2) * (V - R)/(n - V)
} else if (kapa < 2 & kapa > 1) {
Ft <- (1 + 3/(8 * kapa) ) * (n - 2) * (V - R) / (n - V)
} else Ft <- NA
p.value <- NA
if ( !is.na(Ft) ) {
pft <- numeric(B)
for (i in 1:B) {
id <- Rfast2::Sample.int(n, n)
Ci <- Rfast::group(x1[id], ina)
Si <- Rfast::group(x2[id], ina)
Ri <- sqrt(Ci^2 + Si^2)
V <- sum(Ri)
if (kapa > 2) {
pft[i] <- (n - 2) * (V - R)/(n - V)
} else if (kapa < 2 & kapa > 1) {
pft[i] <- (1 + 3/(8 * kapa)) * (n - 2) * (V - R) / (n - V)
} else pft[i] <- NA
}
p.value <- ( sum(pft > Ft) + 1 ) / (B + 1)
}
statistic <- Ft ; names(statistic) <- "hcf test statistic"
parameter <- "NA" ; names(parameter) <- "df"
alternative <- "The 2 circular means differ"
method <- "Permutation ANOVA for 2 circular means using the high concentration test"
data.name <- c("data ", " groups")
result <- list( statistic = statistic, parameter = parameter, p.value = p.value,
alternative = alternative, method = method, data.name = data.name )
class(result) <- "htest"
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.