Generate a Multistate Correlated Random Walk path (preliminary version)

mscrwpath.R

#' @return a two-column matrix of turning angles `"ta"` and step lengths `"len"`.
#' @param Nstates the number of movement states. 
#' @param tp a vector of transition probabilities for `Nstates == 2` or a
#' square transition matrix `[Nstates, Nstates]`. Ignored if `Nstates == 1` 
#' @param mu,rho wrapped Cauchy distribution parameters for the turning angles. Should have a length of `Nstates`.
#' @param wscale,wshape Weibull distribution parameters for step length
# TODO: starting position and initial direction. Currently all are set to 0.
mscrwpath <-
function(N, Nstates = 2, tp = c(.05, .1),
		 mu = c(0, 0), rho = c(.5, .99),
		 wscale = c(10, 5), wshape = c(2, 1)) {
	
	if(is.matrix(tp)) {
		stopifnot(nrow(tp) == ncol(tp))
		stopifnot(tp >= 0 & tp <=1)
		TP <- tp
	} else {
		TP <- matrix(NA, ncol = Nstates, nrow = Nstates)
		if(Nstates == 1) {
		   TP[1L] <- 1
		} else if(Nstates == 2L) {
			TP[c(3L,2L,1L,4L)] <- c(tp[c(1L, 2L)], 1 - tp[c(1L, 2L)])
		}
	}
	# TODO: check lengths of mu,rho,wscale,wshape
	ang <- numeric(N)
	state <- integer(N) 
	state[1L] <- 1L
	for(i in seq.int(2L, N))
		state[i] <- sample.int(Nstates, 1L, prob = TP[state[i - 1L], ])

	rval <- matrix(NA, N, 2L, dimnames = list(NULL, c("ta", "len")))
	rval[, 1L] <- rwrpcauchy(N, location = mu[state], rho = rho[state])
	rval[, 2L] <- rweibull(N, scale = wscale[state], shape = wshape[state])
	
	structure(rval, state = state, Nstates = Nstates, TP = TP,
		mu = mu, rho = rho, wscale = wscale, wshape = wshape,
		angle0 = 0, pos0 = c(0, 0),
		class = c("crw.path", "path", "matrix"))
}

#' Converts the `"path"` object to `x` and `y` coordinates.
coordinates.path <-
function(x) {
	cbind(
		x = attr(x, "pos0")[1L] +
			cumsum(cos(attr(x, "angle0") + cumsum(x[, 1L])) * x[, 2L]),
		y = attr(x, "pos0")[2L] +
			cumsum(sin(attr(x, "angle0") + cumsum(x[, 1L])) * x[, 2L])
		)
}

plot.path <-
function(x, asp = 1, type = "l", ...) {
	a0 <- attr(x, "angle0")
	plot(data.frame(x = cumsum(cos(a0 + cumsum(x[, 1L])) * x[, 2L]),
			y = cumsum(sin(a0 + cumsum(x[, 1L])) * x[, 2L])),
		asp = asp, type = type, ...)
}


rwrpcauchy  <-
function (n, location = 0, rho = exp(-1)) {
	if(n %% length(rho) > 0) warning("'n' is not a multiple of length of 'rho'")
	if(n %% length(location) > 0) warning("'n' is not a multiple of length of 'location'")
	rho <- rep(rho, length.out = n)
    location <- rep(location, length.out = n)
	
	rval <- numeric(n)
	i0 <- rho == 0
	if((m <- sum(i0)) > 0) rval[i0] <- runif(m, 0, 2 * pi)
	i1 <- rho == 1
	if((m <- sum(i1)) > 0) rval[i1] <- location[i1]
	
	ok <- !i0 & !i1 
    rval[ok] <- rcauchy(sum(ok), location[ok], -log(rho[ok])) %% (2 * pi)
	rval
}