# R/change.point.R In circular: Circular Statistics

#### Documented in change.point

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

#############################################################
#                                                           #
#   change.point function                                   #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: May, 26, 2006                                     #
#   Copyright (C) 2006 Claudio Agostinelli                  #
#                                                           #
#   Version 0.2                                             #
#############################################################

change.point <- function(x) {
x <- conversion.circular(x, units="radians", zero=0, rotation="counter", modulo="2pi")
attr(x, "class") <- attr(x, "circularp") <-  NULL
}

phi <- function(x) {
arg <- A1inv(x)
if (besselI(x=arg, nu=0, expon.scaled = FALSE) != Inf)
result <- x * A1inv(x) - log(besselI(x=arg, nu=0, expon.scaled = FALSE))
else
result <- x * A1inv(x) - (arg + log(1/sqrt(2 * pi * arg) * (1 + 1/(8 * arg) + 9/(128 * arg^2) + 225/(1024 * arg^3))))
return(result)
}
n <- length(x)
R1 <- 1:n
R2 <- 1:n
V <- 1:n
for (k in 1:(n - 1)) {
R2[k] <- RhoCircularRad(x[(k + 1):n]) * (n - k)
if (k >= 2 & k <= (n - 2)) {
V[k] <- k/n * phi(R1[k]/k) + (n - k)/n * phi(R2[k]/(n - k))
}
}
R1[n] <- rho * n
R2[n] <- 0
R.diff <- R1 + R2 - rho * n
rmax <- max(R.diff)
rave <- mean(R.diff)
k.r <- (1:n)[R.diff == max(R.diff)]
V <- V[2:(n - 2)]
if (n > 3) {
tmax <- max(V)
tave <- mean(V)
k.t <- (1:(n - 3))[V == max(V)] + 1
} else
stop("Sample size must be at least 4")
return(list(n=n, rho=rho, rmax=rmax, k.r=k.r, rave=rave, tmax=tmax, k.t=k.t, tave=tave))
}
```

## 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.