clean R CMD CHECK as-cran

MultiBD.Rproj

@@ -12,6 +12,9 @@ Encoding: UTF-8
 RnwWeave: knitr
 LaTeX: pdfLaTeX
 
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
 BuildType: Package
 PackageInstallArgs: --no-multiarch --with-keep.source --no-build-vignettes
 PackageBuildArgs: --no-build-vignettes

R/SIR_prob.R

@@ -35,8 +35,8 @@
 #'                log(
 #'                  SIR_prob(  # Compute the forward transition probability matrix
 #'                    t  = data$time[k + 1] - data$time[k], # Time increment
-#'                    alpha = alpha, beta = beta,
-#'                    S0 = data$S[k], I0 = data$I[k],       # From: R(t_k), I(t_k)                                      
+#'                    alpha = alpha, beta = beta, 
+#'                    S0 = data$S[k], I0 = data$I[k],       # From: R(t_k), I(t_k)
 #'                    nSI = data$S[k] - data$S[k + 1], nIR = data$R[k + 1] - data$R[k],
 #'                    computeMode = 4, nblocks = 80         # Compute using 4 threads
 #'                  )[data$S[k] - data$S[k + 1] + 1, 

R/bbd_prob.R

@@ -43,9 +43,9 @@
 #'   sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k
 #'              function(k) {
 #'                log(
-#'                  bbd_prob(  # Compute the transition probability matrix
+#'                  bbd_prob( # Compute the transition probability matrix
 #'                    t  = data$time[k + 1] - data$time[k], # Time increment
-#'                    a0 = data$R[k], b0 = data$I[k],       # From: R(t_k), I(t_k)                                      
+#'                    a0 = data$R[k], b0 = data$I[k],       # From: R(t_k), I(t_k)
 #'                    brates1, brates2, drates2, trans21,
 #'                    A = data$R[k + 1], B = data$R[k + 1] + data$I[k] - data$R[k],
 #'                    computeMode = 4, nblocks = 80         # Compute using 4 threads
@@ -112,7 +112,7 @@ bbd_prob <- function(t, a0, b0, lambda1, lambda2, mu2, gamma, A, B,
     y = l1[a-a0+1,b+1] + l2[a-a0+1,b+1] + m2[a-a0+1,b+1] + g[a-a0+1,b+1]
     return(y)
   }
-  y = matrix(mapply(yf, grid[,1], grid[,2]), nrow=B+1+maxdepth, byrow=TRUE)	
+  y = matrix(mapply(yf, grid[,1], grid[,2]), nrow=B+1+maxdepth, byrow=TRUE) 
 
   #############################
   ### Call C function via Rcpp

R/dbd_prob.R

@@ -4,7 +4,7 @@
 
 #' Transition probabilities of a death/birth-death process
 #'
-#' Computes the transition pobabilities of a death/birth-death process 
+#' Computes the transition pobabilities of a death/birth-death process
 #' using the continued fraction representation of its Laplace transform
 #' @param t time
 #' @param a0 total number of type 1 particles at \code{t = 0}
@@ -22,101 +22,119 @@
 #' @param maxdepth maximum number of iterations for Lentz algorithm
 #' @return a matrix of the transition probabilities
 #' @references Ho LST et al. 2016. "Birth(death)/birth-death processes and their computable transition probabilities with statistical applications". In review.
-#' 
+#'
 #' @examples
 #' \dontrun{
 #' data(Eyam)
-#' 
+#'
 #' loglik_sir <- function(param, data) {
 #'   alpha <- exp(param[1]) # Rates must be non-negative
 #'   beta  <- exp(param[2])
-#'   
+#'
 #'   # Set-up SIR model
 #'   drates1 <- function(a, b) { 0 }
 #'   brates2 <- function(a, b) { 0 }
 #'   drates2 <- function(a, b) { alpha * b     }
 #'   trans12 <- function(a, b) { beta  * a * b }
-#'   
+#'
 #'   sum(sapply(1:(nrow(data) - 1), # Sum across all time steps k
 #'              function(k) {
 #'                log(
 #'                  dbd_prob(  # Compute the transition probability matrix
 #'                    t  = data$time[k + 1] - data$time[k], # Time increment
