# R/watson.test.R In circular: Circular Statistics

#### Documented in print.watson.testwatson.test

```#############################################################
#                                                           #
#       Original Splus: Ulric Lund                          #
#       E-mail: ulund@calpoly.edu                           #
#                                                           #
#############################################################

#############################################################
#                                                           #
#   watson.test function                                    #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio.agostinelli@unive.it                    #
#   Date: April, 18, 2016                                   #
#   Version: 0.3-2                                          #
#                                                           #
#   Copyright (C) 2016 Claudio Agostinelli                  #
#                                                           #
#############################################################

watson.test <- function(x, alpha = 0, dist = c("uniform", "vonmises")) {
# Handling missing values
x <- na.omit(x)
if (length(x)==0) {
warning("No observations (at least after removing missing values)")
return(NULL)
}
dist <- match.arg(dist)
x <- conversion.circular(x, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(x, "circularp") <- attr(x, "class") <- NULL
if (!any(c(0, 0.01, 0.025, 0.05, 0.1)==alpha))
stop("'alpha' must be one of the following values: 0, 0.01, 0.025, 0.05, 0.10")
result <- WatsonTestRad(x, dist)
result\$call <- match.call()
result\$n <- length(x)
result\$alpha <- alpha
result\$dist <- dist
class(result) <-"watson.test"
return(result)
}

WatsonTestRad <- function(x, dist) {
n <- length(x)
if (dist == "uniform") {
u <- sort(x)/(2 * pi)
u.bar <- mean.default(u)
i <- 1:n
sum.terms <- (u - u.bar - (2 * i - 1)/(2 * n) + 0.5)^2
u2 <- sum(sum.terms) + 1/(12 * n)
u2 <- (u2 - 0.1/n + 0.1/(n^2)) * (1 + 0.8/n)
result <- list(statistic=u2, row=NA)
} else {
res <- MlevonmisesRad(x, bias=FALSE)
mu.hat <- res[1]
kappa.hat <- res[4]
x <- (x - mu.hat) %% (2 * pi)
x <- matrix(x, ncol = 1)
z <- apply(x, 1, PvonmisesRad, mu=0, kappa=kappa.hat, tol=1e-020)
z <- sort(z)
z.bar <- mean.default(z)
i <- 1:n
sum.terms <- (z - (2 * i - 1)/(2 * n))^2
Value <- sum(sum.terms) - n * (z.bar - 0.5)^2 + 1/(12 * n)
if (kappa.hat < 0.25)
row <- 1
else if (kappa.hat < 0.75)
row <- 2
else if (kappa.hat < 1.25)
row <- 3
else if (kappa.hat < 1.75)
row <- 4
else if (kappa.hat < 3)
row <- 5
else if (kappa.hat < 5)
row <- 6
else row <- 7
result <- list(statistic=Value, row=row)
}
return(result)
}

#############################################################
#                                                           #
#   print.watson.test function                              #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.watson.test <- function(x, digits=4, ...) {
dist <- x\$dist
n <- x\$n
alpha <- x\$alpha

if (dist == "uniform") {
u2 <- x\$statistic
cat("\n", "      Watson's Test for Circular Uniformity", "\n", "\n")
crits <- c(99, 0.267, 0.221, 0.187, 0.152)
if (n < 8) {
warning("Total Sample Size < 8:  Results may not be valid", "\n", "\n")
}
cat("Test Statistic:", round(u2, digits=digits), "\n")
if (alpha == 0) {
if (u2 > 0.267)
cat("P-value < 0.01", "\n", "\n")
else if (u2 > 0.221)
cat("0.01 < P-value < 0.025", "\n", "\n")
else if (u2 > 0.187)
cat("0.025 < P-value < 0.05", "\n", "\n")
else if (u2 > 0.152)
cat("0.05 < P-value < 0.10", "\n", "\n")
else cat("P-value > 0.10", "\n", "\n")
} else {
index <- (1:5)[alpha == c(0, 0.01, 0.025, 0.05, 0.1)]
Critical <- crits[index]
if (u2 > Critical)
Reject <- "Reject Null Hypothesis"
else Reject <- "Do Not Reject Null Hypothesis"
cat("Level", alpha, "Critical Value:", round(Critical, digits=digits), "\n")
cat(Reject, "\n\n")
}

} else if (dist=="vonmises") {
Value <- x\$statistic
row <- x\$row
cat("\n", "      Watson's Test for the von Mises Distribution \n\n")
u2.crits <- cbind(c(0, 0.5, 1, 1.5, 2, 4, 100), c(0.052, 0.056, 0.066, 0.077, 0.084, 0.093, 0.096), c(0.061, 0.066, 0.079, 0.092, 0.101, 0.113, 0.117), c(0.081, 0.09, 0.11, 0.128, 0.142, 0.158, 0.164))

if (alpha != 0) {
if (alpha == 0.1)
col <- 2
else if (alpha == 0.05)
col <- 3
else if (alpha == 0.01)
col <- 4
Critical <- u2.crits[row, col]
if (Value > Critical)
Reject <- "Reject Null Hypothesis"
else Reject <- "Do Not Reject Null Hypothesis"
cat("Test Statistic:", round(Value, digits=digits), "\n")
cat("Level", alpha, "Critical Value:", round(Critical, digits=digits), "\n")
cat(Reject, "\n\n")
} else {
cat("Test Statistic:", round(Value, digits=digits), "\n")
if (Value < u2.crits[row, 2])
cat("P-value > 0.10", "\n", "\n")
else if ((Value >= u2.crits[row, 2]) && (Value < u2.crits[row, 3]))
cat("0.05 < P-value < 0.10", "\n", "\n")
else if ((Value >= u2.crits[row, 3]) && (Value < u2.crits[row, 4]))
cat("0.01 < P-value < 0.05", "\n", "\n")
else cat("P-value < 0.01", "\n", "\n")
}

}
invisible(x)
}
```

## Try the circular package in your browser

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

circular documentation built on May 1, 2019, 7:57 p.m.