Nothing
### Statistical tests ##########################################################
### Tests of multivariate normality ############################################
##' @title Normality Test based on Distribution of Squared Mahalanobis Distances
##' @param x (n, d) data matrix
##' @param type the type of test to be used
##' @param dist distribution to check
##' @param ... additional arguments passed to the underlying type
##' @return htest object as returned by ad.test() or ks.test()
##' @author Marius Hofert
maha2_test <- function(x, type = c("ad.test", "ks.test"), dist = c("chi2", "beta"),
...)
{
if(!is.matrix(x))
x <- as.matrix(x)
type <- match.arg(type)
dist <- match.arg(dist)
## Build squared Mahalanobis distances
Xbar <- colMeans(x)
S <- cov(x)
D2 <- mahalanobis(x, center = Xbar, cov = S)
## Address different types
d <- ncol(x)
switch(type,
"ad.test" = {
switch(dist,
"chi2" = {
ad.test(D2, distr.fun = pchisq, df = d, ...)
},
"beta" = {
n <- nrow(x)
ad.test(D2 * n * (n-1)^2, distr.fun = pbeta,
shape1 = d/2, shape2 = (n-d-1)/2, ...)
},
stop("Wrong 'dist'"))
},
"ks.test" = {
switch(dist,
"chi2" = {
ks.test(D2, y = "pchisq", df = d, ...)
},
"beta" = {
n <- nrow(x)
ks.test(D2 * n * (n-1)^2, y = "pbeta",
shape1 = d/2, shape2 = (n-d-1)/2, ...)
},
stop("Wrong 'dist'"))
},
stop("Wrong 'type'"))
}
##' @title Mardia's Test of Multivariate Normality
##' @param x (n, d) data matrix
##' @param type the type of test to be used ("skewness" is based on Mahalanobis angles,
##' "kurtosis" on Mahalanobis distances).
##' @param method method for computing the Mahalanobis angles or distances
##' @return htest object
##' @author Marius Hofert
mardia_test <- function(x, type = c("kurtosis", "skewness"), method = c("direct", "chol"))
{
## Basics
if(!is.matrix(x))
x <- as.matrix(x)
n <- nrow(x)
d <- ncol(x)
type <- match.arg(type)
method <- match.arg(method)
## Main
switch(type,
"kurtosis" = {
## Build squared Mahalanobis distances
D2 <- mahalanobis(x, center = colMeans(x), cov = cov(x))
## Compute the test results
k <- mean(D2^2)
T <- (k - d * (d + 2)) / sqrt(8 * d * (d + 2) / n)
p.val <- pnorm(abs(T), lower.tail = FALSE)
alternative <- "two-sided"
strng <- paste0("Mardia's kurtosis test (computed with method = '",method,"')")
},
"skewness" = {
## Build Mahalanobis angles
X.centered <- sweep(x, MARGIN = 2, STATS = colMeans(x)) # X_i - bar(X)
S <- cov(x) # sample covariance matrix S
D <- switch(method,
"direct" = {
S.inv <- solve(S) # S^{-1}; as in mahalanobis()
X.centered %*% S.inv %*% t(X.centered) # (n, d) %*% (d, d) %*% (d, n)
},
"chol" = { # similar to QRM::MardiaTest but more efficient
R <- chol(S)
R.inv <- solve(R)
Z <- X.centered %*% R.inv # (n, d)-matrix
Z %*% t(Z) # (n, d) %*% (d, n)
},
stop("Wrong 'method'"))
## Compute the test statistic results
b <- mean(D^3)
T <- (n/6) * b
p.val <- pchisq(T, d * (d + 1) * (d + 2) / 6, lower.tail = FALSE)
alternative <- "one-sided"
strng <- paste0("Mardia's skewness test (computed with method = '",method,"')")
},
stop("Wrong 'type'"))
## Build and return htest object
structure(class = "htest",
list(statistic = c(statistic = T), # have to give it a name so that it shows up in output
p.value = p.val, # shows up with 'p-value'
alternative = alternative,
method = strng,
data.name = deparse(substitute(x)))) # name of the provided data object
}
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.