-#'                    a0 = data$S[k], b0 = data$I[k],       # From: S(t_k), I(t_k)                                      
+#'                    a0 = data$S[k], b0 = data$I[k],       # From: S(t_k), I(t_k)
 #'                    drates1, brates2, drates2, trans12,
 #'                    a = data$S[k + 1], B = data$S[k] + data$I[k] - data$S[k + 1],
 #'                    computeMode = 4, nblocks = 80         # Compute using 4 threads
 #'                  )[1, data$I[k + 1] + 1]                 # To: S(t_(k+1)), I(t_(k+1))
 #'                )
 #'              }))
+#'   }
+#'
+#'   loglik_sir(log(c(3.204, 0.019)), Eyam) # Evaluate at mode
 #' }
-#' 
-#' loglik_sir(log(c(3.204, 0.019)), Eyam) # Evaluate at mode
-#' }
-#' 
+#'
 #' # Birth-death-shift model for transposable elements
-#' 
+#'
 #' lam = 0.0188; mu = 0.0147; v = 0.00268; # birth, death, shift rates
-#' 
+#'
 #' drates1 <- function(a, b) { mu * a }
 #' brates2 <- function(a, b) { lam * (a + b) }
 #' drates2 <- function(a, b) { mu * b }
 #' trans12 <- function(a, b) { v * a }
-#' 
+#'
 #' # Get transition probabilities
-#' p <- dbd_prob(t = 1, a0 = 10, b0 = 0, 
-#'               drates1, brates2, drates2, trans12, 
-#'               a = 0, B = 50) 
-#' 
+#' p <- dbd_prob(t = 1, a0 = 10, b0 = 0,
+#'               drates1, brates2, drates2, trans12,
+#'               a = 0, B = 50)
+#'
 #' @export
 dbd_prob <- function(t, a0, b0, mu1, lambda2, mu2, gamma, a=0, B,
                     nblocks=256, tol=1e-12, computeMode=0, nThreads=4,
                     maxdepth=400) {
-  
+
   ###################
   ### Input checking
   ###################
-  
-  ## a>=0, a<=a0, B >=a0+b0-a 
+
+  ## a>=0, a<=a0, B >=a0+b0-a
   if(a < 0) stop("a cannot be negative.")
   if (a > a0) stop("a cannot be bigger than a0.")
   if (B < a0+b0-a) stop("B is too small.")
-  
+
   ###########################################################
   ### store rate values in matrices
   ### rates are transformed to fit birth/birth-death process
   ###########################################################
-  
+
   l1 <- function(u,v){
     if (v > B) return(0)
     return(mu1(a0-u,B-v))
   }
-  
+
   l2 <- function(u,v){
     if (v > B) return(0)
     return(mu2(a0-u,B-v))
   }
-  
+
   m2 <- function(u,v){
     if (v > B) return(0)
     return(lambda2(a0-u,B-v))
   }
-  
+
   g <- function(u,v){
     if (v > B) return(0)
     return(gamma(a0-u,B-v))
   }
-  
+
   ###########################
   ### Call bbd_prob function
   ###########################
-  
+
   res = matrix(0, nrow=a0-a+1, ncol=B+1)
   res[(a0-a+1):1,(B+1):1] = bbd_prob(t, 0, B-b0, l1, l2, m2, g, A=a0-a, B,
                                      nblocks, tol, computeMode, nThreads, maxdepth)
-  
+
   colnames(res) = 0:B
   rownames(res) = a:a0
-  
+
   return(res)
 }
+
+# function() {
+#
+#   sum(
+#     sapply(
+#       1:(nrow(data) - 1), # Sum across all time steps k
+#       function(k) {
+#         log(
+#           dbd_prob( # Compute the transition probability matrix
+#             t = data$time[k + 1] - data$time[k], # Time increment
+#             a0 = data$S[k], b0 = data$I[k], # From: S(t_k), I(t_k)
+#             drates1, brates2, drates2, trans12,
+#                  a = data$S[k + 1],
+#                  B = data$S[k] + data$I[k] - data$S[k + 1],
+#                  computeMode = 4, nblocks = 80 # Compute using 4 threads
+#             )[1, data$I[k + 1] + 1] # To: S(t_(k+1)), I(t_(k+1))
+#           )}))
+# }