# inst/statsPureR/statsPureR.R In PoweR: Computation of Power and Level Tables for Hypothesis Tests

```stat2 <- function(x){ # Anderson-Darling
x <- sort(x)
n <- length(x)
logp1 <- pnorm((x - mean(x)) / sd(x), log.p = TRUE)
logp2 <- pnorm(-(x - mean(x)) / sd(x), log.p = TRUE)
h <- (2 * seq(1:n) - 1) * (logp1 + rev(logp2))
A <- -n - mean(h)
AA <- (1 + 0.75 / n + 2.25 / n ^ 2) * A
if (AA < 0.2) {
pval <- 1 - exp(-13.436 + 101.14 * AA - 223.73 * AA ^ 2)
}
else if (AA < 0.34) {
pval <- 1 - exp(-8.318 + 42.796 * AA - 59.938 * AA ^ 2)
}
else if (AA < 0.6) {
pval <- exp(0.9177 - 4.279 * AA - 1.38 * AA ^ 2)
}
else if (AA < 10) {
pval <- exp(1.2937 - 5.709 * AA + 0.0186 * AA ^ 2)
}
else pval <- 3.7e-24
res <- list(stat = A, p.value = pval)
return(res)
}

stat3 <- function(x){ # Zhang Zc
n <- length(x)
phi <- pnorm(sort((x - mean(x)) / sd(x)))
res <- sum((log((1 / phi - 1) / ((n - 0.5) / ((1:n) - 0.75) - 1))) ^ 2)
return(res)
}

stat4 <- function(x) { # Zhang Za
n <- length(x)
phi <- pnorm(sort((x - mean(x)) / sd(x)))
res <- -sum(((log(phi)) / (n - (1:n) + 0.5)) + ((log(1.0 - phi)) / ((1:n) - 0.5)))
return(10.0 * res - 32.0) # See page 711 in their paper
}

stat6 <- function(x) { # D'Agostino-Pearson. Credits to: fBasics package, dagoTest()
n <- length(x)
meanX <- mean(x)
s <- sqrt(mean((x - meanX) ^ 2))
a3 <- mean((x - meanX) ^ 3)/s ^ 3
a4 <- mean((x - meanX) ^ 4) / s ^ 4
SD3 <- sqrt(6 * (n - 2) / ((n + 1) * (n + 3)))
SD4 <- sqrt(24 * (n - 2) * (n - 3) * n / ((n + 1) ^ 2 * (n + 3) *
(n + 5)))
U3 <- a3 / SD3
U4 <- (a4 - 3 + 6 / (n + 1)) / SD4
b <- (3 * (n ^ 2 + 27 * n - 70) * (n + 1) * (n + 3)) / ((n - 2) *
(n + 5) * (n + 7) * (n + 9))
W2 <- sqrt(2 * (b - 1)) - 1
delta <- 1 / sqrt(log(sqrt(W2)))
a <- sqrt(2 / (W2 - 1))
Z3 <- delta * log((U3 / a) + sqrt((U3 / a) ^ 2 + 1))
B <- (6 * (n * n - 5 * n + 2) / ((n + 7) * (n + 9))) * sqrt((6 *
(n + 3) * (n + 5)) / (n * (n - 2) * (n - 3)))
A <- 6 + (8 / B) * ((2 / B) + sqrt(1 + 4 / (B ^ 2)))
jm <- sqrt(2 / (9 * A))
pos <- ((1 - 2 / A) / (1 + U4 * sqrt(2 / (A - 4)))) ^ (1 / 3)
Z4 <- (1 - 2 / (9 * A) - pos) / jm
omni <- Z3 ^ 2 + Z4 ^ 2
pomni <- 1 - pchisq(omni, 2)
res <- list(stat = omni, p.value = pomni)
return(res)
}

stat7 <- function(x){ # Jarque-Bera. Credits to: tseries package, jarque.bera.test()
n <- length(x)
m1 <- sum(x) / n
m2 <- sum((x - m1) ^ 2) / n
m3 <- sum((x - m1) ^ 3) / n
m4 <- sum((x - m1) ^ 4) / n
b1 <- (m3 / m2 ^ (3 / 2)) ^ 2
b2 <- (m4 / m2 ^ 2)
STATISTIC <- n * b1 / 6 + n * (b2 - 3) ^ 2 / 24
PVAL <- 1 - pchisq(STATISTIC, df = 2)
res <- list(stat = STATISTIC, p.value = PVAL)
return(res)
}

stat8 <- function(x){ # Doornik-Hansen. Credits to: normwhn.test package, normality.test1()
n <- length(x)
m1 <- sum(x) / n
v.matrix <- 1 / sqrt(cov.wt(as.matrix(x))\$cov)
xhatprime.matrix <- t(x - m1)
trprime <- t(v.matrix %*% xhatprime.matrix)
v1 <- mean(trprime)
v2 <- (n ^ (-1)) * sum((trprime - v1) ^ 2)
v3 <- (n ^ (-1)) * sum((trprime - v1) ^ 3)
v4 <- (n ^ (-1)) * sum((trprime - v1) ^ 4)
rtb1 <- v3 / (v2 ^ (3 / 2))
b2 <- v4 / v2 ^ 2
beta <- (3 * (n ^ 2 + 27 * n - 70) * (n + 1) * (n + 3)) / ((n -
2) * (n + 5) * (n + 7) * (n + 9))
w2 <- (-1) + sqrt(2 * (beta - 1))
delta <- 1 / sqrt(log(sqrt(w2)))
f <- (w2 - 1) / 2
g <- (n + 1) * (n + 3) / (6 * (n - 2))
h <- sqrt(f * g)
y <- rtb1 * h
z1 <- delta * log(y + sqrt(y ^ 2 + 1))
del <- ((n - 3) * (n + 1) * (n ^ 2 + 15 * n - 4))
aye <- ((n - 2) * (n + 5) * (n + 7) * (n ^ 2 + 27 * n - 70)) / (6 *
del)
cee <- ((n - 7) * (n + 5) * (n + 7) * (n ^ 2 + 2 * n - 5)) / (6 *
del)
alp <- aye + ((rtb1 ^ 2) * cee)
kap <- ((n + 5) * (n + 7) * (n ^ 3 + 37 * n ^ 2 + 11 * n - 313)) / (12 *
del)
chi <- (b2 - 1 - rtb1 ^ 2) * (2 * kap)
chi <- abs(chi)
z2 <- (((chi / (2 * alp)) ^ (1 / 3)) - 1 + (1 / ((9 * alp)))) * sqrt(9 *
alp)
stat <- z1 ^ 2 + z2 ^ 2
pval <- 1 - pchisq(stat, 2)
res <- list(stat = stat, p.value = pval)
return(res)
}

stat17 <-function (x) { # Bonett-Seier. Credits to moments package, bonett.test()
x <- sort(x)
n <- length(x)
rho <- sqrt(sum((x - mean(x)) ^ 2) / n)
tau <- sum(abs(x - mean(x))) / n
omega <- 13.29 * (log(rho) - log(tau))
z <- sqrt(n + 2) * (omega - 3) / 3.54
pval <- 2* pnorm(z, lower.tail = FALSE)
if (pval > 1)
pval <- 2 - pval
res <- list(stat = z, p.value = pval)
return(res)
}

stat26 <- function(x){ # Chen-Shapiro
n <- length(x)
m <- qnorm(((1:n) - 0.375) / (n + 0.25))
xo <- sort(x)
QH <- sum(diff(xo) / diff(m)) / (n - 1) / sd(x);
res <- sqrt(n) * (1 - QH)
return(res)
}

stat29 <- function(x){ # BCMR
n <- length(x)
xo <- sort(x)
dt <- 0.000001
k <- 1:n
sigma <- 0
for (i in k){
u <- seq((k[i] - 1) / n + dt / 2, k[i] / n - dt / 2, dt)
sigma <- sigma +  sum(qnorm(u)) * xo[i] * dt
}
stat <- 1 - sigma ^ 2 / (var(x) * (n - 1) / n)
return(stat)
}

stat30 <- function(x) { # Coin. Credits to http://digilander.libero.it/polynomial/B3test.R
normorder <- function(n){
dt <- 0.01
t <- seq(-10, 10, dt)
int <- function(r) {
rinv <- n - r + 1
logh <- lgamma(n + 1) - lgamma(rinv) - lgamma(n - rinv + 1) +
(rinv - 1) * pnorm(t, low = FALSE, log.p = TRUE) + (n - rinv) * pnorm(t, low = TRUE, log.p = TRUE) +
dnorm(t, log = TRUE)
h <- t * exp(logh)
return(sum(h) * dt)
}
res <- apply(cbind(1:n), 1, int)
return(res)
}
n <- length(x)
a <- normorder(n)
x <- sort(x)
z <- (x - mean(x)) / sd(x)
mod <- lm(z ~ 0 + a + I(a ^ 3))
stat <- mod\$coef[2] ^ 2
return(stat = stat)
}

stat36 <- function(x) { # Xapd
n <- length(x)
euler <- - digamma(1)
z <- (x - mean(x)) / sqrt(var(x) * (n - 1) / n)
z <- z[! (z == 0)]
B2 <- sum(z ^ 2 * sign(z)) / n
K2 <- sum(z ^ 2 * log(abs(z))) / n
var.B2 <- (1 / n) * (3.0 - 8.0 / pi) * (1.0 - 1.9 / n)
Z.B2 <-  B2 / sqrt(var.B2)
esp.net.K2 <- ((2.0 - log(2.0) - euler) / 2.0) ^ (1.0 / 3.0) *
(1 - 1.026 / n)
var.net.K2 <- (1.0 / n) * (1.0 / 72.0) *
((2.0 - log(2.0) - euler) / 2.0) ^ (-4.0 / 3.0) *
(3.0 * pi ^ 2 - 28.0) * (1.0 - 2.25 / n ^ 0.8) ;
Z.net.K2 <- ((K2 - B2 ^ 2) ^ (1.0 / 3.0) - esp.net.K2) /
sqrt(var.net.K2)
Xapd.stat <- Z.B2 ^ 2 + Z.net.K2 ^ 2 ;
p.value <- pchisq(Xapd.stat, 2, lower.tail = FALSE) ;
res <- list(test = "2nd-power skewness and kurtosis omnibus normality test Xapd",
B2 = B2, K2 = K2, Z.B2 = Z.B2, Z.net.K2 = Z.net.K2,
Xapd.stat = Xapd.stat, p.value = p.value)
return(res)
}

stat37 <- function(x) { # Zepd
n <- length(x)
euler <- -digamma(1)
z <- (x - mean(x)) / sqrt(var(x) * (n - 1) / n)
z <- z[! (z == 0)]
K2 <- sum(z ^ 2 * log(abs(z))) / n
alpha.n <- -0.06 + 2.1 / n ^ 0.67
K2.box.cox <- ((2.0 * K2) ^ alpha.n - 1.0) / alpha.n
esp.K2.box.cox <- - ((2.0 - log(2.0) - euler) ^ (-0.06) - 1.0) / 0.06 -
1.32 / n ^ 0.95
var.K2.box.cox <- (1.0 / n) * ((2.0 - log(2.0) - euler) ^ (-2.12) *
(3.0 * pi ^ 2 - 28.0) / 2.0 - 3.78 / n ^ 0.733)
stat.Z.K2 <- (K2.box.cox - esp.K2.box.cox) / sqrt(var.K2.box.cox)
p.value <- 2.0 * pnorm(abs(stat.Z.K2), lower.tail = FALSE)
res <- list(test = "2nd-power kurtosis directional normality test Zepd",
stat.Z.K2 = stat.Z.K2, p.value = p.value)
return(res)
}
```

## Try the PoweR package in your browser

Any scripts or data that you put into this service are public.

PoweR documentation built on Aug. 28, 2018, 5:05 p.m.