R/rNtcondN0.R

In MobilePhoneESSnetBigData/pestim: Population Estimations Using Mobile Phone Data

#' @title Conditioned generation of random deviates of the posterior distribution of population counts.
#'
#' @description Generate random deviates of the posterior distribution of the number of individuals at an
#' arbitrary time instant conditioned upon the initial population.
#'
#' @param n number of values to generate
#'
#' @param N0 initial population in each cell
#'
#' @param nMNOmat transition matrix with the number of individuals displaced from cell to cell
#' detected by the Mobile Network Operator
#'
#' @param distNames character vector with the names of the prior distributions for each cell
#'
#' @param variation list of lists whose components are parameters providing a measure of variation
#' of each prior distribution
#'
#' @return Return a matrix with as many columns as cells and \code{n} rows with the generated values
#'
#' @details The function generates the probabilities according to a Dirichlet distribution with
#' parameters generated by \code{\link{alphaPrior}}. These parameters are generated with
#' distributions whose names are taken from the input parameter \code{distNames} and construct the
#' corresponding prior distribution for each cell \eqn{j} with mode at \eqn{u_{j}^{*}=N_{j}}, where
#' \eqn{N_{j}} is taken from the sum of rows of \code{nMNOmat}. Next the rest of parameters of the
#' distribution are computed according to the dispersion parameters specified in \code{variation}.
#'
#'  As accepted distribution names, currently the user can specify \code{unif}, \code{degen},
#'  \code{triang}, and \code{gamma}.
#'
#'  The dispersion parameters recognised so far are the coefficients of variation only (standard
#'  deviation divided by the mean of the distribution). These dispersion parameters must be
#'  specified by a named component \code{cv} with a numeric value in \eqn{[0, 1]}.
#'
#'  For each distribution the parameters are computed as follows:
#'
#'  \itemize{
#'
#'  \item \code{unif}: This is the uniform distribution with parameters \code{xMax} and \code{xMin}.
#'   Both parameters are computed by \eqn{u_{j}^{*}\cdot(1\pm\sqrt{3}\textrm{cv})}, respectively, in
#'   each cell \eqn{j}.
#'
#'  \item \code{degen}: This is the degenerate distribution with parameter \code{X0} taken as
#'  \eqn{u_{j}^{*}} in each cell \eqn{j}.
#'
#'  \item \code{triang}: This is the triangular distribution \code{\link{triang}} with parameters
#'  \code{xMax}, \code{xMin}, and \code{xMode}. The latter is taken directly from \code{nMNOfrom}.
#'  The distribution is assumed to be symmetrical so that the two former parameters are computed by
#'  \eqn{u_{j}^{*}\cdot(1\pm\sqrt{3}\textrm{cv})}, respectively, in each cell \eqn{j}.
#'
#'  \item \code{gamma}: This is the gamma distribution with parameters \code{shape} and \code{scale}.
#'  The former is computed as \eqn{\frac{1}{\textrm{cv}^2}} and the latter as
#'  \eqn{frac{u_{j}^{*}}{\textrm{scale} - 1}}.
#'
#'  }
#'
#' @examples
#' N0 <- c(93, 123, 130)
#' nMNOmat <- rbind(c(10, 3, 4), c(5, 21, 3), c(3, 9, 18))
#' distNames <- rep('unif', 3)
#' variation <- rep(list(list(cv = 0.20)), 3)
#' rNtcondN0(1e3, N0, nMNOmat, distNames, variation)
#'
#' @include rmatProb.R
#'
#' @export
#'

rNtcondN0 <- function(n, N0, nMNOmat,  distNames, variation){

  nCells <- length(N0)
  if (!is.matrix(nMNOmat)) stop('nMNOmat must be a square matrix.')
  if (dim(nMNOmat)[1] != nCells) stop('The number of rows of nMNOmat must coincide with the number of cells.')
  if (dim(nMNOmat)[2] != nCells) stop('The number of columns of nMNOmat must coincide with the number of cells.')
  if (length(distNames) != nCells) stop('The length of distNames must coincide with the number of cells.')
  if (length(variation) != nCells) stop('The length of variation must coincide with the number of cells.')

  matProb <- rmatProb(n, nMNOmat, distNames, variation)

  output <- lapply(matProb, function(matP){

    outLocal <- N0 %*% matP
    outLocal <- round(outLocal)
    return(outLocal)

  })

  output <- Reduce(rbind, output)
  dimnames(output) <- NULL
  return(output)

}