R/util_funcs.R

In gilesjohnr/genpatch: Tools for individual based epidemic simulations with pathogen evolution

##' Find parameters of Gamma distribution
##'
##' A function that finds the \code{shape} and \code{rate} paramters required by the Gamma distribution given the observed mean \code{mu} and standard deviation \code{sigma}. Parameters are determined using a two-dimensional Nelder-Mead optimization algorithm.
##'
##' @param mu the desired mean of the Gamma distribution
##' @param sigma the desired standard deviation of the Gamma distribution
##'
##' @return a named numeric vector giving the \code{shape} and \code{rate} parameters of the Gamma distribution
##'
##' @author John Giles
##'
##' @example R/examples/gamma_params.R
##'

gamma.params <- function(mu, sigma) {
     suppressWarnings(
          params <- optim(par=c(mu*2, 2),
                          fn=function(x) abs(mu - x[1]/x[2]) + abs(sigma - sqrt(x[1]/(x[2]^2))),
                          method='Nelder-Mead')$par
     )
     names(params) <- c('shape', 'rate')
     return(params)
}

##' Number of secondary cases given range of R values
##'
##' A function that generates random values of the basic reproductive number \emph{R} when given an estimated an estimated confidence interval for \emph{R}. The function draws a point estimate for \emph{R} within the given range and then returns a randomly generated number from a Poisson distribution with this point estimate as the mean.
##'
##' @param n number of random numbers to generate
##' @param R_low the low estimate of \emph{R}. If the \code{R_high} parameter is \code{NULL}, this parameter is taken to be the point estimate of \emph{R}
##' @param R_high the high estimate of \emph{R}
##'
##' @return an \code{n} length numeric vector giving randomly generated \emph{R} values
##'
##' @author John Giles
##'
##' @example R/examples/n_sec_cases.R
##'

n.sec.cases <- function(n=1, R_low, R_high=NULL) {
     if (is.null(R_high) == TRUE) {
          out <- rpois(n, round(R_low))
     } else {
          R <- round(runif(n=n, R_low, R_high))
          out <- rpois(n, R)
     }
     return(out)
}