1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
medianaxis <- function(x, na.rm, ...) UseMethod("medianaxis")
medianaxis.default <- function(x, na.rm, ...) .NotYetImplemented()
#############################################################
#
# medianaxis.circular function
# Author: Claudio Agostinelli and Alessandro Gagliardi
# E-mail: claudio@unive.it
# Date: August, 3, 2011
# Version: 0.1
#
# Copyright (C) 2012 Claudio Agostinelli and Alessandro Gagliardi
#
#############################################################
medianaxis.circular <- function(x, na.rm=FALSE, ...) {
if (na.rm)
x <- x[!is.na(x)]
if (length(x)==0) {
warning("No observations (at least after removing missing values)")
return(NULL)
}
if (is.circular(x)) {
dc <- circularp(x)
} else {
dc <- list(type="angles", units="radians", template="none", modulo="asis", zero=0, rotation="counter")
}
x <- conversion.circular(x, units="radians", zero=0, rotation="counter")
attr(x, "class") <- attr(x, "circularp") <- NULL
circmedianaxis <- MedianAxisCircularRad(x)
circmedianaxis <- conversion.circular(circular(circmedianaxis$medianaxis), dc$units, dc$type, dc$template, dc$modulo, dc$zero, dc$rotation)
return(circmedianaxis)
}
MedianAxisCircularRad <- function(x) {
n <- length(x)
y <- c(x, x+pi)
y <- sort(y %% (2*pi))
if (odd <- as.logical(n%%2)) {
## odd number of observations
ismedianaxis <- rep(FALSE, n)
for (i in 1:length(x)) {
z <- MinusPiPlusPiRad((x-x[i])%%(2*pi))
pos <- sum(z > 0 & z != pi)
neg <- sum(z < 0 & z !=-pi)
zero <- sum(z==0)
atpi <- sum(z==pi | z==-pi)
### ismedianaxis[i] <- pos == (neg)
}
} else {
## even number of observations
eps <- min(diff(y), y[1]+(2*pi-y[2*n]))/4
}
}
|
