# R/density.circular.R In circular: Circular Statistics

#### Documented in density.circularlines.density.circularplot.density.circularprint.density.circular

```#############################################################
#                                                           #
#   density.circular function                               #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   date: February, 14, 2013                                #
#   Copyright (C) 2013 Claudio Agostinelli                  #
#                                                           #
#   Version 0.3-1                                           #
#                                                           #
#############################################################

density.circular <- function(x, z=NULL, bw, adjust = 1, type = c("K", "L"), kernel= c("vonmises", "wrappednormal"), na.rm = FALSE, from=circular(0), to=circular(2*pi), n=512, K=NULL, min.k=10, control.circular=list(), ...) {

name <- deparse(substitute(x))
data <- x

if (!is.numeric(from))
stop("argument 'from' must be numeric")
if (!is.numeric(to))
stop("argument 'to' must be numeric")
if (!is.finite(from))
stop("non-finite `from'")
if (!is.finite(to))
stop("non-finite `to'")
if (!is.numeric(n))
stop("argument 'n' must be numeric")
n <- round(n)
if (n <=0)
stop("argument 'n' must be integer and positive")
if (!is.numeric(x))
stop("argument 'x' must be numeric")
if (!is.null(z) && is.circular(z)) {
datacircularp <- circularp(z)
} else if (is.circular(x))
datacircularp <- circularp(x)
else {
datacircularp <- list(type="angles", units="radians", template="none", modulo="asis", zero=0, rotation="counter")
}
dc <- control.circular
if (is.null(dc\$type))
dc\$type <- datacircularp\$type
if (is.null(dc\$units))
dc\$units <- datacircularp\$units
if (is.null(dc\$template))
dc\$template <- datacircularp\$template
if (is.null(dc\$modulo))
dc\$modulo <- datacircularp\$modulo
if (is.null(dc\$zero))
dc\$zero <- datacircularp\$zero
if (is.null(dc\$rotation))
dc\$rotation <- datacircularp\$rotation
if (dc\$modulo=="pi")
stop("The function does not work yet for modulo='pi'")

x <- conversion.circular(x, units="radians", zero=0, rotation="counter")
attr(x, "class") <- attr(x, "circularp") <- NULL
from <- conversion.circular(from, units="radians", zero=0, rotation="counter")
attr(from, "class") <- attr(from, "circularp") <- NULL
to <- conversion.circular(to, units="radians", zero=0, rotation="counter")
attr(to, "class") <- attr(to, "circularp") <- NULL

kernel <- match.arg(kernel)

x <- as.vector(x)
x.na <- is.na(x)
if (any(x.na)) {
if (na.rm)
x <- x[!x.na]
else stop("x contains missing values")
}

x.finite <- is.finite(x)
if (any(!x.finite)) {
x <- x[x.finite]
}
nx <- length(x)

if (is.null(z)) {
z <- circular(seq(from=from, to=to, length=n))
} else {
if (!is.numeric(z))
stop("argument 'z' must be numeric")
namez <- deparse(substitute(z))
z.na <- is.na(z)
if (any(z.na)) {
if (na.rm) {
z <- z[!z.na]
} else {
stop("z contains missing values")
}
}
z.finite <- is.finite(z)
if (any(!z.finite)) {
z <- z[z.finite]
}
}
zz <- conversion.circular(z, dc\$units, dc\$type, dc\$template, dc\$modulo, dc\$zero, dc\$rotation)
z <- conversion.circular(z, units="radians", zero=0, rotation="counter")
attr(z, "class") <- attr(z, "circularp") <- NULL
z <- as.vector(z)

if (!is.numeric(bw))
stop("argument 'bw' and 'adjust' must be numeric")
if (!is.finite(bw))
stop("non-finite `bw'")
if (bw <= 0)
stop("`bw' is not positive.")

y <- DensityCircularRad(x=x, z=z, bw=bw, kernel=kernel, K=K, min.k=min.k)

structure(list(data = data, x = zz, y = y, bw = bw, n = nx, kernel=kernel, call = match.call(), data.name=name, has.na = FALSE), class = "density.circular")
}

DensityCircularRad <- function(x, z, bw, kernel, K=NULL, min.k=10) {
nx <- length(x)
if (kernel=="vonmises") {
y <- sapply(z, DvonmisesRad, mu=x, kappa=bw, log=FALSE)
} else if (kernel=="wrappednormal") {
rho <- exp(-bw^2/2)
y <- sapply(z, DwrappednormalRad, mu=x, rho=rho, K=K, min.k=min.k)
} else {
stop("other kernels not implemented yet")
}
y <- apply(y, 2, sum)/nx
return(y)
}

