From 3cb37599b5c0a47c885a2f7a0fe86c681af072ad Mon Sep 17 00:00:00 2001 From: Devin Johnson Date: Fri, 21 Aug 2015 15:48:48 -0700 Subject: [PATCH] Added drift models and namespace corrections --- DESCRIPTION | 6 +- NAMESPACE | 1 - R/RcppExports.R | 40 +++++++--- R/crawl-package.R | 6 +- R/crwMLE.R | 2 +- R/crwN2ll.R | 16 ++-- R/crwPostIS.R | 14 +--- R/crwPredict.R | 43 +++++----- demo/harborSeal.R | 19 ++--- demo/harborSealSmooth.R | 20 ++--- demo/northernFurSeal.R | 28 +++---- man/crawl-package.Rd | 4 +- src/CTCRWN2LL.cpp | 73 +++-------------- src/CTCRWN2LL_DRIFT.cpp | 66 ++++++++++++++++ src/CTCRWPREDICT.cpp | 24 ++---- src/CTCRWPREDICT_DRIFT.cpp | 86 ++++++++++++++++++++ src/CTCRWSAMPLE.cpp | 33 ++------ src/CTCRWSAMPLE_DRIFT.cpp | 99 +++++++++++++++++++++++ src/RcppExports.cpp | 156 ++++++++++++++++++++++++++++++------- src/SMM_MATS.cpp | 75 ++++++++++++++++++ 20 files changed, 577 insertions(+), 234 deletions(-) create mode 100644 src/CTCRWN2LL_DRIFT.cpp create mode 100644 src/CTCRWPREDICT_DRIFT.cpp create mode 100644 src/CTCRWSAMPLE_DRIFT.cpp create mode 100644 src/SMM_MATS.cpp diff --git a/DESCRIPTION b/DESCRIPTION index 1cdbdec..12c4f4c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: crawl Type: Package Title: Fit Continuous-Time Correlated Random Walk Models to Animal Movement Data -Version: 1.913 -Date: 2015-08-20 +Version: 1.920 +Date: 2015-08-21 Author: Devin S. Johnson Maintainer: Devin S. Johnson Imports: @@ -10,9 +10,9 @@ Imports: sp, raster, Rcpp (>= 0.11.1), - RcppArmadillo (>= 0.4.200.0), methods LinkingTo: Rcpp, RcppArmadillo +Suggests: ggplot2 Description: The Correlated RAndom Walk Library of R functions was designed for fitting continuous-time correlated random walk (CTCRW) models with time indexed diff --git a/NAMESPACE b/NAMESPACE index 5bc5620..738df58 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,7 +19,6 @@ export(flatten) export(intToPOSIX) export(mergeTrackStop) import(Rcpp) -import(RcppArmadillo) import(mvtnorm) import(raster) import(sp) diff --git a/R/RcppExports.R b/R/RcppExports.R index 48443f8..85d111c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,23 +1,43 @@ # This file was generated by Rcpp::compileAttributes # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -makeT <- function(b, delta) { - .Call('crawl_makeT', PACKAGE = 'crawl', b, delta) +CTCRWNLL_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) { + .Call('crawl_CTCRWNLL_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) } -makeQ <- function(b, sig2, delta) { - .Call('crawl_makeQ', PACKAGE = 'crawl', b, sig2, delta) +CTCRWNLL <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) { + .Call('crawl_CTCRWNLL', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P) } -CTCRWNLL <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) { - .Call('crawl_CTCRWNLL', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P) +CTCRWPREDICT_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) { + .Call('crawl_CTCRWPREDICT_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) } -CTCRWPREDICT <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) { - .Call('crawl_CTCRWPREDICT', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P) +CTCRWPREDICT <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) { + .Call('crawl_CTCRWPREDICT', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P) } -CTCRWSAMPLE <- function(y, Hmat, Qmat, Tmat, noObs, active, a, P) { - .Call('crawl_CTCRWSAMPLE', PACKAGE = 'crawl', y, Hmat, Qmat, Tmat, noObs, active, a, P) +CTCRWSAMPLE_DRIFT <- function(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) { + .Call('crawl_CTCRWSAMPLE_DRIFT', PACKAGE = 'crawl', y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P) +} + +CTCRWSAMPLE <- function(y, Hmat, beta, sig2, delta, noObs, active, a, P) { + .Call('crawl_CTCRWSAMPLE', PACKAGE = 'crawl', y, Hmat, beta, sig2, delta, noObs, active, a, P) +} + +makeT <- function(b, delta, active) { + .Call('crawl_makeT', PACKAGE = 'crawl', b, delta, active) +} + +makeQ <- function(b, sig2, delta, active) { + .Call('crawl_makeQ', PACKAGE = 'crawl', b, sig2, delta, active) +} + +makeT_drift <- function(b, b_drift, delta, active) { + .Call('crawl_makeT_drift', PACKAGE = 'crawl', b, b_drift, delta, active) +} + +makeQ_drift <- function(b, b_drift, sig2, sig2_drift, delta, active) { + .Call('crawl_makeQ_drift', PACKAGE = 'crawl', b, b_drift, sig2, sig2_drift, delta, active) } diff --git a/R/crawl-package.R b/R/crawl-package.R index 8580989..af5a97d 100644 --- a/R/crawl-package.R +++ b/R/crawl-package.R @@ -10,8 +10,8 @@ #' \tabular{ll}{ #' Package: \tab crawl\cr #' Type: \tab Package\cr -#' Version: \tab 1.913\cr -#' Date: \tab August 20, 2015\cr +#' Version: \tab 1.920\cr +#' Date: \tab August 21, 2015\cr #' License: \tab file LICENSE \cr #' LazyLoad: \tab yes\cr #' } @@ -25,7 +25,7 @@ #' @references Johnson, D., J. London, M. -A. Lea, and J. Durban (2008) #' Continuous-time correlated random walk model for animal telemetry data. #' Ecology 89(5) 1208-1215. -#' @import Rcpp RcppArmadillo +#' @import Rcpp #' @importFrom graphics layout #' @importFrom methods as slot #' @importFrom stats approx filter model.frame model.matrix diff --git a/R/crwMLE.R b/R/crwMLE.R index caa293a..015961b 100644 --- a/R/crwMLE.R +++ b/R/crwMLE.R @@ -169,7 +169,7 @@ crwMLE = function(mov.model=~1, err.model=NULL, activity=NULL, drift=FALSE, initial.state, theta, fixPar, method="L-BFGS-B", control=NULL, constr=list(lower=-Inf, upper=Inf), prior=NULL, need.hess=TRUE, initialSANN=list(maxit=200), attempts=1) { - if(drift) stop("At this time drift models are not supported with this function. Use 'crwMLE' for now.\n") + #if(drift) stop("At this time drift models are not supported with this function. Use 'crwMLE' for now.\n") st <- Sys.time() # if (missing(Time.name) & !inherits(data,"STIDF")) stop("Argument 'Time.name' missing and NOT a spacetime STIDF object. Please specify") diff --git a/R/crwN2ll.R b/R/crwN2ll.R index 7f3dc33..352091f 100644 --- a/R/crwN2ll.R +++ b/R/crwN2ll.R @@ -71,21 +71,17 @@ crwN2ll = function(theta, fixPar, y, noObs, delta, a, active <- ifelse(b==Inf, 0, 1) b <- ifelse(b==Inf, 0, b) } else {active=rep(1,N)} - #if (driftMod) { - ### Change back for drift model - if(as.logical(FALSE)){ + if (driftMod) { theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1): (n.errX + n.errY + 2 * n.mov + 2)] b.drift <- exp(log(b) - log(1+exp(theta.drift[2]))) - sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) # - call.lik <- "crwdriftn2ll" + sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) + ll <- CTCRWNLL_DRIFT( y=as.matrix(y), Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a, P)$ll } else { - b.drift <- sig2.drift <- 0.0 - call.lik <- CTCRWNLL + ll <- CTCRWNLL( y=as.matrix(y), Hmat, b, sig2, delta, noObs, active, a, P)$ll } - movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod=FALSE) - - ll <- CTCRWNLL( y=as.matrix(y), Hmat, movMats[["Qmat"]], movMats[["Tmat"]], noObs, active, a, P)$ll + #movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod=FALSE) + if(is.null(prior)) return(-2 * ll) else return(-2 * (ll + prior(theta))) } diff --git a/R/crwPostIS.R b/R/crwPostIS.R index 3dd43aa..e797943 100644 --- a/R/crwPostIS.R +++ b/R/crwPostIS.R @@ -121,22 +121,16 @@ active <- ifelse(b==Inf, 0, 1) b <- ifelse(b==Inf, 0, b) } else active = rep(1,N) - #if (driftMod) { - ### Change back for drift model - if(as.logical(FALSE)){ + if (driftMod) { theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1): (n.errX + n.errY + 2 * n.mov + 2)] b.drift <- exp(log(b) - log(1+exp(theta.drift[2]))) - sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) # - call.lik <- "crwdriftn2ll" + sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) + out=CTCRWSAMPLE_DRIFT(y, Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a, P) } else { - b.drift <- sig2.drift <- 0.0 - #call.lik <- CTCRWNLL + out=CTCRWSAMPLE(y, Hmat, b, sig2, delta, noObs, active, a, P) } - movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod) - out=CTCRWSAMPLE(y, Hmat, movMats$Qmat, movMats$Tmat, noObs, active, a, P) - if(driftMod){ colnames(out$sim) <- apply(expand.grid(c("mu","theta","gamma"), c("x","y")), 1, paste, collapse=".") } else { diff --git a/R/crwPredict.R b/R/crwPredict.R index 5a1c235..032e077 100644 --- a/R/crwPredict.R +++ b/R/crwPredict.R @@ -115,29 +115,21 @@ crwPredict=function(object.crwFit, predTime=NULL, flat=TRUE, ...) theta.mov <- par[(n.errX + n.errY + 1):(n.errX + n.errY + 2 * n.mov)] sig2 <- exp(2 * (mov.mf %*% theta.mov[1:n.mov])) b <- exp(mov.mf %*% theta.mov[(n.mov + 1):(2 * n.mov)]) - #activity <- rep(1, N) if (!is.null(activity)) { theta.activ <- par[(n.errX + n.errY + 2 * n.mov + 1)] b <- b / ((activity) ^ exp(theta.activ)) active <- ifelse(b==Inf, 0, 1) b <- ifelse(b==Inf, 0, b) } else active = rep(1,N) - #if (driftMod) { - ### Change back for drift model - if(as.logical(FALSE)){ + if (driftMod) { theta.drift <- par[(n.errX + n.errY + 2 * n.mov + 1): (n.errX + n.errY + 2 * n.mov + 2)] b.drift <- exp(log(b) - log(1+exp(theta.drift[2]))) - sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) #rep(exp(2*theta.drift[1]), length(sig2)) # - call.lik <- "crwdriftn2ll" + sig2.drift <- exp(log(sig2) + 2 * theta.drift[1]) + out = CTCRWPREDICT_DRIFT(y, Hmat, b, b.drift, sig2, sig2.drift, delta, noObs, active, a, P) } else { - b.drift <- sig2.drift <- 0.0 - #call.lik <- CTCRWNLL + out=CTCRWPREDICT(y, Hmat, b, sig2, delta, noObs, active, a, P) } - movMats <- getQT(sig2, b, sig2.drift, b.drift, delta, driftMod) - - out=CTCRWPREDICT(y, Hmat, movMats$Qmat, movMats$Tmat, noObs, active, a, P) - pred <- data.frame(t(out$pred)) if (driftMod) { @@ -149,22 +141,23 @@ crwPredict=function(object.crwFit, predTime=NULL, flat=TRUE, ...) obsFit$outlier.chisq <- out$chisq obsFit$naive.p.val <- 1 - pchisq(obsFit$outlier.chisq, 2) if(getUseAvail){ - warning("'getUseAvail' not implemented yet in this version of 'crawl' contact Devin to fix this! ") -# idx <- data$locType=="p" -# movMatsPred <- getQT(sig2[idx], b[idx], sig2.drift[idx], b.drift[idx], delta=c(diff(data[idx,tn]),1), driftMod) -# TmatP <- movMatsPred$Tmat -# QmatP <- movMatsPred$Qmat -# avail <- t(sapply(1:(nrow(TmatP)-1), makeAvail, Tmat=TmatP, Qmat=QmatP, predx=predx[idx,], predy=predy[idx,], -# vary=vary[,,idx], varx=varx[,,idx], driftMod=driftMod, lonadj=lonAdjVals[idx])) -# avail <- cbind(data[idx,tn][-1], avail) -# colnames(avail) <- c(tn, "meanAvail.x", "meanAvail.y", "varAvail.x", "varAvail.y") -# use <- cbind(data[idx,tn], predx[idx,1], predy[idx,1], varx[1,1,idx], vary[1,1,idx])[-1,] -# colnames(use) <- c(tn, "meanUse.x", "meanUse.y", "varUse.x", "varUse.y") -# UseAvail.lst <- list(use=use, avail=avail) + warning("'getUseAvail' not implemented yet in this version of 'crawl' contact maintainer to fix this! ") + # idx <- data$locType=="p" + # movMatsPred <- getQT(sig2[idx], b[idx], sig2.drift[idx], b.drift[idx], delta=c(diff(data[idx,tn]),1), driftMod) + # TmatP <- movMatsPred$Tmat + # QmatP <- movMatsPred$Qmat + # avail <- t(sapply(1:(nrow(TmatP)-1), makeAvail, Tmat=TmatP, Qmat=QmatP, predx=predx[idx,], predy=predy[idx,], + # vary=vary[,,idx], varx=varx[,,idx], driftMod=driftMod, lonadj=lonAdjVals[idx])) + # avail <- cbind(data[idx,tn][-1], avail) + # colnames(avail) <- c(tn, "meanAvail.x", "meanAvail.y", "varAvail.x", "varAvail.y") + # use <- cbind(data[idx,tn], predx[idx,1], predy[idx,1], varx[1,1,idx], vary[1,1,idx])[-1,] + # colnames(use) <- c(tn, "meanUse.x", "meanUse.y", "varUse.x", "varUse.y") + # UseAvail.lst <- list(use=use, avail=avail) UseAvail.lst=NULL } else UseAvail.lst=NULL - speed = sqrt(apply(matrix(pred[,2:(2+driftMod)]), 1, sum)^2 + apply(matrix(pred[,(4+driftMod):(4+2*driftMod)]), 1, sum)^2) + speed = sqrt(apply(as.matrix(pred[,2:(2+driftMod)]), 1, sum)^2 + + apply(as.matrix(pred[,(4+driftMod):(4+2*driftMod)]), 1, sum)^2) out <- list(originalData=fillCols(data), alpha.hat=pred, V.hat=var, speed=speed, loglik=out$ll, useAvail=UseAvail.lst) if (flat) { diff --git a/demo/harborSeal.R b/demo/harborSeal.R index 832e62a..97c6461 100644 --- a/demo/harborSeal.R +++ b/demo/harborSeal.R @@ -7,17 +7,18 @@ harborSeal$Argos_loc_class = factor(harborSeal$Argos_loc_class, levels=c("3","2" ## Project data ## library(rgdal) -toProj = harborSeal[!is.na(harborSeal$latitude),] +toProj = harborSeal[!is.na(harborSeal$latitude),c("Time","latitude","longitude")] coordinates(toProj) = ~longitude+latitude proj4string(toProj) <- CRS("+proj=longlat") toProj <- spTransform(toProj, CRS("+init=epsg:3338")) toProj = as.data.frame(toProj) -harborSeal = merge(toProj, harborSeal, all=TRUE) +colnames(toProj)[2:3] = c("x","y") +harborSeal = merge(toProj, harborSeal, by="Time", all=TRUE) harborSeal = harborSeal[order(harborSeal$Time),] -initial.cpp = list( +initial = list( a=c(harborSeal$x[1],0,harborSeal$y[1],0), - P=diag(c(10000^2,54000^2,10000^2,5400^2)) + P=diag(c(10000^2,5400^2,10000^2,5400^2)) ) ##Fit model as given in Johnson et al. (2008) Ecology 89:1208-1215 @@ -34,10 +35,10 @@ set.seed(123) fit1 <- crwMLE( mov.model=~1, err.model=list(x=~Argos_loc_class-1), activity=~I(1-DryTime), data=harborSeal, coord=c("x","y"), Time.name="Time", - initial.state=initial.cpp, fixPar=fixPar, theta=c(rep(log(5000),3),log(3*3600), 3), + initial.state=initial, fixPar=fixPar, theta=c(rep(log(5000),3),log(3*3600), 0), constr=constr, control=list(maxit=2000, trace=1, REPORT=1), - initialSANN=list(maxit=200, trace=1, temp=20, tmax=20, REPORT=1) + initialSANN=list(maxit=200, trace=1, REPORT=1) ) print(fit1) @@ -47,9 +48,9 @@ p2=ggplot(aes(x=Time, y=mu.x), data=pred1) + geom_ribbon(aes(ymin=mu.x-2*se.mu.x geom_path(, col="red") + geom_point(aes(x=Time, y=x), col="blue", size=1) p3=ggplot(aes(x=Time, y=mu.y), data=pred1) + geom_ribbon(aes(ymin=mu.y-2*se.mu.y, ymax=mu.y+2*se.mu.y), fill="green", alpha=0.5) + geom_path(, col="red") + geom_point(aes(x=Time, y=y), col="blue", size=1) -print(p1) -print(p2) -print(p3) +suppressWarnings(print(p1)) +suppressWarnings(print(p2)) +suppressWarnings(print(p3)) # ggsave("map.pdf", p1) # ggsave("xaxis.pdf", p2, width=10, height=2) # ggsave("yaxis.pdf", p3, width=10, height=2) diff --git a/demo/harborSealSmooth.R b/demo/harborSealSmooth.R index 158e3d7..2c7e102 100644 --- a/demo/harborSealSmooth.R +++ b/demo/harborSealSmooth.R @@ -1,23 +1,23 @@ -require(crawl) -require(splines) +# library(crawl) +library(ggplot2) +library(splines) +library(rgdal) data(harborSeal) head(harborSeal) harborSeal$Argos_loc_class = factor(harborSeal$Argos_loc_class, levels=c("3","2","1","0","A","B")) ## Project data ## -library(sp) -library(rgdal) -toProj = harborSeal[!is.na(harborSeal$latitude),] +toProj = harborSeal[!is.na(harborSeal$latitude),c("Time","latitude","longitude")] coordinates(toProj) = ~longitude+latitude proj4string(toProj) <- CRS("+proj=longlat") toProj <- spTransform(toProj, CRS("+init=epsg:3338")) toProj = as.data.frame(toProj) -harborSeal = merge(toProj, harborSeal, all=TRUE) +colnames(toProj)[2:3] = c("x","y") +harborSeal = merge(toProj, harborSeal, by="Time", all=TRUE) harborSeal = harborSeal[order(harborSeal$Time),] -#harborSeal$Time = harborSeal$Time*3600 -initial.cpp = list( +initial = list( a=c(harborSeal$x[1],0,harborSeal$y[1],0), P=diag(c(10000^2,54000^2,10000^2,5400^2)) ) @@ -43,8 +43,8 @@ fit1 <- crwMLE( constr=constr, theta = theta.start, prior=prior, - control=list(maxit=2000, trace=1, REPORT=1), - initialSANN=list(maxit=10000, temp=100, tmax=100, trace=1, REPORT=1) + control=list(maxit=2000, trace=1, REPORT=1)#, + #initialSANN=list(maxit=10000, temp=100, tmax=100, trace=1, REPORT=1) ) print(fit1) diff --git a/demo/northernFurSeal.R b/demo/northernFurSeal.R index a1544b8..72d7a4f 100644 --- a/demo/northernFurSeal.R +++ b/demo/northernFurSeal.R @@ -10,34 +10,34 @@ proj4string(northernFurSeal) <- CRS("+proj=longlat") northernFurSeal <- spTransform(northernFurSeal, CRS("+proj=aea +lat_1=30 +lat_2=70 +lat_0=52 +lon_0=-170 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")) initial = list( - a=c(coordinates(northernFurSeal)[1,1],0,coordinates(northernFurSeal)[1,2],0), - P=diag(c(10000^2,54000^2,10000^2,5400^2)) + a=c(coordinates(northernFurSeal)[1,1],0,0,coordinates(northernFurSeal)[1,2],0,0), + P=diag(c(10000^2,5400^2,5400^2,10000^2,5400^2,5400^2)) ) -fixPar = c(log(250), log(500), log(1500), rep(NA,3), NA) -displayPar( mov.model=~1, err.model=list(x=~Argos_loc_class-1), +fixPar = c(log(250), log(500), log(1500), rep(NA,6)) +displayPar( mov.model=~1, err.model=list(x=~Argos_loc_class-1), drift = TRUE, data=northernFurSeal,fixPar=fixPar) constr=list( - lower=c(rep(log(1500),2), rep(-Inf,2)), - upper=rep(Inf,4) + lower=c(rep(log(1500),2), rep(-Inf,4)), + upper=rep(Inf,6) ) -ln.prior = function(theta){-abs(theta[4]-3)/0.5} +ln.prior = function(theta){-abs(theta[4]+3)/0.5} #set.seed(123) fit1 <- crwMLE( - mov.model=~1, err.model=list(x=~Argos_loc_class-1), + mov.model=~1, err.model=list(x=~Argos_loc_class-1), drift=TRUE, data=northernFurSeal, Time.name="Time", - initial.state=initial, fixPar=fixPar, constr=constr, prior=ln.prior, - control=list(maxit=30, trace=1, REPORT=1), - initialSANN=list(maxit=200, trace=1, REPORT=1) + initial.state=initial, fixPar=fixPar, constr=constr,# prior=ln.prior, + control=list(maxit=300, trace=1, REPORT=1), + initialSANN=list(maxit=200, temp=100, tmax=20, trace=1, REPORT=1) ) fit1 ##Make hourly location predictions predTime <- seq(ceiling(min(northernFurSeal$Time)), floor(max(northernFurSeal$Time)), 1) predObj <- crwPredict(object.crwFit=fit1, predTime, speedEst=TRUE, flat=TRUE) head(predObj) -crwPredictPlot(predObj, "map") +crwPredictPlot(predObj, "map", asp=TRUE) ##Create simulation object with 100 parameter draws set.seed(123) @@ -54,13 +54,13 @@ round(100/(1+(sd(w)/mean(w))^2)) jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) -crwPredictPlot(predObj, 'map') +crwPredictPlot(predObj, 'map', asp=TRUE) ## Sample 20 tracks from posterior predictive distribution iter <- 20 cols <- jet.colors(iter) for(i in 1:iter){ - samp <- crwPostIS(simObj) + samp <- crwPostIS(simObj, fullPost = FALSE) lines(samp$alpha.sim[,'mu.x'], samp$alpha.sim[,'mu.y'],col=cols[i]) } diff --git a/man/crawl-package.Rd b/man/crawl-package.Rd index 6daf24a..03498f4 100644 --- a/man/crawl-package.Rd +++ b/man/crawl-package.Rd @@ -16,8 +16,8 @@ version of the continuous-time staochistic movement process. \tabular{ll}{ Package: \tab crawl\cr Type: \tab Package\cr -Version: \tab 1.913\cr -Date: \tab August 20, 2015\cr +Version: \tab 1.920\cr +Date: \tab August 21, 2015\cr License: \tab file LICENSE \cr LazyLoad: \tab yes\cr } diff --git a/src/CTCRWN2LL.cpp b/src/CTCRWN2LL.cpp index d670c7d..48cda7a 100644 --- a/src/CTCRWN2LL.cpp +++ b/src/CTCRWN2LL.cpp @@ -4,61 +4,23 @@ using namespace Rcpp; using namespace arma; -// Below is a simple example of exporting a C++ function to R. You can -// source this function into an R session using the Rcpp::sourceCpp -// function (or via the Source button on the editor toolbar) - -// For more on using Rcpp click the Help button on the editor toolbar - -//const double log2pi = std::log(2.0 * M_PI); -//double ln_dmvnrm(const arma::vec& x, const arma::vec& mean, const arma::mat& sigma) { -// int xdim = x.n_elem; -// double out; -// arma::mat rooti = arma::trans(arma::inv(trimatu(arma::chol(sigma)))); -// double rootisum = arma::sum(log(rooti.diag())); -// double constants = -(static_cast(xdim)/2.0) * log2pi; -// arma::vec z = rooti * ( x - mean) ; -// out = - 0.5 * arma::sum(z%z) + rootisum; -// return(out); -//} - -// [[Rcpp::export]] -arma::mat makeT(const double& b, const double& delta){ - arma::mat T(4,4, fill::zeros); - T(0,0) = 1; T(2,2) = 1; - T(0,1) = exp(R::pexp(delta,1/b,1,1) - log(b)); - T(1,1) = exp(-b*delta); - T(2,3) = exp(R::pexp(delta,1/b,1,1) - log(b)); - T(3,3) = exp(-b*delta); - return T; -} - -// [[Rcpp::export]] -arma::mat makeQ(const double& b, const double& sig2, const double& delta){ - arma::mat Q(4,4, fill::zeros); - Q(0,0) = sig2*(delta - 2*exp(R::pexp(delta,1/b,1,1)-log(b)) + exp(R::pexp(delta,1/(2*b),1,1)-log(2*b))); - Q(1,1) = sig2*exp(log(b) + R::pexp(delta,1/(2*b),1,1))/2; - Q(0,1) = sig2*(1-2*exp(-b*delta)+exp(-2*b*delta))/2; - Q(1,0) = Q(0,1); - Q.submat(2,2,3,3) = Q.submat(0,0,1,1); - return Q; -} - +// Function prototypes +arma::mat makeQ(const double& b, const double& sig2, const double& delta, const double& active); +arma::mat makeT(const double& b, const double& delta, const double& active); // [[Rcpp::export]] Rcpp::List CTCRWNLL(const arma::mat& y, const arma::mat& Hmat, -const arma::mat& Qmat, const arma::mat& Tmat, -const arma::vec& noObs,const arma::vec& active, const arma::colvec& a, -const arma::mat& P) + const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta, + const arma::vec& noObs,const arma::vec& active, const arma::colvec& a, + const arma::mat& P) { // Define fixed matrices int N = y.n_rows; double detF; arma::mat Z(2,4, fill::zeros); Z(0,0) = 1; Z(1,2) = 1; - arma::mat T(4,4, fill::zeros); - T(0,0) = 1; T(2,2) = 1; - arma::mat Q(4,4, fill::zeros); + arma::mat T(4,4); + arma::mat Q(4,4); arma::mat F(2,2, fill::zeros); arma::mat H(2,2, fill::zeros); arma::mat K(4,2, fill::zeros); @@ -71,23 +33,8 @@ const arma::mat& P) double ll=0; //Begin Kalman filter for(int i=0; i0; j--){ + if(noObs(j-1)==1 || F.slice(j-1)(0,0)*F.slice(j-1)(1,1)==0){ + r = L.slice(j-1).t() * r; + N = L.slice(j-1).t() * N * L.slice(j-1); + } else{ + u.col(j-1) = solve(F.slice(j-1),v.col(j-1))-K.slice(j-1).t()*r; + M.slice(j-1) = F.slice(j-1).i() + K.slice(j-1).t()*N*K.slice(j-1); + chisq(j-1) = dot(u.col(j-1),solve(M.slice(j-1),u.col(j-1))); + jk.col(j-1) = y.row(j-1).t() - solve(M.slice(j-1),u.col(j-1)); + r = Z.t()*solve(F.slice(j-1),v.col(j-1)) + L.slice(j-1).t() * r; + N = Z.t() * solve(F.slice(j-1),Z) + L.slice(j-1).t()*N*L.slice(j-1); + } + pred.col(j-1) = aest.col(j-1) + Pest.slice(j-1)*r; + predVar.slice(j-1) = Pest.slice(j-1) - Pest.slice(j-1)*N*Pest.slice(j-1); + } + return Rcpp::List::create( + Rcpp::Named("ll") = ll, + Rcpp::Named("pred") = pred, Rcpp::Named("predVar")=predVar, + Rcpp::Named("chisq")=chisq, Rcpp::Named("predObs")=jk + ); +} diff --git a/src/CTCRWSAMPLE.cpp b/src/CTCRWSAMPLE.cpp index a98956a..6d3ecf1 100644 --- a/src/CTCRWSAMPLE.cpp +++ b/src/CTCRWSAMPLE.cpp @@ -4,23 +4,16 @@ using namespace Rcpp; using namespace arma; - -arma::vec armaNorm(int n){ - NumericVector x = Rcpp::rnorm(n,0,1); - arma::vec out(x.begin(), x.size(), false); - return out; -} - -arma::vec mvn(const arma::vec& mu, const arma::mat& V){ - arma::mat out = mu + chol(V).t()*armaNorm(mu.n_elem); - return out; -} +// Function prototypes +arma::mat makeQ(const double& beta, const double& sig2, const double& delta, const double& active); +arma::mat makeT(const double& beta, const double& delta, const double& active); +arma::vec mvn(const arma::vec& mu, const arma::mat& V); // [[Rcpp::export]] Rcpp::List CTCRWSAMPLE(const arma::mat& y, const arma::mat& Hmat, -const arma::mat& Qmat, const arma::mat& Tmat, +const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta, const arma::vec& noObs,const arma::vec& active, const arma::colvec& a, const arma::mat& P) { @@ -52,23 +45,11 @@ const arma::mat& P) double ll=0; //Forward filter and simulation for(int i=0; i0; j--){ + if(noObs(j-1)==1 || F.slice(j-1)(0,0)*F.slice(j-1)(1,1)==0){ + r = L.slice(j-1).t() * r; + r_plus = L.slice(j-1).t() * r_plus; + } else{ + r = Z.t()*solve(F.slice(j-1),v.col(j-1)) + L.slice(j-1).t() * r; + r_plus = Z.t()*solve(F.slice(j-1),v_plus.col(j-1)) + L.slice(j-1).t() * r_plus; + } + alpha_hat.col(j-1) += P_hat.slice(j-1)*r; + alpha_plus_hat.col(j-1) += P_hat.slice(j-1)*r_plus; + } + sim = (alpha_hat.cols(0,I-1) - (alpha_plus.cols(0,I-1)-alpha_plus_hat.cols(0,I-1))).t(); + return Rcpp::List::create( + Rcpp::Named("ll") = ll, + Rcpp::Named("sim") = sim + ); +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8f26e6e..ea621b2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,82 +6,180 @@ using namespace Rcpp; -// makeT -arma::mat makeT(const double& b, const double& delta); -RcppExport SEXP crawl_makeT(SEXP bSEXP, SEXP deltaSEXP) { +// CTCRWNLL_DRIFT +Rcpp::List CTCRWNLL_DRIFT(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& beta_drift, const arma::vec& sig2, const arma::vec& sig2_drift, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWNLL_DRIFT(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP beta_driftSEXP, SEXP sig2SEXP, SEXP sig2_driftSEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject __result; Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< const double& >::type b(bSEXP); - Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); - __result = Rcpp::wrap(makeT(b, delta)); + Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta_drift(beta_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2_drift(sig2_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); + Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); + __result = Rcpp::wrap(CTCRWNLL_DRIFT(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)); return __result; END_RCPP } -// makeQ -arma::mat makeQ(const double& b, const double& sig2, const double& delta); -RcppExport SEXP crawl_makeQ(SEXP bSEXP, SEXP sig2SEXP, SEXP deltaSEXP) { +// CTCRWNLL +Rcpp::List CTCRWNLL(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWNLL(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP sig2SEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject __result; Rcpp::RNGScope __rngScope; - Rcpp::traits::input_parameter< const double& >::type b(bSEXP); - Rcpp::traits::input_parameter< const double& >::type sig2(sig2SEXP); - Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); - __result = Rcpp::wrap(makeQ(b, sig2, delta)); + Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); + Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); + __result = Rcpp::wrap(CTCRWNLL(y, Hmat, beta, sig2, delta, noObs, active, a, P)); return __result; END_RCPP } -// CTCRWNLL -Rcpp::List CTCRWNLL(const arma::mat& y, const arma::mat& Hmat, const arma::mat& Qmat, const arma::mat& Tmat, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); -RcppExport SEXP crawl_CTCRWNLL(SEXP ySEXP, SEXP HmatSEXP, SEXP QmatSEXP, SEXP TmatSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { +// CTCRWPREDICT_DRIFT +Rcpp::List CTCRWPREDICT_DRIFT(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& beta_drift, const arma::vec& sig2, const arma::vec& sig2_drift, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWPREDICT_DRIFT(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP beta_driftSEXP, SEXP sig2SEXP, SEXP sig2_driftSEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject __result; Rcpp::RNGScope __rngScope; Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Qmat(QmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Tmat(TmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta_drift(beta_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2_drift(sig2_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); - __result = Rcpp::wrap(CTCRWNLL(y, Hmat, Qmat, Tmat, noObs, active, a, P)); + __result = Rcpp::wrap(CTCRWPREDICT_DRIFT(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)); return __result; END_RCPP } // CTCRWPREDICT -Rcpp::List CTCRWPREDICT(const arma::mat& y, const arma::mat& Hmat, const arma::mat& Qmat, const arma::mat& Tmat, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); -RcppExport SEXP crawl_CTCRWPREDICT(SEXP ySEXP, SEXP HmatSEXP, SEXP QmatSEXP, SEXP TmatSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { +Rcpp::List CTCRWPREDICT(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWPREDICT(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP sig2SEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject __result; Rcpp::RNGScope __rngScope; Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Qmat(QmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Tmat(TmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); - __result = Rcpp::wrap(CTCRWPREDICT(y, Hmat, Qmat, Tmat, noObs, active, a, P)); + __result = Rcpp::wrap(CTCRWPREDICT(y, Hmat, beta, sig2, delta, noObs, active, a, P)); + return __result; +END_RCPP +} +// CTCRWSAMPLE_DRIFT +Rcpp::List CTCRWSAMPLE_DRIFT(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& beta_drift, const arma::vec& sig2, const arma::vec& sig2_drift, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWSAMPLE_DRIFT(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP beta_driftSEXP, SEXP sig2SEXP, SEXP sig2_driftSEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { +BEGIN_RCPP + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta_drift(beta_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2_drift(sig2_driftSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); + Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); + __result = Rcpp::wrap(CTCRWSAMPLE_DRIFT(y, Hmat, beta, beta_drift, sig2, sig2_drift, delta, noObs, active, a, P)); return __result; END_RCPP } // CTCRWSAMPLE -Rcpp::List CTCRWSAMPLE(const arma::mat& y, const arma::mat& Hmat, const arma::mat& Qmat, const arma::mat& Tmat, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); -RcppExport SEXP crawl_CTCRWSAMPLE(SEXP ySEXP, SEXP HmatSEXP, SEXP QmatSEXP, SEXP TmatSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { +Rcpp::List CTCRWSAMPLE(const arma::mat& y, const arma::mat& Hmat, const arma::vec& beta, const arma::vec& sig2, const arma::vec& delta, const arma::vec& noObs, const arma::vec& active, const arma::colvec& a, const arma::mat& P); +RcppExport SEXP crawl_CTCRWSAMPLE(SEXP ySEXP, SEXP HmatSEXP, SEXP betaSEXP, SEXP sig2SEXP, SEXP deltaSEXP, SEXP noObsSEXP, SEXP activeSEXP, SEXP aSEXP, SEXP PSEXP) { BEGIN_RCPP Rcpp::RObject __result; Rcpp::RNGScope __rngScope; Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Hmat(HmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Qmat(QmatSEXP); - Rcpp::traits::input_parameter< const arma::mat& >::type Tmat(TmatSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type beta(betaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type delta(deltaSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type noObs(noObsSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type active(activeSEXP); Rcpp::traits::input_parameter< const arma::colvec& >::type a(aSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type P(PSEXP); - __result = Rcpp::wrap(CTCRWSAMPLE(y, Hmat, Qmat, Tmat, noObs, active, a, P)); + __result = Rcpp::wrap(CTCRWSAMPLE(y, Hmat, beta, sig2, delta, noObs, active, a, P)); + return __result; +END_RCPP +} +// makeT +arma::mat makeT(const double& b, const double& delta, const double& active); +RcppExport SEXP crawl_makeT(SEXP bSEXP, SEXP deltaSEXP, SEXP activeSEXP) { +BEGIN_RCPP + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< const double& >::type b(bSEXP); + Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const double& >::type active(activeSEXP); + __result = Rcpp::wrap(makeT(b, delta, active)); + return __result; +END_RCPP +} +// makeQ +arma::mat makeQ(const double& b, const double& sig2, const double& delta, const double& active); +RcppExport SEXP crawl_makeQ(SEXP bSEXP, SEXP sig2SEXP, SEXP deltaSEXP, SEXP activeSEXP) { +BEGIN_RCPP + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< const double& >::type b(bSEXP); + Rcpp::traits::input_parameter< const double& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const double& >::type active(activeSEXP); + __result = Rcpp::wrap(makeQ(b, sig2, delta, active)); + return __result; +END_RCPP +} +// makeT_drift +arma::mat makeT_drift(const double& b, const double& b_drift, const double& delta, const double& active); +RcppExport SEXP crawl_makeT_drift(SEXP bSEXP, SEXP b_driftSEXP, SEXP deltaSEXP, SEXP activeSEXP) { +BEGIN_RCPP + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< const double& >::type b(bSEXP); + Rcpp::traits::input_parameter< const double& >::type b_drift(b_driftSEXP); + Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const double& >::type active(activeSEXP); + __result = Rcpp::wrap(makeT_drift(b, b_drift, delta, active)); + return __result; +END_RCPP +} +// makeQ_drift +arma::mat makeQ_drift(const double& b, const double& b_drift, const double& sig2, const double& sig2_drift, const double& delta, const double& active); +RcppExport SEXP crawl_makeQ_drift(SEXP bSEXP, SEXP b_driftSEXP, SEXP sig2SEXP, SEXP sig2_driftSEXP, SEXP deltaSEXP, SEXP activeSEXP) { +BEGIN_RCPP + Rcpp::RObject __result; + Rcpp::RNGScope __rngScope; + Rcpp::traits::input_parameter< const double& >::type b(bSEXP); + Rcpp::traits::input_parameter< const double& >::type b_drift(b_driftSEXP); + Rcpp::traits::input_parameter< const double& >::type sig2(sig2SEXP); + Rcpp::traits::input_parameter< const double& >::type sig2_drift(sig2_driftSEXP); + Rcpp::traits::input_parameter< const double& >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< const double& >::type active(activeSEXP); + __result = Rcpp::wrap(makeQ_drift(b, b_drift, sig2, sig2_drift, delta, active)); return __result; END_RCPP } diff --git a/src/SMM_MATS.cpp b/src/SMM_MATS.cpp new file mode 100644 index 0000000..90d3307 --- /dev/null +++ b/src/SMM_MATS.cpp @@ -0,0 +1,75 @@ +// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*- +#include "RcppArmadillo.h" +// [[Rcpp::depends(RcppArmadillo)]] +using namespace Rcpp; +using namespace arma; + +// Random normal draws for posterior sampling +arma::vec armaNorm(int n){ + NumericVector x = Rcpp::rnorm(n,0,1); + arma::vec out(x.begin(), x.size(), false); + return out; +} +arma::vec mvn(const arma::vec& mu, const arma::mat& V){ + arma::mat out = mu + chol(V).t()*armaNorm(mu.n_elem); + return out; +} + +// [[Rcpp::export]] +arma::mat makeT(const double& b, const double& delta, const double& active){ + arma::mat T(4,4, fill::zeros); + T(0,0) = 1; + T(2,2) = 1; + if(active > 0){ + T(0,1) = exp(R::pexp(delta,1/b,1,1) - log(b)); + T(1,1) = exp(-b*delta); + T(2,3) = exp(R::pexp(delta,1/b,1,1) - log(b)); + T(3,3) = exp(-b*delta); + } + return T; +} + +// [[Rcpp::export]] +arma::mat makeQ(const double& b, const double& sig2, const double& delta, const double& active){ + arma::mat Q(4,4, fill::zeros); + if(active > 0){ + Q(0,0) = sig2*(delta - 2*exp(R::pexp(delta,1/b,1,1)-log(b)) + exp(R::pexp(delta,1/(2*b),1,1)-log(2*b))); + Q(1,1) = sig2*exp(log(b) + R::pexp(delta,1/(2*b),1,1))/2; + Q(0,1) = sig2*(1-2*exp(-b*delta)+exp(-2*b*delta))/2; + Q(1,0) = Q(0,1); + Q.submat(2,2,3,3) = Q.submat(0,0,1,1); + } + return Q; +} + +// [[Rcpp::export]] +arma::mat makeT_drift(const double& b, const double& b_drift, const double& delta, const double& active){ + arma::mat T(6,6, fill::zeros); + T(0,0) = 1; + if(active > 0){ + T(0,1) = exp(R::pexp(delta,1/b,1,1) - log(b)); + T(0,2) = exp(R::pexp(delta,1/b_drift,1,1) - log(b_drift)); + T(1,1) = exp(-b*delta); + T(2,2) = exp(-b_drift*delta); + } + T.submat(3,3,5,5) = T.submat(0,0,2,2); + return T; +} + +// [[Rcpp::export]] +arma::mat makeQ_drift(const double& b, const double& b_drift, const double& sig2, const double& sig2_drift, + const double& delta, const double& active){ + arma::mat Q(6,6, fill::zeros); + if(active > 0){ + Q(0,0) = sig2*(delta - 2*exp(R::pexp(delta,1/b,1,1)-log(b)) + exp(R::pexp(delta,1/(2*b),1,1)-log(2*b))) + + sig2_drift*(delta - 2*exp(R::pexp(delta,1/b_drift,1,1)-log(b_drift)) + exp(R::pexp(delta,1/(2*b_drift),1,1)-log(2*b_drift))); + Q(1,1) = sig2*exp(log(b) + R::pexp(delta,1/(2*b),1,1))/2; + Q(2,2) = sig2_drift*exp(log(b_drift) + R::pexp(delta,1/(2*b_drift),1,1))/2; + Q(0,1) = sig2*(1-2*exp(-b*delta)+exp(-2*b*delta))/2; + Q(1,0) = Q(0,1); + Q(0,2) = sig2_drift*(1-2*exp(-b_drift*delta)+exp(-2*b_drift*delta))/2; + Q(2,0) = Q(0,2); + Q.submat(3,3,5,5) = Q.submat(0,0,2,2); + } + return Q; +}