From 49da34a4086c12242b3a1af7f9152e099b41da1c Mon Sep 17 00:00:00 2001 From: David Lawrence Miller Date: Wed, 15 Jul 2015 20:48:56 +0100 Subject: [PATCH] minor speed-up to logistic code when distance is a covariate --- NEWS | 1 + R/g0.R | 2 +- R/integratelogisticdup.R | 26 +++++++++++++++++++++++--- R/logisticdupbyx.R | 10 +++++++--- R/logisticdupbyx_fast.R | 24 ++++++++++++++++++++++++ R/predict.io.fi.R | 10 +++++----- man/logisticdupbyx_fast.Rd | 30 ++++++++++++++++++++++++++++++ 7 files changed, 91 insertions(+), 12 deletions(-) create mode 100644 R/logisticdupbyx_fast.R create mode 100644 man/logisticdupbyx_fast.Rd diff --git a/NEWS b/NEWS index 1251f1b..dab747d 100644 --- a/NEWS +++ b/NEWS @@ -3,6 +3,7 @@ mrds 2.1.14 * updated initialvalues calculation for hazard-rate -- now uses Beavers & Ramsay method to scale parameters for hazard-rate * automatic parameter rescaling for covariate models when covariates are poorly scaled. Now default for nlminb method + * minor speed-up to logistic code when distance is a covariate mrds 2.1.13 diff --git a/R/g0.R b/R/g0.R index 8ed9bb7..6f0be73 100644 --- a/R/g0.R +++ b/R/g0.R @@ -6,5 +6,5 @@ #' @return vector of p(0) values #' @author Jeff Laake g0 <- function(beta, z){ - exp(as.matrix(z) %*% beta)/(1 + exp(as.matrix(z) %*% beta)) + exp(z %*% beta)/(1 + exp(z %*% beta)) } diff --git a/R/integratelogisticdup.R b/R/integratelogisticdup.R index f04d955..a1aec4e 100644 --- a/R/integratelogisticdup.R +++ b/R/integratelogisticdup.R @@ -1,7 +1,27 @@ integratelogisticdup <- function(x1, x2, models, beta, lower=0, width, point){ # numerical integral of product of logistic detection functions - integrate(logisticdupbyx, lower=lower, upper=width, - subdivisions=10, rel.tol=0.01, abs.tol=0.01, - x1=x1, x2=x2, models=models, beta=beta, point=point)$value + # computation speed-up when there is distance in the formula + # but only if there are no interactions with distance + if(sum(grepl("distance",names(beta)))==1){ + + # set parameter to be zero for distance + beta_distance <- beta[grepl("distance",names(beta))] + beta[grepl("distance",names(beta))] <- 0 + + # calculate the rest of the linear predictor + x1 <- setcov(x1,models$g0model)%*%beta + x2 <- setcov(x2,models$g0model)%*%beta + + # do some integration + integrate(logisticdupbyx_fast, lower=lower, upper=width, + subdivisions=10, rel.tol=0.01, abs.tol=0.01, + x1=x1, x2=x2, models=models, beta=beta, point=point, + beta_distance=beta_distance)$value + }else{ + # Otherwise just go ahead and do the numerical integration + integrate(logisticdupbyx, lower=lower, upper=width, + subdivisions=10, rel.tol=0.01, abs.tol=0.01, + x1=x1, x2=x2, models=models, beta=beta, point=point)$value + } } diff --git a/R/logisticdupbyx.R b/R/logisticdupbyx.R index 2b480fc..42ab4c3 100644 --- a/R/logisticdupbyx.R +++ b/R/logisticdupbyx.R @@ -14,20 +14,24 @@ #' @return vector of probabilities #' @author Jeff Laake logisticdupbyx <- function(distance, x1, x2, models, beta, point){ + + # avoid using g0 which calls exp and matrix multiplication twice + ologit <- function(p) p/(1+p) + # Functions used: g0, setcov xlist <- as.list(x1) xlist$distance <- distance xmat <- expand.grid(xlist) - gx1 <- g0(beta, setcov(xmat, models$g0model)) + gx1 <- ologit(exp(setcov(xmat, models$g0model) %*% beta)) xlist <- as.list(x2) xlist$distance <- distance xmat <- expand.grid(xlist) if(!point){ - return(gx1 * g0(beta, setcov(xmat, models$g0model))) + return(gx1 * ologit(exp(setcov(xmat, models$g0model) %*% beta))) }else{ - return(gx1 * g0(beta, setcov(xmat, models$g0model))*2*distance) + return(gx1 * ologit(exp(setcov(xmat, models$g0model) %*% beta))*2*distance) } } diff --git a/R/logisticdupbyx_fast.R b/R/logisticdupbyx_fast.R new file mode 100644 index 0000000..36738fb --- /dev/null +++ b/R/logisticdupbyx_fast.R @@ -0,0 +1,24 @@ +#' Logistic for duplicates as a function of covariates (fast) +#' +#' As \code{\link{logisticdupbyx}}, but faster when distance is a covariate (but no interactions with distance occur. +#' +#' @inheritParams logisticdupbyx +#' @param beta_distance parameter for distance +#' @param x1 linear predictor for 1, without distance +#' @param x2 linear predictor for 2, without distance +#' @author David L Miller +logisticdupbyx_fast <- function(distance, x1, x2, models, beta, point, beta_distance){ + + # function to calculate p/(1+p) + ologit <- function(p) p/(1+p) + + # first part of the function + gx1 <- ologit(exp(x1 + distance*beta_distance)) + + # calculate second and return + if(!point){ + return(gx1 * ologit(exp(x2 + distance*beta_distance))) + }else{ + return(gx1 * ologit(exp(x2 + distance*beta_distance))*2*distance) + } +} diff --git a/R/predict.io.fi.R b/R/predict.io.fi.R index 986ee00..7d82a97 100644 --- a/R/predict.io.fi.R +++ b/R/predict.io.fi.R @@ -58,19 +58,19 @@ predict.io.fi <- function(object,newdata=NULL,compute=FALSE, int.range=NULL, # now int.range is a vector with lower and upper bounds if(is.null(int.range)){ - pdot.list <- pdot.dsr.integrate.logistic(width,width, model$mr$coef, - newdata,integral.numeric, FALSE, models,GAM, point=point) + pdot.list <- pdot.dsr.integrate.logistic(width, width, model$mr$coef, + newdata, integral.numeric, FALSE, models,GAM, point=point) }else{ pdot.list <- pdot.dsr.integrate.logistic(int.range,width, model$mr$coef, - newdata,integral.numeric, FALSE, models,GAM, point=point) + newdata, integral.numeric, FALSE, models,GAM, point=point) } # if there is left truncation, take that off the integral if(left !=0){ pdot.list$pdot <- pdot.list$pdot - pdot.dsr.integrate.logistic(left, width, model$mr$coef, - newdata,integral.numeric, FALSE, models,GAM, - point=point)$pdot + newdata, integral.numeric, FALSE, models, + GAM, point=point)$pdot } fitted <- pdot.list$pdot diff --git a/man/logisticdupbyx_fast.Rd b/man/logisticdupbyx_fast.Rd new file mode 100644 index 0000000..75e9e58 --- /dev/null +++ b/man/logisticdupbyx_fast.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/logisticdupbyx_fast.R +\name{logisticdupbyx_fast} +\alias{logisticdupbyx_fast} +\title{Logistic for duplicates as a function of covariates (fast)} +\usage{ +logisticdupbyx_fast(distance, x1, x2, models, beta, point, beta_distance) +} +\arguments{ +\item{distance}{vector of distance values} + +\item{x1}{linear predictor for 1, without distance} + +\item{x2}{linear predictor for 2, without distance} + +\item{models}{model list} + +\item{beta}{logistic parameters} + +\item{point}{\code{TRUE} for point transect data} + +\item{beta_distance}{parameter for distance} +} +\description{ +As \code{\link{logisticdupbyx}}, but faster when distance is a covariate (but no interactions with distance occur. +} +\author{ +David L Miller +} +