#############################################################
#                                                           #
#   plot.density.circular function                          #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 22, 2011                                     #
#   Copyright (C) 2011 Claudio Agostinelli                  #
#                                                           #
#   Version 0.5-12                                           #
#                                                           #
#############################################################

plot.density.circular <- function(x, main = NULL, sub=NULL, xlab = NULL, ylab ="Density circular", type = "l", zero.line = TRUE, points.plot=FALSE, points.col=1, points.pch=1, points.cex=1, plot.type = c("circle", "line"), axes=TRUE, ticks=FALSE, bins=NULL, offset=1, shrink=1, tcl=0.025,  tcl.text=0.125, sep=0.025, tol = 0.04, digits=2, cex=1, uin=NULL, xlim=NULL, ylim=NULL, join=FALSE, nosort=FALSE, units=NULL, template=NULL, zero=NULL, rotation=NULL, control.circle=circle.control(), ...) {
xcircularp <- attr(x\$x, "circularp")
if (is.null(xcircularp))
stop("the component 'x' of the object must be of class circular")
##   type <- xcircularp\$type
modulo <- xcircularp\$modulo
if (is.null(units))
units <- xcircularp\$units
if (is.null(template))
template <- xcircularp\$template
if (template=="geographics" | template=="clock24") {
zero <- pi/2
rotation <- "clock"
} else if (template=="clock12") {
zero <- pi/2
rotation <- "clock"
modulo <- "pi"
} else {
if (is.null(zero))
zero <- xcircularp\$zero
if (is.null(rotation))
rotation <- xcircularp\$rotation
}
next.points <- 0
if (template=="clock12") {
x\$x <- 2*x\$x
x\$data <- 2*x\$data
}
attr(x\$x, "class") <- attr(x\$x, "circularp") <- NULL
attr(x\$data, "class") <- attr(x\$data, "circularp") <- NULL

plot.type <- match.arg(plot.type)

if (is.null(xlab))
xlab <- paste("N =", x\$n, "  Bandwidth =", formatC(x\$bw), " Unit =", units)
if (is.null(main))
main <- deparse(x\$call)
#### as scatter plot
if (plot.type == "line") {
if (units=='degrees') {
x\$x <- x\$x/pi*180
x\$data <- x\$data/pi*180
}
if (units=='hours') {
x\$x <- x\$x/pi*12
x\$data <- x\$data/pi*12
}
if (is.null(xlim))
xlim <- range(c(x\$x, x\$data))
if (is.null(ylim)) {
ylim <- range(x\$y)
if (points.plot)
ylim[1] <- ylim[1]-0.04*points.cex
}

xorder <- order(x\$x)
x\$x <- x\$x[xorder]
x\$y <- x\$y[xorder]

plot.default(x, type = type, xlim=xlim, ylim=ylim, main=main, xlab=xlab, ylab=ylab, axes=axes, ...)
if (zero.line)
abline(h = 0, lwd = 0.2, col = "gray")
if (points.plot)
points(x\$data, rep(ylim[1]+0.02*points.cex, length(x\$data)), col=points.col, pch=points.pch, cex=points.cex)
return(NULL)
} else {
#### as circular plot
if (is.null(xlim))
xlim <- c(-1, 1)
if (is.null(ylim))
ylim <- c(-1, 1)
if (is.null(bins)) {
bins <- NROW(x)
} else {
bins <- round(bins)
if (bins<=0)
stop("'bins' must be non negative")
}

CirclePlotRad(xlim=xlim, ylim=ylim, uin=uin, shrink=shrink, tol=tol, main=main, sub=sub, xlab=xlab, ylab=ylab, control.circle=control.circle)

if (!is.logical(ticks))
stop("ticks must be logical")

if (axes) {
axis.circular(at=NULL, labels=NULL, units=units, template=template, modulo=modulo, zero=zero, rotation=rotation, tick=ticks, cex=cex, tcl=tcl, tcl.text=tcl.text, digits=digits)
}
if (axes==FALSE & ticks) {
at <- circular((0:bins)/bins*2*pi, zero=zero, rotation=rotation)
ticks.circular(at, tcl=tcl)
}

if (rotation=="clock")
x\$x <- -x\$x
x\$x <- x\$x+zero
x\$x <- x\$x%%(2*pi)
ll <- LinesCircularRad(x=x\$x, y=x\$y, join=join, nosort=nosort, offset=offset, shrink=shrink, ...)

if (points.plot) {
next.points <- sep
if (rotation=="clock")
x\$data <- -x\$data
x\$data <- x\$data+zero
x\$data <- x\$data%%(2*pi)
PointsCircularRad(x\$data, bins, FALSE, points.col, points.pch, 1, 1, 0, sep, next.points, shrink, points.cex)
}
return(invisible(list(x=ll\$x, y=ll\$y, zero=zero, rotation=rotation, next.points=next.points)))
}
}

#############################################################
#                                                           #
#   lines.density.circular function                         #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 5, 2011                                     #
#   Copyright (C) 2010 Claudio Agostinelli                  #
#                                                           #
#   Version 0.5-1                                           #
#                                                           #
#############################################################

lines.density.circular <- function(x, type = "l", zero.line = TRUE, points.plot=FALSE, points.col=1, points.pch=1, points.cex=1, plot.type = c("circle", "line"), bins=NULL, offset=1, shrink=1, tcl=0.025, sep=0.025, join=TRUE, nosort=FALSE, plot.info=NULL, zero=NULL, rotation=NULL, ...) {
xcircularp <- attr(x\$x, "circularp")
if (is.null(xcircularp))
stop("the component 'x' of the object must be of class circular")
##   type <- xcircularp\$type
modulo <- xcircularp\$modulo
units <- xcircularp\$units
template <- xcircularp\$template
if (is.null(plot.info)) {
if (is.null(zero)) {
if (template=="geographics" | template=="clock24" | template=="clock12")
zero <- pi/2
else
zero <- xcircularp\$zero
}
if (is.null(rotation)) {
if (template=="geographics" | template=="clock24" | template=="clock12")
rotation <- "clock"
else
rotation <- xcircularp\$rotation
}
next.points <- 0
} else {
zero <- plot.info\$zero
rotation <- plot.info\$rotation
next.points <- plot.info\$next.points
}

if (template=="clock12") {
x\$x <- 2*x\$x
x\$data <- 2*x\$data
}

attr(x\$x, "circularp") <- attr(x\$x, "class") <- NULL
attr(x\$data, "circularp") <- attr(x\$data, "class") <- NULL

ll <- list()

plot.type <- match.arg(plot.type)
if (is.null(bins)) {
bins <- NROW(x)
} else {
bins <- round(bins)
if (bins<=0)
stop("bins must be non negative")
}

if (plot.type == "line") {
if (units=='degrees') {
x\$x <- x\$x/pi*180
x\$data <- x\$data/pi*180
}
if (units=='hours') {
x\$x <- x\$x/pi*12
x\$data <- x\$data/pi*12
}
xorder <- order(x\$x)
x\$x <- x\$x[xorder]
x\$y <- x\$y[xorder]
lines.default(x, type = type, ...)
if (zero.line)
abline(h = 0, lwd = 0.2, col = "gray")
if (points.plot)
points.default(x\$data, rep(min(x\$y)-0.02*points.cex, length(x\$data)), col=points.col, pch=points.pch)
} else {
if (rotation=="clock")
x\$x <- -x\$x
x\$x <- x\$x+zero
x\$x <- x\$x%%(2*pi)
ll <- LinesCircularRad(x=x\$x, y=x\$y, join=join, nosort=nosort, offset=offset, shrink=shrink, ...)
if (points.plot) {
if (rotation=="clock")
x\$data <- -x\$data
x\$data <- x\$data+zero
x\$data <- x\$data%%(2*pi)
next.points <- next.points+sep
PointsCircularRad(x\$data, bins, FALSE, points.col, points.pch, 1, 1, 0, sep, next.points, shrink, points.cex)
}
}

return(invisible(list(x=ll\$x, y=ll\$y, zero=zero, rotation=rotation, next.points=next.points)))
}

#############################################################
#                                                           #
#   print.density.circular function                         #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 23, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#                                                           #
#############################################################

print.density.circular <- function(x, digits=NULL, ...)
{
cat("\nCall:\n\t",deparse(x\$call),
"\n\nData: ",x\$data.name," (",x\$n," obs.);",
"\tBandwidth 'bw' = ",formatC(x\$bw,digits=digits), "\n\n",sep="")
print(summary(as.data.frame(x[c("x","y")])), digits=digits, ...)
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.